From bef3e1112ebcbe6f9952ea9dbf660577c184f11a Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Tue, 10 Dec 2024 03:25:15 -0500 Subject: [PATCH] improvements; fix unit tests; reviewer comments --- .coveragerc | 1 + drizzle/resample.py | 114 ++++++++++++++++++++++++--------- drizzle/tests/test_resample.py | 35 ++++++++-- src/cdrizzlebox.c | 79 +++++++++++------------ 4 files changed, 151 insertions(+), 78 deletions(-) diff --git a/.coveragerc b/.coveragerc index 644d554..ecee275 100644 --- a/.coveragerc +++ b/.coveragerc @@ -3,6 +3,7 @@ source = drizzle omit = drizzle/__init__* + drizzle/util.py drizzle/**/conftest.py drizzle/**/setup* drizzle/**/tests/* diff --git a/drizzle/resample.py b/drizzle/resample.py index 2dc85b1..26b78e7 100644 --- a/drizzle/resample.py +++ b/drizzle/resample.py @@ -315,17 +315,6 @@ def __init__(self, kernel="square", fillval=None, fillval2=None, out_img = np.asarray(out_img, dtype=np.float32) shapes.add(out_img.shape) - if out_img2 is not None: - if isinstance(out_img2, np.ndarray): - out_img2 = np.asarray(out_img2, dtype=np.float32) - shapes.add(out_img2.shape) - else: - out_img2 = [ - np.asarray(img, dtype=np.float32) for img in out_img2 - ] - for img in out_img2: - shapes.add(img.shape) - if out_wht is not None: out_wht = np.asarray(out_wht, dtype=np.float32) shapes.add(out_wht.shape) @@ -350,7 +339,6 @@ def __init__(self, kernel="square", fillval=None, fillval2=None, out_wht=out_wht, out_ctx=out_ctx, ) - self._alloc_output_arrays2_init(out_img2=out_img2) elif len(shapes) > 1: raise ValueError( @@ -363,11 +351,36 @@ def __init__(self, kernel="square", fillval=None, fillval2=None, self._out_wht = None self._out_ctx = None + if out_img2 is not None: + if self._out_shape is not None: + shapes.add(self._out_shape) + if isinstance(out_img2, np.ndarray): + out_img2 = np.asarray(out_img2, dtype=np.float32) + shapes.add(out_img2.shape) + else: + for img in out_img2: + if img is not None: + shapes.add(np.shape(img)) + if len(shapes) > 1: + raise ValueError( + "Inconsistent data shapes specified: 'out_shape' " + "and/or out_img, out_img2, out_wht, out_ctx have " + "different shapes." + ) + self._output_shapes = shapes + self._alloc_output_arrays2_init(out_img2=out_img2) + @property def fillval(self): - """Fill value for output pixels without contributions from input images.""" + """ Fill value for output pixels without contributions from input images. """ return self._fillval + @property + def fillval2(self): + """ Fill value for output pixels in ``out_img2`` without contributions + from input images. """ + return self._fillval2 + @property def kernel(self): """Resampling kernel.""" @@ -375,27 +388,27 @@ def kernel(self): @property def ctx_id(self): - """Context image "ID" (0-based ) of the next image to be resampled.""" + """ Context image "ID" (0-based ) of the next image to be resampled. """ return self._ctx_id @property def out_img(self): - """Output resampled image.""" + """ Output resampled image. """ return self._out_img @property def out_wht(self): - """Output weight image.""" + """ Output weight image. """ return self._out_wht @property def out_ctx(self): - """Output "context" image.""" + """ Output "context" image. """ return self._out_ctx @property def out_img2(self): - """Output resampled image(s) obtained with squared weights. + """ Output resampled image(s) obtained with squared weights. It is always a list of one or more 2D arrays. """ @@ -458,14 +471,35 @@ def _alloc_output_arrays2_init(self, out_img2=None): out_img2 = [out_img2] self._out_img2 = [] + + if (isinstance(self._fillval2, str) and + self._fillval2.strip().upper() == "INDEF"): + fv = np.nan + else: + fv = np.float32(self._fillval2) + for i2 in out_img2: if i2 is None: - arr = np.zeros(self._out_shape, dtype=np.float32) + if self._out_shape is None: + if len(self._output_shapes) == 1: + shape = list(self._output_shapes)[0] + else: + self._out_img2.append(None) + continue + else: + shape = self._out_shape + arr = np.full(shape, fill_value=fv, dtype=np.float32) else: arr = np.asarray(i2, dtype=np.float32) self._out_img2.append(arr) + del arr def _alloc_output_arrays2_add(self, ninputs2=None): + if (isinstance(self._fillval2, str) and + self._fillval2.strip().upper() == "INDEF"): + fv = np.nan + else: + fv = np.float32(self._fillval2) if self._out_img2 is None: if ninputs2 is None or ninputs2 < 1: # nothing to do @@ -476,13 +510,22 @@ def _alloc_output_arrays2_add(self, ninputs2=None): "and the number of inputs." ) self._out_img2 = [ - np.zeros(self._out_shape, dtype=np.float32) + np.full(self._out_shape, fill_value=fv, dtype=np.float32) for _ in range(ninputs2) ] else: nimg2 = len(self._out_img2) + # replace None values with arrays of _out_shape: + for k, img in enumerate(self._out_img2): + if img is None: + self._out_img2[k] = np.full( + self._out_shape, + fill_value=fv, + dtype=np.float32 + ) + if ((ninputs2 is not None and ninputs2 != nimg2) or (ninputs2 is None and nimg2 > 0)): raise ValueError( @@ -627,17 +670,26 @@ def add_image(self, data, exptime, pixmap, data2=None, scale=1.0, # set output image shape based on output coordinates from the pixmap. # if self._out_shape is None: - pmap_xmin = int(np.floor(np.nanmin(pixmap[:, :, 0]))) - pmap_xmax = int(np.ceil(np.nanmax(pixmap[:, :, 0]))) - pmap_ymin = int(np.floor(np.nanmin(pixmap[:, :, 1]))) - pmap_ymax = int(np.ceil(np.nanmax(pixmap[:, :, 1]))) - pixmap = pixmap.copy() - pixmap[:, :, 0] -= pmap_xmin - pixmap[:, :, 1] -= pmap_ymin - self._out_shape = ( - pmap_xmax - pmap_xmin + 1, - pmap_ymax - pmap_ymin + 1 - ) + nshapes = len(self._output_shapes) + if nshapes == 0: + pmap_xmin = int(np.floor(np.nanmin(pixmap[:, :, 0]))) + pmap_xmax = int(np.ceil(np.nanmax(pixmap[:, :, 0]))) + pmap_ymin = int(np.floor(np.nanmin(pixmap[:, :, 1]))) + pmap_ymax = int(np.ceil(np.nanmax(pixmap[:, :, 1]))) + pixmap = pixmap.copy() + pixmap[:, :, 0] -= pmap_xmin + pixmap[:, :, 1] -= pmap_ymin + self._out_shape = ( + pmap_xmax - pmap_xmin + 1, + pmap_ymax - pmap_ymin + 1 + ) + elif nshapes == 1: + self._out_shape = list(self._output_shapes)[0] + else: + raise ValueError( + "Inconsistent data shapes: 'out_shape' and/or " + "out_img, out_img2, out_wht, out_ctx have different shapes." + ) # pragma: no cover self._alloc_output_arrays( out_shape=self._out_shape, diff --git a/drizzle/tests/test_resample.py b/drizzle/tests/test_resample.py index 935328c..9d4345a 100644 --- a/drizzle/tests/test_resample.py +++ b/drizzle/tests/test_resample.py @@ -1477,8 +1477,6 @@ def test_drizzle_weights_squared(kernel, fc): out_shape = (int(y.max()) + 1, int(x.max()) + 1) - assert np.min(x) > -0.5 and np.min(y) > -0.5 - if fc: # create a Drizzle object using all default parameters # (except for 'kernel', 'out_shape') @@ -1512,15 +1510,14 @@ def test_drizzle_weights_squared(kernel, fc): assert len(driz.out_img2) == 1 else: - # create a Drizzle object using all default parameters - # (except for 'kernel') + # create a Drizzle object using mostly default parameters driz = resample.Drizzle( kernel=kernel, + out_img2=[None], fillval2=-99.0, ) assert driz.out_img is None - assert driz.out_img2 is None assert driz.total_exptime == 0.0 with pytest.warns(Warning): @@ -1553,6 +1550,7 @@ def test_drizzle_weights_squared(kernel, fc): rtol=1.0e-6, atol=0.0 ) + assert abs(float(driz.fillval2) + 99.0) < 1e-7 def test_drizzle_weights_squared_bad_inputs(): @@ -1695,6 +1693,7 @@ def test_drizzle_weights_squared_array_shape_mismatch(): y, x = np.indices(in_shape, dtype=np.float64) in_sci1 = np.zeros(in_shape, dtype=np.float32) + in_sci1[n // 2, n // 2] = 2.222222222222 in_sci1_sq = np.zeros(in_shape, dtype=np.float32) in_wht2 = np.zeros(in_shape1, dtype=np.float32) @@ -1719,7 +1718,8 @@ def test_drizzle_weights_squared_array_shape_mismatch(): driz = resample.Drizzle( kernel=kernel, - out_img2=[out_img2, out_img2], + out_img=out_img2.copy(), + out_img2=[out_img2, out_img2, None], ) with pytest.raises(ValueError) as err_info: driz.add_image( @@ -1727,7 +1727,7 @@ def test_drizzle_weights_squared_array_shape_mismatch(): exptime=1.0, pixmap=pixmap, weight_map=in_wht2, - data2=[in_sci1_sq, in_sci2_sq], + data2=[in_sci1_sq, in_sci2_sq, None], ) assert str(err_info.value).startswith( "'data2' shape(s) is not consistent with 'data' shape." @@ -1772,3 +1772,24 @@ def test_drizzle_weights_squared_array_shape_mismatch(): assert str(err_info.value).startswith( "'weight_map' shape is not consistent with 'data' shape." ) + + # zero-sized variance array + driz = resample.Drizzle( + kernel=kernel, + out_img2=[out_img2, out_img2.copy(), out_img2.copy(), None] + ) + driz.add_image( + data=in_sci1, + exptime=1.0, + pixmap=pixmap, + data2=[in_sci1, in_sci1, np.array([]), None] + ) + driz.add_image( + data=in_sci1, + exptime=1.0, + pixmap=pixmap, + data2=[in_sci1, None, in_sci1, None] + ) + assert np.allclose(np.nansum(driz.out_img2[0]), 2.0 * np.nansum(driz.out_img2[1])) + assert np.allclose(np.nansum(driz.out_img2[0]), 2.0 * np.nansum(driz.out_img2[2])) + assert np.allclose(0.0, np.nansum(driz.out_img2[3])) diff --git a/src/cdrizzlebox.c b/src/cdrizzlebox.c index 5fbe151..8dfa505 100644 --- a/src/cdrizzlebox.c +++ b/src/cdrizzlebox.c @@ -269,7 +269,7 @@ do_kernel_point(struct driz_param_t *p) { struct scanner s; integer_t i, j, ii, jj, k; integer_t osize[2]; - float scale2, d, dow, *d2 = NULL; + float scale2, scale4, d, dow, *d2 = NULL; integer_t bv; int xmin, xmax, ymin, ymax, n; int ndata2; @@ -288,6 +288,8 @@ do_kernel_point(struct driz_param_t *p) { get_dimensions(p->output_data, osize); if (ndata2 > 0) { + scale4 = scale2 * scale2; + if (!(d2 = (float *)malloc(p->ndata2 * sizeof(float)))) { driz_error_set(p->error, PyExc_MemoryError, "Memory allocation failed."); @@ -346,13 +348,11 @@ do_kernel_point(struct driz_param_t *p) { } else { /* Allow for stretching because of scale change */ d = get_pixel(p->data, i, j) * scale2; - if (d2) { - for (k = 0; k < ndata2; ++k) { - if (p->data2[k]) { - d2[k] = get_pixel(p->data2[k], i, j) * scale2; - } else { - d2[k] = 0.0f; - } + for (k = 0; k < ndata2; ++k) { + if (p->data2[k]) { + d2[k] = get_pixel(p->data2[k], i, j) * scale4; + } else { + d2[k] = 0.0f; } } @@ -398,7 +398,7 @@ do_kernel_gaussian(struct driz_param_t *p) { integer_t osize[2]; float d, dow, *d2 = NULL; double gaussian_efac, gaussian_es; - double pfo, ac, scale2, xxi, xxa, yyi, yya, w, ddx, ddy, r2, dover; + double pfo, ac, scale2, scale4, xxi, xxa, yyi, yya, w, ddx, ddy, r2, dover; const double nsig = 2.5; int xmin, xmax, ymin, ymax, n; int ndata2; @@ -429,6 +429,7 @@ do_kernel_gaussian(struct driz_param_t *p) { get_dimensions(p->output_data, osize); if (ndata2 > 0) { + scale4 = scale2 * scale2; if (!(d2 = (float *)malloc(p->ndata2 * sizeof(float)))) { driz_error_set(p->error, PyExc_MemoryError, "Memory allocation failed."); @@ -491,13 +492,11 @@ do_kernel_gaussian(struct driz_param_t *p) { /* Allow for stretching because of scale change */ d = get_pixel(p->data, i, j) * scale2; - if (d2) { - for (k = 0; k < ndata2; ++k) { - if (p->data2[k]) { - d2[k] = get_pixel(p->data2[k], i, j) * scale2; - } else { - d2[k] = 0.0f; - } + for (k = 0; k < ndata2; ++k) { + if (p->data2[k]) { + d2[k] = get_pixel(p->data2[k], i, j) * scale4; + } else { + d2[k] = 0.0f; } } @@ -560,7 +559,7 @@ do_kernel_lanczos(struct driz_param_t *p) { struct scanner s; integer_t bv, i, j, ii, jj, k, nxi, nxa, nyi, nya, nhit, ix, iy; integer_t osize[2]; - float scale2, d, dow, *d2 = NULL; + float scale2, scale4, d, dow, *d2 = NULL; double pfo, xx, yy, xxi, xxa, yyi, yya, w, dx, dy, dover; int kernel_order; struct lanczos_param_t lanczos; @@ -600,6 +599,8 @@ do_kernel_lanczos(struct driz_param_t *p) { get_dimensions(p->output_data, osize); if (ndata2 > 0) { + scale4 = scale2 * scale2; + if (!(d2 = (float *)malloc(p->ndata2 * sizeof(float)))) { driz_error_set(p->error, PyExc_MemoryError, "Memory allocation failed."); @@ -659,13 +660,11 @@ do_kernel_lanczos(struct driz_param_t *p) { /* Allow for stretching because of scale change */ d = get_pixel(p->data, i, j) * scale2; - if (d2) { - for (k = 0; k < ndata2; ++k) { - if (p->data2[k]) { - d2[k] = get_pixel(p->data2[k], i, j) * scale2; - } else { - d2[k] = 0.0f; - } + for (k = 0; k < ndata2; ++k) { + if (p->data2[k]) { + d2[k] = get_pixel(p->data2[k], i, j) * scale4; + } else { + d2[k] = 0.0f; } } @@ -736,7 +735,7 @@ do_kernel_turbo(struct driz_param_t *p) { integer_t bv, i, j, ii, jj, k, nxi, nxa, nyi, nya, nhit, iis, iie, jjs, jje; integer_t osize[2]; float d, dow, *d2 = NULL; - double pfo, scale2, ac; + double pfo, scale2, scale4, ac; double xxi, xxa, yyi, yya, w, dover; int xmin, xmax, ymin, ymax, n; int ndata2; @@ -759,6 +758,8 @@ do_kernel_turbo(struct driz_param_t *p) { get_dimensions(p->output_data, osize); if (ndata2 > 0) { + scale4 = scale2 * scale2; + if (!(d2 = (float *)malloc(ndata2 * sizeof(float)))) { driz_error_set(p->error, PyExc_MemoryError, "Memory allocation failed."); @@ -828,13 +829,11 @@ do_kernel_turbo(struct driz_param_t *p) { /* Allow for stretching because of scale change */ d = get_pixel(p->data, i, j) * scale2; - if (d2) { - for (k = 0; k < ndata2; ++k) { - if (p->data2[k]) { - d2[k] = get_pixel(p->data2[k], i, j) * scale2; - } else { - d2[k] = 0.0f; - } + for (k = 0; k < ndata2; ++k) { + if (p->data2[k]) { + d2[k] = get_pixel(p->data2[k], i, j) * scale4; + } else { + d2[k] = 0.0f; } } @@ -899,7 +898,7 @@ int do_kernel_square(struct driz_param_t *p) { integer_t bv, i, j, ii, jj, k, min_ii, max_ii, min_jj, max_jj, nhit; integer_t osize[2]; - float scale2, d, dow, *d2 = NULL; + float scale2, scale4, d, dow, *d2 = NULL; double dh, jaco, tem, dover, w; double xin[4], yin[4], xout[4], yout[4]; int ndata2; @@ -926,6 +925,8 @@ do_kernel_square(struct driz_param_t *p) { get_dimensions(p->output_data, osize); if (ndata2 > 0) { + scale4 = scale2 * scale2; + if (!(d2 = (float *)malloc(p->ndata2 * sizeof(float)))) { driz_error_set(p->error, PyExc_MemoryError, "Memory allocation failed."); @@ -1004,13 +1005,11 @@ do_kernel_square(struct driz_param_t *p) { /* Allow for stretching because of scale change */ d = get_pixel(p->data, i, j) * scale2; - if (d2) { - for (k = 0; k < ndata2; ++k) { - if (p->data2[k]) { - d2[k] = get_pixel(p->data2[k], i, j) * scale2; - } else { - d2[k] = 0.0f; - } + for (k = 0; k < ndata2; ++k) { + if (p->data2[k]) { + d2[k] = get_pixel(p->data2[k], i, j) * scale4; + } else { + d2[k] = 0.0f; } }