0
# High-Performance C Extensions
1

2
Low-level C functions providing optimized implementations of core drizzling algorithms for maximum performance with large astronomical datasets. These functions are typically called internally by the Python interface but can be used directly for specialized applications.
3

4
```python { .api }
5
from typing import Optional, Tuple, Any
6
import numpy as np
7
```
8

9
## Capabilities
10

11
### Core Drizzling Function
12

13
The primary C implementation of the drizzling algorithm, handling pixel redistribution with various kernels and comprehensive parameter control.
14

15
```python { .api }
16
def tdriz(
17
    input: np.ndarray,
18
    weights: np.ndarray,
19
    pixmap: np.ndarray,
20
    output: np.ndarray,
21
    counts: np.ndarray,
22
    context: Optional[np.ndarray],
23
    uniqid: int,
24
    xmin: int,
25
    xmax: int,
26
    ymin: int,
27
    ymax: int,
28
    scale: float,
29
    pixfrac: float,
30
    kernel: str,
31
    in_units: str,
32
    expscale: float,
33
    wtscale: float,
34
    fillstr: str
35
) -> Tuple[str, float, float]:
36
    """
37
    Core drizzling algorithm implementation in C for maximum performance.
38

39
    This function directly implements the drizzle algorithm at the C level,
40
    providing the computational engine for the Python Drizzle class.
41

42
    Parameters:
43
    - input: Input image array (float32, 2D)
44
    - weights: Pixel weighting array (float32, 2D, same shape as input)
45
    - pixmap: Coordinate mapping array (float64, 3D with shape (Ny, Nx, 2))
46
    - output: Output image array (float32, 2D, modified in-place)
47
    - counts: Output weight array (float32, 2D, same shape as output, modified in-place)
48
    - context: Context bit-field array (int32, 2D, same shape as output, modified in-place)
49
               Can be None to disable context tracking
50
    - uniqid: Unique identifier for this input image (1-based, used for context bit field)
51
    - xmin: Minimum X coordinate of bounding box in input image (0-based)
52
    - xmax: Maximum X coordinate of bounding box in input image (0-based)
53
    - ymin: Minimum Y coordinate of bounding box in input image (0-based)
54
    - ymax: Maximum Y coordinate of bounding box in input image (0-based)
55
    - scale: Pixel scale factor for flux scaling
56
    - pixfrac: Pixel fraction (0-1) controlling flux distribution area
57
    - kernel: Drizzling kernel name ("square", "gaussian", "point", "turbo", "lanczos2", "lanczos3")
58
    - in_units: Input units ("cps" for counts per second, "counts" for total counts)
59
    - expscale: Exposure scaling factor
60
    - wtscale: Weight scaling factor
61
    - fillstr: Fill value string for unsampled pixels ("INDEF", "NaN", or numeric string)
62

63
    Returns:
64
    Tuple of (version_string, nmiss, nskip):
65
    - version_string: C library version identifier
66
    - nmiss: Number of input pixels that did not contribute to output (float)
67
    - nskip: Number of input lines that were skipped (float)
68

69
    Notes:
70
    - All arrays must be C-contiguous for optimal performance
71
    - Output, counts, and context arrays are modified in-place
72
    - Context tracking uses bit fields; bit (uniqid-1) is set for contributing pixels
73
    - Kernel choice affects flux distribution and conservation properties
74
    - Performance-critical function optimized for large astronomical datasets
75
    """
76
```
77

78
### Core Blotting Function
79

80
C implementation of the blotting (inverse drizzling) algorithm using various interpolation methods.
81

82
```python { .api }
83
def tblot(
84
    image: np.ndarray,
85
    pixmap: np.ndarray,
86
    output: np.ndarray,
87
    scale: float,
88
    kscale: float,
89
    interp: str,
90
    exptime: float,
91
    misval: float,
92
    sinscl: float
93
) -> None:
94
    """
95
    Core blotting algorithm implementation in C.
96

97
    Performs inverse drizzling by resampling input image onto output grid
98
    using interpolation methods. Optimized for performance with large images.
99

100
    Parameters:
101
    - image: Input image to be blotted (float32, 2D)
102
    - pixmap: Coordinate mapping array (float64, 3D with shape (Ny, Nx, 2))
103
              Maps from output pixel coordinates to input pixel coordinates
104
    - output: Output image array (float32, 2D, modified in-place)
105
              Must be pre-allocated with correct size
106
    - scale: Pixel ratio scale factor
107
    - kscale: Kernel scale factor (typically 1.0)
108
    - interp: Interpolation method string:
109
              "nearest" - nearest neighbor
110
              "linear" - bilinear interpolation
111
              "poly3" - cubic polynomial
112
              "poly5" - quintic polynomial
113
              "sinc" - sinc interpolation
114
              "lan3" - 3rd order Lanczos
115
              "lan5" - 5th order Lanczos
116
    - exptime: Exposure time for flux scaling
117
    - misval: Value to use for missing/invalid pixels
118
    - sinscl: Scaling factor for sinc interpolation (used when interp="sinc")
119

120
    Returns:
121
    None (modifies output array in-place)
122

123
    Notes:
124
    - Output array must be pre-allocated and zeroed
125
    - Not flux-conserving (unlike tdriz)
126
    - Best used with well-sampled input images
127
    - pixmap maps from output coordinates to input coordinates (inverse of tdriz)
128
    - Performance optimized for large astronomical image processing
129
    """
130
```
131

132
### Testing Function
133

134
C-level testing function for validating drizzling operations and performance benchmarking.
135

136
```python { .api }
137
def test_cdrizzle(
138
    data: np.ndarray,
139
    weights: np.ndarray,
140
    pixmap: np.ndarray,
141
    output_data: np.ndarray,
142
    output_counts: np.ndarray,
143
    output_context: np.ndarray
144
) -> None:
145
    """
146
    Test function for C drizzling implementation.
147

148
    Runs C-level unit tests for the cdrizzle implementation, providing validation
149
    and performance testing of the core drizzling algorithms.
150

151
    Parameters:
152
    - data: Input test data array (float32, 2D)
153
    - weights: Input weight array (float32, 2D, same shape as data)
154
    - pixmap: Pixel mapping array (float64, 3D with shape (height, width, 2))
155
    - output_data: Output data array (float32, 2D, modified in-place)
156
    - output_counts: Output counts array (float32, 2D, same shape as output_data, modified in-place)
157
    - output_context: Context bit-field array (int32, 2D, same shape as output_data, modified in-place)
158

159
    Returns:
160
    None
161

162
    Notes:
163
    - Primarily for internal testing and validation
164
    - Calls internal C function utest_cdrizzle() for validation
165
    - All arrays must have correct dtypes and be C-contiguous
166
    - Not typically used in production workflows
167
    """
168
```
169

170
### Pixel Map Inversion
171

172
Utility function for inverting coordinate transformation arrays using iterative search algorithms.
173

174
```python { .api }
175
def invert_pixmap(
176
    pixmap: np.ndarray,
177
    xyout: np.ndarray,
178
    bbox: Optional[np.ndarray]
179
) -> np.ndarray:
180
    """
181
    Invert pixel mapping transformation using golden ratio search.
182

183
    Takes output coordinates and finds corresponding input coordinates by
184
    iteratively searching through the pixmap transformation array.
185

186
    Parameters:
187
    - pixmap: 3D coordinate mapping array (float64, shape (height, width, 2))
188
              Contains forward transformation from input to output coordinates
189
    - xyout: 1D array (float64, 2 elements) containing output coordinates [x_out, y_out]
190
    - bbox: 2D bounding box array (float64, shape (2, 2)) as [[xmin, xmax], [ymin, ymax]]
191
            or None to use full pixmap extent
192

193
    Returns:
194
    1D array (float64, 2 elements) containing inverted input coordinates [x_in, y_in]
195

196
    Notes:
197
    - Used internally for coordinate system conversions
198
    - Helps convert between drizzling and blotting coordinate mappings
199
    - Uses golden ratio search algorithm for iterative inversion
200
    - Essential for computing reverse transformations
201
    """
202
```
203

204
### Polygon Clipping
205

206
Geometric utility for clipping polygons using the Sutherland-Hodgman algorithm.
207

208
```python { .api }
209
def clip_polygon(
210
    p: np.ndarray,
211
    q: np.ndarray
212
) -> List[Tuple[float, float]]:
213
    """
214
    Clip polygon p to the window defined by polygon q.
215

216
    Implements the Sutherland-Hodgman polygon clipping algorithm for geometric
217
    operations on pixel regions and image boundaries.
218

219
    Parameters:
220
    - p: 2D array (float64, shape (n_vertices, 2)) representing first polygon vertices
221
         Each row contains [x, y] coordinates of a vertex
222
    - q: 2D array (float64, shape (m_vertices, 2)) representing clipping window polygon
223
         Each row contains [x, y] coordinates of a clipping window vertex
224

225
    Returns:
226
    List of tuples, each containing (x, y) coordinates as floats representing
227
    vertices of the clipped polygon
228

229
    Notes:
230
    - Both input polygons are automatically oriented counter-clockwise before clipping
231
    - Uses Sutherland-Hodgman clipping algorithm for robust polygon intersection
232
    - Primarily used for geometric operations on image pixel regions
233
    - Returns empty list if polygons do not intersect
234
    """
235
```
236

237
## Usage Examples
238

239
### Direct C Function Usage
240

241
```python
242
import numpy as np
243
from drizzle import cdrizzle
244

245
# Prepare input arrays
246
input_data = np.random.random((100, 100)).astype(np.float32)
247
weights = np.ones((100, 100), dtype=np.float32)
248

249
# Create pixel mapping (identity transformation for example)
250
y, x = np.indices((100, 100), dtype=np.float64)
251
pixmap = np.dstack([x, y])
252

253
# Prepare output arrays
254
output_shape = (120, 120)
255
output_data = np.zeros(output_shape, dtype=np.float32)
256
output_counts = np.zeros(output_shape, dtype=np.float32)
257
context = np.zeros(output_shape, dtype=np.int32)
258

259
# Call C drizzling function directly
260
version, nmiss, nskip = cdrizzle.tdriz(
261
    input=input_data,
262
    weights=weights,
263
    pixmap=pixmap,
264
    output=output_data,
265
    counts=output_counts,
266
    context=context,
267
    uniqid=1,  # First image
268
    xmin=0,
269
    xmax=99,
270
    ymin=0,
271
    ymax=99,
272
    scale=1.0,
273
    pixfrac=1.0,
274
    kernel="square",
275
    in_units="cps",
276
    expscale=1.0,
277
    wtscale=1.0,
278
    fillstr="INDEF"
279
)
280

281
print(f"C library version: {version}")
282
print(f"Missed pixels: {nmiss}, Skipped lines: {nskip}")
283
print(f"Output data range: {output_data.min():.3f} to {output_data.max():.3f}")
284
```
285

286
### Direct Blotting Example
287

288
```python
289
import numpy as np
290
from drizzle import cdrizzle
291

292
# Create well-sampled input image (e.g., drizzled result)
293
input_image = np.random.random((150, 150)).astype(np.float32)
294

295
# Create coordinate mapping from output to input
296
output_shape = (100, 100)
297
y, x = np.indices(output_shape, dtype=np.float64)
298

299
# Scale coordinates to map to input image
300
x_scaled = x * 1.5  # Scale factor
301
y_scaled = y * 1.5
302
pixmap = np.dstack([x_scaled, y_scaled])
303

304
# Prepare output array
305
blotted_image = np.zeros(output_shape, dtype=np.float32)
306

307
# Call C blotting function directly
308
cdrizzle.tblot(
309
    image=input_image,
310
    pixmap=pixmap,
311
    output=blotted_image,
312
    scale=1.5,
313
    kscale=1.0,
314
    interp="poly5",
315
    exptime=1.0,
316
    misval=0.0,
317
    sinscl=1.0
318
)
319

320
print(f"Blotted image range: {blotted_image.min():.3f} to {blotted_image.max():.3f}")
321
```
322

323
### Performance Testing
324

325
```python
326
import numpy as np
327
import time
328
from drizzle import cdrizzle
329

330
# Create large test arrays for performance testing
331
size = 2048
332
input_data = np.random.random((size, size)).astype(np.float32)
333
weights = np.ones_like(input_data)
334

335
# Create pixel mapping
336
y, x = np.indices((size, size), dtype=np.float64)
337
pixmap = np.dstack([x, y])
338

339
# Prepare outputs
340
output_data = np.zeros_like(input_data)
341
output_counts = np.zeros_like(input_data)
342

343
# Time the C function
344
start_time = time.time()
345

346
version, nmiss, nskip = cdrizzle.tdriz(
347
    input=input_data,
348
    weights=weights,
349
    pixmap=pixmap,
350
    output=output_data,
351
    counts=output_counts,
352
    context=None,  # Disable context for speed
353
    uniqid=1,
354
    xmin=0, xmax=size-1,
355
    ymin=0, ymax=size-1,
356
    scale=1.0,
357
    pixfrac=1.0,
358
    kernel="square",
359
    in_units="cps",
360
    expscale=1.0,
361
    wtscale=1.0,
362
    fillstr="0"
363
)
364

365
elapsed = time.time() - start_time
366
pixels_per_sec = (size * size) / elapsed
367

368
print(f"Processed {size}x{size} image in {elapsed:.2f} seconds")
369
print(f"Performance: {pixels_per_sec/1e6:.2f} Mpixels/sec")
370
```
371

372
### Kernel Comparison
373

374
```python
375
import numpy as np
376
from drizzle import cdrizzle
377

378
def test_kernel(kernel_name, input_data, pixmap):
379
    """Test different drizzling kernels."""
380
    output_shape = (150, 150)
381
    output_data = np.zeros(output_shape, dtype=np.float32)
382
    output_counts = np.zeros(output_shape, dtype=np.float32)
383
    weights = np.ones_like(input_data)
384

385
    version, nmiss, nskip = cdrizzle.tdriz(
386
        input=input_data,
387
        weights=weights,
388
        pixmap=pixmap,
389
        output=output_data,
390
        counts=output_counts,
391
        context=None,
392
        uniqid=1,
393
        xmin=0, xmax=input_data.shape[1]-1,
394
        ymin=0, ymax=input_data.shape[0]-1,
395
        scale=1.0,
396
        pixfrac=1.0,
397
        kernel=kernel_name,
398
        in_units="cps",
399
        expscale=1.0,
400
        wtscale=1.0,
401
        fillstr="0"
402
    )
403

404
    return output_data, nmiss, nskip
405

406
# Test data
407
input_data = np.random.random((100, 100)).astype(np.float32)
408
y, x = np.indices((100, 100), dtype=np.float64)
409
pixmap = np.dstack([x * 1.2, y * 1.2])  # Slight magnification
410

411
# Test different kernels
412
kernels = ["square", "gaussian", "point", "turbo", "lanczos2", "lanczos3"]
413

414
for kernel in kernels:
415
    try:
416
        result, nmiss, nskip = test_kernel(kernel, input_data, pixmap)
417
        print(f"Kernel '{kernel}': output range [{result.min():.3f}, {result.max():.3f}], "
418
              f"missed={nmiss}, skipped={nskip}")
419
    except Exception as e:
420
        print(f"Kernel '{kernel}' failed: {e}")
421
```
422

423
## Error Handling
424

425
```python
426
import numpy as np
427
from drizzle import cdrizzle
428

429
# Handle array type mismatches
430
try:
431
    bad_input = np.ones((50, 50), dtype=np.int32)  # Wrong dtype
432
    weights = np.ones((50, 50), dtype=np.float32)
433
    y, x = np.indices((50, 50), dtype=np.float64)
434
    pixmap = np.dstack([x, y])
435
    output = np.zeros((50, 50), dtype=np.float32)
436
    counts = np.zeros((50, 50), dtype=np.float32)
437

438
    cdrizzle.tdriz(
439
        input=bad_input,  # int32 instead of float32
440
        weights=weights,
441
        pixmap=pixmap,
442
        output=output,
443
        counts=counts,
444
        context=None,
445
        uniqid=1,
446
        xmin=0, xmax=49, ymin=0, ymax=49,
447
        scale=1.0, pixfrac=1.0,
448
        kernel="square", in_units="cps",
449
        expscale=1.0, wtscale=1.0, fillstr="0"
450
    )
451
except (TypeError, ValueError) as e:
452
    print(f"Array type error: {e}")
453

454
# Handle invalid kernel names
455
try:
456
    input_data = np.ones((50, 50), dtype=np.float32)
457
    # ... other arrays setup ...
458
    cdrizzle.tdriz(
459
        input=input_data, weights=weights, pixmap=pixmap,
460
        output=output, counts=counts, context=None,
461
        uniqid=1, xmin=0, xmax=49, ymin=0, ymax=49,
462
        scale=1.0, pixfrac=1.0,
463
        kernel="invalid_kernel",  # Bad kernel name
464
        in_units="cps", expscale=1.0, wtscale=1.0, fillstr="0"
465
    )
466
except Exception as e:
467
    print(f"Invalid kernel error: {e}")
468

469
# Handle array shape mismatches
470
try:
471
    input_data = np.ones((50, 50), dtype=np.float32)
472
    bad_weights = np.ones((40, 40), dtype=np.float32)  # Wrong shape
473
    # ... rest of setup ...
474
    cdrizzle.tdriz(
475
        input=input_data,
476
        weights=bad_weights,  # Shape mismatch
477
        # ... other parameters ...
478
    )
479
except Exception as e:
480
    print(f"Shape mismatch error: {e}")
481
```
482

483
## Notes
484

485
- These C functions provide the computational core of the drizzle package
486
- Direct usage requires careful array preparation and parameter validation
487
- The Python `Drizzle` class provides a safer, more convenient interface
488
- C functions modify output arrays in-place for memory efficiency
489
- All input arrays must be C-contiguous for optimal performance
490
- Context tracking uses bit fields and requires proper unique ID management
491
- Error handling is minimal at the C level; validation should be done in Python
492
- These functions are optimized for large astronomical datasets and can process gigapixel images efficiently