From c4730d2ac47a1ab3f1be2ffa6729a73ff431f60a Mon Sep 17 00:00:00 2001 From: Fabian Jankowski Date: Thu, 7 Sep 2023 10:58:22 +0100 Subject: [PATCH 1/2] Implement new combine function to compute the two dimensional maximum map over all co-added arrays. This is useful for PSF calculations, for instance. --- reproject/mosaicking/coadd.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index 3f2bdc9a9..b54fe5bc8 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -73,7 +73,7 @@ def reproject_and_coadd( `~astropy.io.fits.HDUList` instance, specifies the HDU to use. reproject_function : callable The function to use for the reprojection - combine_function : { 'mean', 'sum', 'median', 'first', 'last' } + combine_function : { 'mean', 'sum', 'median', 'first', 'last', 'max' } The type of function to use for combining the values into the final image. For 'first' and 'last', respectively, the reprojected images are simply overlaid on top of each other. With respect to the order of the @@ -97,8 +97,8 @@ def reproject_and_coadd( # Validate inputs - if combine_function not in ("mean", "sum", "median", "first", "last"): - raise ValueError("combine_function should be one of mean/sum/median/first/last") + if combine_function not in ("mean", "sum", "median", "first", "last", "max"): + raise ValueError("combine_function should be one of mean/sum/median/first/last/max") if reproject_function is None: raise ValueError( @@ -253,6 +253,16 @@ def reproject_and_coadd( final_array[array.view_in_original_array] = np.where( array.footprint > 0, array.array, final_array[array.view_in_original_array] ) + elif combine_function == "max": + for array in arrays: + array.array[array.footprint == 0] = 0 + + old_vals = final_array[array.view_in_original_array] + new_vals = array.array * array.footprint + + final_array[array.view_in_original_array] = np.maximum(old_vals, new_vals) + final_footprint[array.view_in_original_array] += array.footprint + elif combine_function == "median": # Here we need to operate in chunks since we could otherwise run # into memory issues From 8e3c6e743427124fe74c2175736b1fb6b290744e Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 7 Sep 2023 11:41:10 +0100 Subject: [PATCH 2/2] Added combine_function=="min", added test for "min" and "max", and fixed implementation --- reproject/mosaicking/coadd.py | 48 +++++++++++++----------- reproject/mosaicking/tests/test_coadd.py | 21 ++++++++--- 2 files changed, 41 insertions(+), 28 deletions(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index b54fe5bc8..8df9809b5 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -73,7 +73,7 @@ def reproject_and_coadd( `~astropy.io.fits.HDUList` instance, specifies the HDU to use. reproject_function : callable The function to use for the reprojection - combine_function : { 'mean', 'sum', 'median', 'first', 'last', 'max' } + combine_function : { 'mean', 'sum', 'median', 'first', 'last', 'min', 'max' } The type of function to use for combining the values into the final image. For 'first' and 'last', respectively, the reprojected images are simply overlaid on top of each other. With respect to the order of the @@ -97,8 +97,8 @@ def reproject_and_coadd( # Validate inputs - if combine_function not in ("mean", "sum", "median", "first", "last", "max"): - raise ValueError("combine_function should be one of mean/sum/median/first/last/max") + if combine_function not in ("mean", "sum", "median", "first", "last", "min", "max"): + raise ValueError("combine_function should be one of mean/sum/median/first/last/min/max") if reproject_function is None: raise ValueError( @@ -223,6 +223,11 @@ def reproject_and_coadd( final_array = np.zeros(shape_out) final_footprint = np.zeros(shape_out) + if combine_function == "min": + final_array[...] = np.inf + elif combine_function == "max": + final_array[...] = -np.inf + if combine_function in ("mean", "sum"): for array in arrays: # By default, values outside of the footprint are set to NaN @@ -236,32 +241,28 @@ def reproject_and_coadd( if combine_function == "mean": with np.errstate(invalid="ignore"): final_array /= final_footprint - elif combine_function == "first": + + elif combine_function in ("first", "last", "min", "max"): for array in arrays: - mask = final_footprint[array.view_in_original_array] == 0 + if combine_function == "first": + mask = final_footprint[array.view_in_original_array] == 0 + elif combine_function == "last": + mask = array.footprint > 0 + elif combine_function == "min": + mask = (array.footprint > 0) & ( + array.array < final_array[array.view_in_original_array] + ) + elif combine_function == "max": + mask = (array.footprint > 0) & ( + array.array > final_array[array.view_in_original_array] + ) + final_footprint[array.view_in_original_array] = np.where( mask, array.footprint, final_footprint[array.view_in_original_array] ) final_array[array.view_in_original_array] = np.where( mask, array.array, final_array[array.view_in_original_array] ) - elif combine_function == "last": - for array in arrays: - final_footprint[array.view_in_original_array] = np.where( - array.footprint, array.footprint, final_footprint[array.view_in_original_array] - ) - final_array[array.view_in_original_array] = np.where( - array.footprint > 0, array.array, final_array[array.view_in_original_array] - ) - elif combine_function == "max": - for array in arrays: - array.array[array.footprint == 0] = 0 - - old_vals = final_array[array.view_in_original_array] - new_vals = array.array * array.footprint - - final_array[array.view_in_original_array] = np.maximum(old_vals, new_vals) - final_footprint[array.view_in_original_array] += array.footprint elif combine_function == "median": # Here we need to operate in chunks since we could otherwise run @@ -269,4 +270,7 @@ def reproject_and_coadd( raise NotImplementedError("combine_function='median' is not yet implemented") + if combine_function in ("min", "max"): + final_array[final_footprint == 0] = 0.0 + return final_array, final_footprint diff --git a/reproject/mosaicking/tests/test_coadd.py b/reproject/mosaicking/tests/test_coadd.py index e069dd900..1d668992a 100644 --- a/reproject/mosaicking/tests/test_coadd.py +++ b/reproject/mosaicking/tests/test_coadd.py @@ -108,14 +108,17 @@ def test_coadd_with_overlap(self, reproject_function): assert_allclose(array, self.array, atol=ATOL) - @pytest.mark.parametrize("combine_function", ["first", "last"]) + @pytest.mark.parametrize("combine_function", ["first", "last", "min", "max"]) def test_coadd_with_overlap_first_last(self, reproject_function, combine_function): views = self._overlapping_views input_data = self._get_tiles(views) # Make each of the overlapping tiles different for i, (array, wcs) in enumerate(input_data): - input_data[i] = (np.full_like(array, i), wcs) + # We give each tile integer values that range from 0 to 19 but we + # deliberately don't make the first one 0 and the last one 19 so + # that min/max differs from first/last. + input_data[i] = (np.full_like(array, (i + 7) % 20), wcs) array, footprint = reproject_and_coadd( input_data, @@ -127,17 +130,23 @@ def test_coadd_with_overlap_first_last(self, reproject_function, combine_functio # Test that either the correct tile sets the output value in the overlap regions test_sequence = list(enumerate(views)) + if combine_function == "last": test_sequence = test_sequence[::-1] + elif combine_function == "min": + test_sequence = test_sequence[13:] + test_sequence[:13] + elif combine_function == "max": + test_sequence = (test_sequence[13:] + test_sequence[:13])[::-1] + for i, view in test_sequence: - # Each tile in test_sequence should overwrite teh following tiles - # in the overlap regions. We'll use nans to mark pixels in the - # output array that have already been set by a preceeding tile, so + # Each tile in test_sequence should overwrite the following tiles + # in the overlap regions. We'll use NaNs to mark pixels in the + # output array that have already been set by a preceding tile, so # we'll go through, check that each tile matches the non-nan pixels # in its region, and then set that whole region to nan. output_tile = array[view] output_values = output_tile[np.isfinite(output_tile)] - assert_equal(output_values, i) + assert_equal(output_values, (i + 7) % 20) array[view] = np.nan def test_coadd_background_matching(self, reproject_function):