modifications while getting the weighting working

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/interpolation.py

@@ -1,5 +1,6 @@
 import numba
 import numpy as np
+from numpy.typing import ArrayLike
 
 from dolphin._log import get_log
 
@@ -9,15 +10,20 @@
 
 
 def interpolate(
-    ifg: np.ndarray,
-    weights: np.ndarray,
+    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,
-    weight_cutoff: float = 0.0,
 ) -> np.ndarray:
-    r"""Interpolation of persistent scatterers.
+    """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
     ----------
@@ -30,6 +36,12 @@ def interpolate(
             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
@@ -43,49 +55,23 @@ def interpolate(
         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
-    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.
-        If `weight_cutoff = 0` (default),  All pixels are replaced with a
-        smoothed version of the surrounding pixels.
-        If `weight_cutoff = 1`, only pixels with exactly weight=1
-        are kept, and the rest are replaced with an interpolated value.
 
     Returns
     -------
     interpolated_ifg : 2D complex array
         interpolated interferogram with the same amplitude, but different
         wrapped phase at non-ps pixels.
 
-    References
-    ----------
-    "A persistent scatterer interpolation for retrieving accurate ground
-    deformation over InSAR\-decorrelated agricultural fields"
-    Chen et al., 2015, https://doi.org/10.1002/2015GL065031
-
     """
     nrow, ncol = weights.shape
 
-    if np.all(
-        np.logical_or(
-            np.logical_or(weights == 0, weights == 1),
-            np.logical_or(weights is True, weights is False),
-        )
-    ):
-        logger.info("Binary weights, using PS-like interpolation.")
-    else:
-        logger.info(
-            f"Range of values as weights, using weight_cutoff = {weight_cutoff}"
-        )
-
+    weights_float = np.clip(weights.astype(np.float32), 0, 1)
     # Ensure weights are between 0 and 1
-    if np.any(weights.astype(np.float32) > 1):
+    if np.any(weights_float > 1):
         logger.warning("weights array has values greater than 1. Clipping to 1.")
-    if np.any(weights.astype(np.float32) < 0):
+    if np.any(weights_float < 0):
         logger.warning("weights array has negative values. Clipping to 0.")
-    # Make shared versions of the input arrays to avoid copying in each thread
-    weights_float = np.clip(weights.astype(np.float32), 0, 1)
+    weights_float = np.clip(weights_float, 0, 1)
 
     interpolated_ifg = np.zeros((nrow, ncol), dtype=np.complex64)
 
@@ -110,7 +96,6 @@ def _interp_loop(
     nrow, ncol = weights.shape
     nindices = len(indices)
     for r0 in numba.prange(nrow):
-        # convert linear idx to row, col
         for c0 in range(ncol):
             if weights[r0, c0] >= weight_cutoff:
                 interpolated_ifg[r0, c0] = ifg[r0, c0]
@@ -141,7 +126,8 @@ def _interp_loop(
                     if counter >= num_neighbors:
                         break
 
-            # TODO : why use the "counter - 1" here to normalize?
+            # `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]

src/dolphin/unwrap/_unwrap.py

@@ -51,6 +51,7 @@ def run(
     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.
 
@@ -109,6 +110,10 @@ def run(
         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
     -------
@@ -175,6 +180,7 @@ def run(
                 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)
         ]
@@ -210,6 +216,7 @@ def unwrap(
     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.
 
@@ -272,6 +279,10 @@ def unwrap(
         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
     -------
@@ -354,11 +365,16 @@ def unwrap(
 
         ifg = io.load_gdal(ifg_filename)
         corr = io.load_gdal(corr_filename)
-        logger.info("Masking pixels with correlation above 0.7")
-        coherent_pixel_mask = corr[:] >= 0.7
+        logger.info(
+            f"Masking pixels with correlation below {interpolation_cor_threshold}"
+        )
+        coherent_pixel_mask = corr[:] >= interpolation_cor_threshold
         logger.info(f"Interpolating {ifg_filename} -> {interp_ifg_filename}")
         modified_ifg = interpolate(
-            ifg=ifg, weights=coherent_pixel_mask, max_radius=max_radius
+            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(

src/dolphin/workflows/config/_common.py

@@ -222,8 +222,17 @@ class UnwrapOptions(BaseModel, extra="forbid"):
     )
     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,7 @@ def run(
         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: