Add interpolation option for unwrapping

docs/references.bib

@@ -55,8 +55,21 @@ @article{Chen2012IonosphericArtifactsSimultaneous
   keywords = {Advanced Land Observation Satellite (ALOS) PhasedArray L-band Synthetic Aperture Radar (PALSAR),Azimuth,California,correlation signatures,Delay,dual-frequency GPS carrier phase data,global positioning system,Global Positioning System,global positioning system (GPS),GPS observations,Hawaii,Iceland,InSAR images,interferograms,interferometric synthetic aperture radar,interpretability,Ionosphere,ionospheric artifacts,ionospheric delay,ionospheric electromagnetic wave propagation,ionospheric propagation delays,ionospheric total electron content,ionospheric variability,L-band InSAR data,L-band InSAR observations,L-band SAR interferometry,neutral atmospheric delays,phase artifacts,pixel misregistration,radar imaging,radar interferometry,radiowave propagation,SAR acquisition times,Satellites,Spaceborne radar,synthetic aperture length scales,synthetic aperture radar,terrain,total electron content (TEC)}
 }
 
+@article{Chen2015PersistentScattererInterpolation,
+  title = {A Persistent Scatterer Interpolation for Retrieving Accurate Ground Deformation over {{InSAR-decorrelated}} Agricultural Fields},
+  author = {Chen, Jingyi and Zebker, Howard A. and Knight, Rosemary},
+  year = {2015},
+  journal = {Geophysical Research Letters},
+  volume = {42},
+  number = {21},
+  pages = {9294--9301},
+  issn = {19448007},
+  doi = {10.1002/2015GL065031},
+  keywords = {decorrelation,groundwater,InSAR deformation map,persistent scatterers}
+}
+
 @article{Fornaro2015CAESARApproachBased,
-  title = {{{CAESAR}}: {{An Approach Based}} on {{Covariance Matrix Decomposition}} to {{Improve Multibaseline}}{\textendash}{{Multitemporal Interferometric SAR Processing}}},
+  title = {{{CAESAR}}: {{An Approach Based}} on {{Covariance Matrix Decomposition}} to {{Improve Multibaseline}}--{{Multitemporal Interferometric SAR Processing}}},
   shorttitle = {{{CAESAR}}},
   author = {Fornaro, Gianfranco and Verde, Simona and Reale, Diego and Pauciullo, Antonio},
   year = {2015},

src/dolphin/__init__.py

@@ -6,3 +6,4 @@
 from ._show_versions import *
 from ._types import *
 from .goldstein import goldstein
+from .interpolation import interpolate

src/dolphin/interpolation.py

@@ -0,0 +1,137 @@
+import numba
+import numpy as np
+from numpy.typing import ArrayLike
+
+from dolphin._log import get_log
+
+from .similarity import get_circle_idxs
+
+logger = get_log(__name__)
+
+
+def interpolate(
+    ifg: ArrayLike,
+    weights: ArrayLike,
+    weight_cutoff: float = 0.5,
+    num_neighbors: int = 20,
+    max_radius: int = 51,
+    min_radius: int = 0,
+    alpha: float = 0.75,
+) -> np.ndarray:
+    """Interpolate a complex interferogram based on pixel weights.
+
+    Build upon persistent scatterer interpolation used in
+    [@Chen2015PersistentScattererInterpolation] and
+    [@Wang2022AccuratePersistentScatterer] by allowing floating-point weights
+    instead of 0/1 PS weights.
+
+    Parameters
+    ----------
+    ifg : np.ndarray, 2D complex array
+        wrapped interferogram to interpolate
+    weights : 2D float array
+        Array of weights from 0 to 1 indicating how strongly to weigh
+        the ifg values when interpolating.
+        A special case of this is a PS mask where
+            weights[i,j] = True if radar pixel (i,j) is a PS
+            weights[i,j] = False if radar pixel (i,j) is not a PS
+        Can also pass a coherence image to use as weights.
+    weight_cutoff: float
+        Threshold to use on `weights` so that pixels where
+        `weight[i, j] < weight_cutoff` have phase values replaced by
+        an interpolated value.
+        The default is 0.5: pixels with weight less than 0.5 are replaced with a
+        smoothed version of the surrounding pixels.
+    num_neighbors: int (optional)
+        number of nearest PS pixels used for interpolation
+        num_neighbors = 20 by default
+    max_radius : int (optional)
+        maximum radius (in pixels) for PS searching
+        max_radius = 51 by default
+    min_radius : int (optional)
+        minimum radius (in pixels) for PS searching
+        max_radius = 0 by default
+    alpha : float (optional)
+        hyperparameter controlling the weight of PS in interpolation: smaller
+        alpha means more weight is assigned to PS closer to the center pixel.
+        alpha = 0.75 by default
+
+    Returns
+    -------
+    interpolated_ifg : 2D complex array
+        interpolated interferogram with the same amplitude, but different
+        wrapped phase at non-ps pixels.
+
+    """
+    nrow, ncol = weights.shape
+
+    weights_float = np.clip(weights.astype(np.float32), 0, 1)
+    # Ensure weights are between 0 and 1
+    if np.any(weights_float > 1):
+        logger.warning("weights array has values greater than 1. Clipping to 1.")
+    if np.any(weights_float < 0):
+        logger.warning("weights array has negative values. Clipping to 0.")
+    weights_float = np.clip(weights_float, 0, 1)
+
+    interpolated_ifg = np.zeros((nrow, ncol), dtype=np.complex64)
+
+    indices = np.array(
+        get_circle_idxs(max_radius, min_radius=min_radius, sort_output=False)
+    )
+
+    _interp_loop(
+        ifg,
+        weights_float,
+        weight_cutoff,
+        num_neighbors,
+        alpha,
+        indices,
+        interpolated_ifg,
+    )
+    return interpolated_ifg
+
+
+@numba.njit(parallel=True)
+def _interp_loop(
+    ifg, weights, weight_cutoff, num_neighbors, alpha, indices, interpolated_ifg
+):
+    nrow, ncol = weights.shape
+    nindices = len(indices)
+    for r0 in numba.prange(nrow):
+        for c0 in range(ncol):
+            if weights[r0, c0] >= weight_cutoff:
+                interpolated_ifg[r0, c0] = ifg[r0, c0]
+                continue
+
+            csum = 0.0 + 0j
+            counter = 0
+            r2 = np.zeros(num_neighbors, dtype=np.float64)
+            cphase = np.zeros(num_neighbors, dtype=np.complex128)
+
+            for i in range(nindices):
+                idx = indices[i]
+                r = r0 + idx[0]
+                c = c0 + idx[1]
+
+                if (
+                    (r >= 0)
+                    and (r < nrow)
+                    and (c >= 0)
+                    and (c < ncol)
+                    and weights[r, c] >= weight_cutoff
+                ):
+                    # calculate the square distance to the center pixel
+                    r2[counter] = idx[0] ** 2 + idx[1] ** 2
+
+                    cphase[counter] = np.exp(1j * np.angle(ifg[r, c]))
+                    counter += 1
+                    if counter >= num_neighbors:
+                        break
+
+            # `counter` got up to one more than the number of elements
+            # The last one will be the largest radius
+            r2_norm = (r2[counter - 1] ** alpha) / 2
+            for i in range(counter):
+                csum += np.exp(-r2[i] / r2_norm) * cphase[i]
+
+            interpolated_ifg[r0, c0] = np.abs(ifg[r0, c0]) * np.exp(1j * np.angle(csum))

src/dolphin/unwrap/_unwrap.py

@@ -7,7 +7,7 @@
 import numpy as np
 from tqdm.auto import tqdm
 
-from dolphin import goldstein, io
+from dolphin import goldstein, interpolate, io
 from dolphin._log import get_log, log_runtime
 from dolphin._types import Filename
 from dolphin.utils import DummyProcessPoolExecutor, full_suffix
@@ -53,6 +53,9 @@ def run(
     overwrite: bool = False,
     run_goldstein: bool = False,
     alpha: float = 0.5,
+    run_interpolation: bool = False,
+    max_radius: int = 51,
+    interpolation_cor_threshold: float = 0.5,
 ) -> tuple[list[Path], list[Path]]:
     """Run snaphu on all interferograms in a directory.
 
@@ -111,6 +114,14 @@ def run(
         Whether to run Goldstein filtering on interferogram
     alpha : float, optional, default = 0.5
         Alpha parameter for Goldstein filtering
+    run_interpolation : bool, optional, default = False
+        Whether to run interpolation on interferogram
+    max_radius : int, optional, default = 51
+        maximum radius (in pixel) for scatterer searching for interpolation
+    interpolation_cor_threshold : float, optional, default = 0.5
+        Threshold on the correlation raster to use for interpolation.
+        Pixels with less than this value are replaced by a weighted
+        combination of neighboring pixels.
 
     Returns
     -------
@@ -176,6 +187,9 @@ def run(
                 scratchdir=scratchdir,
                 run_goldstein=run_goldstein,
                 alpha=alpha,
+                run_interpolation=run_interpolation,
+                max_radius=max_radius,
+                interpolation_cor_threshold=interpolation_cor_threshold,
             )
             for ifg_file, out_file, cor_file in zip(in_files, out_files, cor_filenames)
         ]
@@ -220,6 +234,9 @@ def unwrap(
     scratchdir: Optional[Filename] = None,
     run_goldstein: bool = False,
     alpha: float = 0.5,
+    run_interpolation: bool = False,
+    max_radius: int = 51,
+    interpolation_cor_threshold: float = 0.5,
 ) -> tuple[Path, Path]:
     """Unwrap a single interferogram using snaphu, isce3, or tophu.
 
@@ -278,6 +295,14 @@ def unwrap(
         Whether to run Goldstein filtering on interferogram
     alpha : float, optional, default = 0.5
         Alpha parameter for Goldstein filtering
+    run_interpolation : bool, optional, default = False
+        Whether to run interpolation on interferogram
+    max_radius : int, optional, default = 51
+        maximum radius (in pixel) for scatterer searching for interpolation
+    interpolation_cor_threshold : float, optional, default = 0.5
+        Threshold on the correlation raster to use for interpolation.
+        Pixels with less than this value are replaced by a weighted
+        combination of neighboring pixels.
 
     Returns
     -------
@@ -306,6 +331,10 @@ def unwrap(
             output_filename=combined_mask_file,
         )
 
+    unwrapper_ifg_filename = Path(ifg_filename)
+    unwrapper_unw_filename = Path(unw_filename)
+    name_change = "."
+
     if run_goldstein:
         suf = Path(unw_filename).suffix
         if suf == ".tif":
@@ -315,29 +344,74 @@ def unwrap(
             driver = "ENVI"
             opts = list(io.DEFAULT_ENVI_OPTIONS)
 
+        name_change = ".filt" + name_change
         # If we're running Goldstein filtering, the intermediate
         # filtered/unwrapped rasters are temporary rasters in the scratch dir.
-        filt_ifg_filename = (
-            Path(scratchdir or ".") / Path(ifg_filename).with_suffix(".filt" + suf).name
+        filt_ifg_filename = Path(scratchdir or ".") / (
+            Path(ifg_filename).stem.split(".")[0] + (name_change + "int" + suf)
+        )
+        filt_unw_filename = Path(
+            str(unw_filename).split(".")[0] + (name_change + "unw" + suf)
         )
-        scratch_unw_filename = Path(unw_filename).with_suffix(".filt.unw" + suf)
 
         ifg = io.load_gdal(ifg_filename)
         logger.info(f"Goldstein filtering {ifg_filename} -> {filt_ifg_filename}")
-        filt_ifg = goldstein(ifg, alpha=alpha)
+        modified_ifg = goldstein(ifg, alpha=alpha)
         logger.info(f"Writing filtered output to {filt_ifg_filename}")
         io.write_arr(
-            arr=filt_ifg,
+            arr=modified_ifg,
             output_name=filt_ifg_filename,
             like_filename=ifg_filename,
             driver=driver,
             options=opts,
         )
         unwrapper_ifg_filename = filt_ifg_filename
-        unwrapper_unw_filename = scratch_unw_filename
-    else:
-        unwrapper_ifg_filename = Path(ifg_filename)
-        unwrapper_unw_filename = Path(unw_filename)
+        unwrapper_unw_filename = filt_unw_filename
+
+    if run_interpolation:
+        suf = Path(ifg_filename).suffix
+        if suf == ".tif":
+            driver = "GTiff"
+            opts = list(io.DEFAULT_TIFF_OPTIONS)
+        else:
+            driver = "ENVI"
+            opts = list(io.DEFAULT_ENVI_OPTIONS)
+
+        pre_interp_ifg_filename = unwrapper_ifg_filename
+        pre_interp_unw_filename = unwrapper_unw_filename
+        name_change = ".interp" + name_change
+
+        # temporarily storing the intermediate interpolated rasters in the scratch dir.
+        interp_ifg_filename = Path(scratchdir or ".") / (
+            pre_interp_ifg_filename.stem.split(".")[0] + (name_change + "int" + suf)
+        )
+        interp_unw_filename = Path(
+            str(pre_interp_unw_filename).split(".")[0] + (name_change + "unw" + suf)
+        )
+
+        ifg = io.load_gdal(pre_interp_ifg_filename)
+        corr = io.load_gdal(corr_filename)
+        logger.info(
+            f"Masking pixels with correlation below {interpolation_cor_threshold}"
+        )
+        coherent_pixel_mask = corr[:] >= interpolation_cor_threshold
+        logger.info(f"Interpolating {pre_interp_ifg_filename} -> {interp_ifg_filename}")
+        modified_ifg = interpolate(
+            ifg=ifg,
+            weights=coherent_pixel_mask,
+            weight_cutoff=interpolation_cor_threshold,
+            max_radius=max_radius,
+        )
+        logger.info(f"Writing interpolated output to {interp_ifg_filename}")
+        io.write_arr(
+            arr=modified_ifg,
+            output_name=interp_ifg_filename,
+            like_filename=ifg_filename,
+            driver=driver,
+            options=opts,
+        )
+        unwrapper_ifg_filename = interp_ifg_filename
+        unwrapper_unw_filename = interp_unw_filename
 
     if unwrap_method == UnwrapMethod.SNAPHU:
         # Pass everything to snaphu-py
@@ -388,15 +462,16 @@ def unwrap(
         filename=conncomp_path, output_nodata=ccl_nodata, like_filename=ifg_filename
     )
 
-    # Transfer ambiguity numbers from filtered unwrapped interferogram
+    # Transfer ambiguity numbers from filtered/interpolated unwrapped interferogram
     # back to original interferogram
-    if run_goldstein:
+    if run_goldstein or run_interpolation:
         logger.info(
-            f"Transferring ambiguity numbers from filtered ifg {scratch_unw_filename}"
+            "Transferring ambiguity numbers from filtered/interpolated"
+            "ifg {scratch_unw_filename}"
         )
-        unw_arr = io.load_gdal(scratch_unw_filename)
+        unw_arr = io.load_gdal(unwrapper_unw_filename)
 
-        final_arr = np.angle(ifg) + (unw_arr - np.angle(filt_ifg))
+        final_arr = np.angle(ifg) + (unw_arr - np.angle(modified_ifg))
 
         io.write_arr(
             arr=final_arr,

src/dolphin/workflows/config/_common.py

@@ -170,6 +170,10 @@ class UnwrapOptions(BaseModel, extra="forbid"):
             "Whether to run Goldstein filtering step on wrapped interferogram."
         ),
     )
+    run_interpolation: bool = Field(
+        False,
+        description=("Whether to run interpolation step on wrapped interferogram."),
+    )
     _directory: Path = PrivateAttr(Path("unwrapped"))
     unwrap_method: UnwrapMethod = UnwrapMethod.SNAPHU
     n_parallel_jobs: int = Field(
@@ -216,13 +220,25 @@ class UnwrapOptions(BaseModel, extra="forbid"):
             "Set wrapped phase/correlation to 0 where mask is 0 before unwrapping. "
         ),
     )
-
     alpha: float = Field(
         0.5,
         description=(
             "(for Goldstein filtering) Power parameter for Goldstein algorithm."
         ),
     )
+    max_radius: int = Field(
+        51,
+        ge=0.0,
+        description=("(for interpolation) maximum radius to find scatterers."),
+    )
+    interpolation_cor_threshold: float = Field(
+        0.5,
+        description=" Threshold on the correlation raster to use for interpolation. "
+        "Pixels with less than this value are replaced by a weighted "
+        "combination of neighboring pixels.",
+        ge=0.0,
+        le=1.0,
+    )
 
     @field_validator("ntiles", "downsample_factor", mode="before")
     @classmethod

src/dolphin/workflows/unwrapping.py

@@ -91,6 +91,9 @@ def run(
         scratchdir=unwrap_scratchdir,
         run_goldstein=unwrap_options.run_goldstein,
         alpha=unwrap_options.alpha,
+        run_interpolation=unwrap_options.run_interpolation,
+        max_radius=unwrap_options.max_radius,
+        interpolation_cor_threshold=unwrap_options.interpolation_cor_threshold,
     )
 
     if add_overviews: