Setup gaussian_filter_nan in utils, use in estimate_correlation_from_phase

src/dolphin/filtering.py

@@ -2,6 +2,7 @@
 from concurrent.futures import ProcessPoolExecutor
 from itertools import repeat
 from pathlib import Path
+from typing import Sequence
 
 import numpy as np
 from numpy.typing import ArrayLike, NDArray
@@ -259,3 +260,48 @@ def _filter_and_save(
     create_image_overviews(output_name, resampling=Resampling.AVERAGE)
 
     return output_name
+
+
+def gaussian_filter_nan(
+    image: ArrayLike, sigma: float | Sequence[float], mode="constant", **kwargs
+) -> np.ndarray:
+    """Apply a gaussian filter to an image with NaNs (avoiding all nans).
+
+    The scipy.ndimage `gaussian_filter` will make the output all NaNs if
+    any of the pixels in the input that touches the kernel is NaN
+
+    Source:
+    https://stackoverflow.com/a/36307291
+
+    Parameters
+    ----------
+    image : ndarray
+        Image with nans to filter
+    sigma : float
+        Size of filter kernel. passed into `gaussian_filter`
+    mode : str, default = "constant"
+        Boundary mode for `[scipy.ndimage.gaussian_filter][]`
+    **kwargs : Any
+        Passed into `[scipy.ndimage.gaussian_filter][]`
+
+    Returns
+    -------
+    ndarray
+        Filtered version of `image`.
+
+    """
+    from scipy.ndimage import gaussian_filter
+
+    if np.sum(np.isnan(image)) == 0:
+        return gaussian_filter(image, sigma=sigma, mode=mode, **kwargs)
+
+    V = image.copy()
+    nan_idxs = np.isnan(image)
+    V[nan_idxs] = 0
+    V_filt = gaussian_filter(V, sigma, **kwargs)
+
+    W = np.ones(image.shape)
+    W[nan_idxs] = 0
+    W_filt = gaussian_filter(W, sigma, **kwargs)
+
+    return V_filt / W_filt

src/dolphin/interferogram.py

@@ -14,11 +14,11 @@
 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 tqdm.contrib.concurrent import thread_map
 
 from dolphin import io, utils
 from dolphin._types import DateOrDatetime, Filename, T
+from dolphin.filtering import gaussian_filter_nan
 
 gdal.UseExceptions()
 
@@ -608,7 +608,7 @@ def __eq__(self, other):
 
 
 def estimate_correlation_from_phase(
-    ifg: Union[VRTInterferogram, ArrayLike], window_size: Union[int, tuple[int, int]]
+    ifg: Union[VRTInterferogram, ArrayLike], window_size: int | tuple[int, int]
 ) -> np.ndarray:
     """Estimate correlation from only an interferogram (no SLCs/magnitudes).
 
@@ -622,10 +622,9 @@ def estimate_correlation_from_phase(
         Interferogram to estimate correlation from.
         If a VRTInterferogram, will load and take the phase.
         If `ifg` is complex, will normalize to unit magnitude before estimating.
-    window_size : Union[int, tuple[int, int]]
+    window_size : int | tuple[int, int]
         Size of window to use for correlation estimation.
-        If int, will use a square window of that size.
-        If tuple, the rectangular window has shape  `size=(row_size, col_size)`.
+        If int, will use a symmetrical Gaussian sigma.
 
     Returns
     -------
@@ -643,10 +642,18 @@ def estimate_correlation_from_phase(
     else:
         # If they passed complex, normalize to unit magnitude
         inp = np.exp(1j * np.nan_to_num(np.angle(ifg)))
+    # Now set the 0s to nans (so that either 0 input, or nan input gets ignored)
+    inp[inp == 0] = np.nan
 
-    # Note: the clipping is from possible partial windows producing correlation
+    # The gaussian window becomes negligible at about "window_size/3":
+    sigma: float | tuple[float, float]
+    if isinstance(window_size, (float, int)):
+        sigma = window_size / 3
+    else:
+        sigma = (window_size[0] / 3), (window_size[1] / 3)
+    # Note: clipping is from possible partial windows producing correlation
     # above 1
-    cor = np.clip(np.abs(uniform_filter(inp, window_size, mode="nearest")), 0, 1)
+    cor = np.clip(np.abs(gaussian_filter_nan(inp, sigma, mode="nearest")), 0, 1)
     # Return the input nans to nan
     cor[nan_mask] = np.nan
     # If the input was 0, the correlation is 0
@@ -666,13 +673,14 @@ def estimate_interferometric_correlations(
     """Estimate correlations for a sequence of interferograms.
 
     Will use the same filename base as inputs with a new suffix.
+    Calls `estimate_correlation_from_phase` on each of `ifg_filenames`.
 
     Parameters
     ----------
     ifg_filenames : Sequence[Filename]
         Paths to complex interferogram files.
-    window_size : tuple[int, int]
-        (row, column) size of window to use for estimate.
+    window_size : tuple[int, int] | int
+        sigma, or (row sigma, column sigma) to use for Gaussian filtering.
     out_driver : str, optional
         Name of output GDAL driver, by default "GTiff".
     out_suffix : str, optional

src/dolphin/unwrap/_post_process.py

@@ -1,7 +1,6 @@
 import numpy as np
 from numba import njit, stencil
 from numpy.typing import ArrayLike
-from scipy.ndimage import gaussian_filter
 
 TWOPI = 2 * np.pi
 
@@ -55,49 +54,6 @@ def rewrap_to_twopi(arr: ArrayLike) -> np.ndarray:
     return np.mod(np.pi + arr, TWOPI) - np.pi
 
 
-def gaussian_filter_nan(
-    image: ArrayLike, sigma: float, mode="constant", **kwargs
-) -> np.ndarray:
-    """Apply a gaussian filter to an image with NaNs (avoiding all nans).
-
-    The scipy.ndimage `gaussian_filter` will make the output all NaNs if
-    any of the pixels in the input that touches the kernel is NaN
-
-    Source:
-    https://stackoverflow.com/a/36307291
-
-    Parameters
-    ----------
-    image : ndarray
-        Image with nans to filter
-    sigma : float
-        Size of filter kernel. passed into `gaussian_filter`
-    mode : str, default = "constant"
-        Boundary mode for `[scipy.ndimage.gaussian_filter][]`
-    **kwargs : Any
-        Passed into `[scipy.ndimage.gaussian_filter][]`
-
-    Returns
-    -------
-    ndarray
-        Filtered version of `image`.
-
-    """
-    if np.sum(np.isnan(image)) == 0:
-        return gaussian_filter(image, sigma=sigma, mode=mode, **kwargs)
-
-    V = image.copy()
-    nan_idxs = np.isnan(image)
-    V[nan_idxs] = 0
-    V_filt = gaussian_filter(V, sigma, **kwargs)
-
-    W = np.ones(image.shape)
-    W[nan_idxs] = 0
-    W_filt = gaussian_filter(W, sigma, **kwargs)
-
-    return V_filt / W_filt
-
-
 def _get_ambiguities(unw: ArrayLike, round_decimals: int = 4) -> np.ndarray:
     mod_2pi_image = np.mod(np.pi + unw, TWOPI) - np.pi
     re_wrapped = np.round(mod_2pi_image, round_decimals)
@@ -107,6 +63,8 @@ def _get_ambiguities(unw: ArrayLike, round_decimals: int = 4) -> np.ndarray:
 def _fill_masked_ambiguities(
     amb_image: ArrayLike, mask: ArrayLike, filter_sigma: int = 60
 ) -> np.ndarray:
+    from dolphin.filtering import gaussian_filter_nan
+
     masked_ambs = amb_image.copy()
     masked_ambs[mask] = np.nan
     ambs_filled = np.round(gaussian_filter_nan(amb_image, filter_sigma))