From 80f83fe30ac256e06a0433db52c491c7885c1cb9 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Wed, 31 Jan 2024 14:33:30 -0500 Subject: [PATCH 1/2] Deprecate `moving_window_mean` in favor of `uniform_filter` See https://github.com/OPERA-Cal-Val/OPERA_Applications/pull/35 `uniform_filter` seems to run faster for all input/filter sizes. --- src/dolphin/interferogram.py | 3 +- src/dolphin/utils.py | 53 ++++-------------------------------- tests/test_utils.py | 40 --------------------------- 3 files changed, 8 insertions(+), 88 deletions(-) diff --git a/src/dolphin/interferogram.py b/src/dolphin/interferogram.py index 38e428624..f297eb7cd 100644 --- a/src/dolphin/interferogram.py +++ b/src/dolphin/interferogram.py @@ -12,6 +12,7 @@ from opera_utils import get_dates from osgeo import gdal from pydantic import BaseModel, Field, ValidationInfo, field_validator, model_validator +from scipy.ndimage import uniform_filter from dolphin import io, utils from dolphin._log import get_log @@ -623,7 +624,7 @@ def estimate_correlation_from_phase( # Note: the clipping is from possible partial windows producing correlation # above 1 - cor = np.clip(np.abs(utils.moving_window_mean(inp, window_size)), 0, 1) + cor = np.clip(np.abs(uniform_filter(inp, window_size)), 0, 1) # Return the input nans to nan cor[nan_mask] = np.nan # If the input was 0, the correlation is 0 diff --git a/src/dolphin/utils.py b/src/dolphin/utils.py index a0b3ecc6b..cec39ca7d 100644 --- a/src/dolphin/utils.py +++ b/src/dolphin/utils.py @@ -431,55 +431,14 @@ def get_mem(process): def moving_window_mean( image: ArrayLike, size: Union[int, tuple[int, int]] ) -> np.ndarray: - """Calculate the mean of a moving window of size `size`. + """DEPRECATED: use `scipy.ndimage.uniform_filter` directly.""" # noqa: D401 + from scipy.ndimage import uniform_filter - Parameters - ---------- - image : ndarray - input image - size : int or tuple of int - Window size. If a single int, the window is square. - If a tuple of (row_size, col_size), the window can be rectangular. - - Returns - ------- - ndarray - image the same size as `image`, where each pixel is the mean - of the corresponding window. - """ - if isinstance(size, int): - size = (size, size) - if len(size) != 2: - msg = "size must be a single int or a tuple of 2 ints" - raise ValueError(msg) - if size[0] % 2 == 0 or size[1] % 2 == 0: - msg = "size must be odd in both dimensions" - raise ValueError(msg) - - row_size, col_size = size - row_pad = row_size // 2 - col_pad = col_size // 2 - - # Pad the image with zeros - image_padded = np.pad( - image, ((row_pad + 1, row_pad), (col_pad + 1, col_pad)), mode="constant" - ) - - # Calculate the cumulative sum of the image - integral_img = np.cumsum(np.cumsum(image_padded, axis=0), axis=1) - if not np.iscomplexobj(integral_img): - integral_img = integral_img.astype(float) - - # Calculate the mean of the moving window - # Uses the algorithm from https://en.wikipedia.org/wiki/Summed-area_table - window_mean = ( - integral_img[row_size:, col_size:] - - integral_img[:-row_size, col_size:] - - integral_img[row_size:, :-col_size] - + integral_img[:-row_size, :-col_size] + msg = ( + "`moving_window_mean` is deprecated. Please use `scipy.ndimage.uniform_filter`." ) - window_mean /= row_size * col_size - return window_mean + warnings.warn(msg, category=DeprecationWarning, stacklevel=2) + return uniform_filter(image, size=size) def set_num_threads(num_threads: int): diff --git a/tests/test_utils.py b/tests/test_utils.py index c0e134639..03649d0c0 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,6 +1,5 @@ import numpy as np import numpy.testing as npt -import pytest from dolphin import utils @@ -76,42 +75,3 @@ def test_upsample_nearest(): assert upsampled3d.shape == (3, 4, 4) for img in upsampled3d: npt.assert_array_equal(img, upsampled) - - -def test_moving_window_mean_basic(): - image = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) - result = utils.moving_window_mean(image, 3) - expected = np.array( - [ - [1.33333333, 2.33333333, 1.77777778], - [3.0, 5.0, 3.66666667], - [2.66666667, 4.33333333, 3.11111111], - ] - ) - assert np.allclose(result, expected) - - -def test_moving_window_mean_single_pixel(): - image = np.array([[5]]) - result = utils.moving_window_mean(image, 1) - expected = np.array([[5]]) - assert np.allclose(result, expected) - - -def test_moving_window_mean_even_size(): - image = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) - with pytest.raises(ValueError): - utils.moving_window_mean(image, (2, 2)) - - -def test_moving_window_mean_invalid_size_type(): - image = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) - with pytest.raises(ValueError): - utils.moving_window_mean(image, (1, 2, 3)) - - -def test_moving_window_mean_empty_image(): - image = np.array([[]]) - result = utils.moving_window_mean(image, 1) - expected = np.array([[]]) - assert np.allclose(result, expected) From e4df5c11a8588fc20c1ca77a180d3cffa0cc045e Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Wed, 31 Jan 2024 14:58:18 -0500 Subject: [PATCH 2/2] call `nearest` for mode, now there are not as obvious edge effects --- src/dolphin/interferogram.py | 2 +- tests/test_interferogram.py | 9 +-------- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/dolphin/interferogram.py b/src/dolphin/interferogram.py index f297eb7cd..84db2b3b5 100644 --- a/src/dolphin/interferogram.py +++ b/src/dolphin/interferogram.py @@ -624,7 +624,7 @@ def estimate_correlation_from_phase( # Note: the clipping is from possible partial windows producing correlation # above 1 - cor = np.clip(np.abs(uniform_filter(inp, window_size)), 0, 1) + cor = np.clip(np.abs(uniform_filter(inp, window_size, mode="nearest")), 0, 1) # Return the input nans to nan cor[nan_mask] = np.nan # If the input was 0, the correlation is 0 diff --git a/tests/test_interferogram.py b/tests/test_interferogram.py index 3d8c7594e..832994513 100644 --- a/tests/test_interferogram.py +++ b/tests/test_interferogram.py @@ -253,14 +253,7 @@ def test_manual_indexes(tmp_path, four_slc_files): @pytest.fixture() def expected_3x3_cor(): - # the edges will be less than 1 because of the windowing - return np.array( - [ - [0.44444444, 0.66666667, 0.44444444], - [0.66666667, 1.0, 0.66666667], - [0.44444444, 0.66666667, 0.44444444], - ] - ) + return np.ones((3, 3)) @pytest.mark.parametrize("window_size", [3, (3, 3)])