From 7115a57e8861766b04f949eb0bd8fa1ceaa02267 Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Tue, 12 Mar 2024 15:29:09 -0700
Subject: [PATCH 01/15] add interpolation option for unwrapping

---
 src/dolphin/__init__.py                 |   1 +
 src/dolphin/interpolation.py            | 149 ++++++++++++++++++++++++
 src/dolphin/unwrap/_unwrap.py           |  73 ++++++++++--
 src/dolphin/workflows/config/_common.py |   9 +-
 src/dolphin/workflows/unwrapping.py     |   2 +
 tests/test_unwrap.py                    |  15 +++
 6 files changed, 238 insertions(+), 11 deletions(-)
 create mode 100644 src/dolphin/interpolation.py

diff --git a/src/dolphin/__init__.py b/src/dolphin/__init__.py
index 730d52974..3ba98cbb7 100644
--- a/src/dolphin/__init__.py
+++ b/src/dolphin/__init__.py
@@ -6,3 +6,4 @@
 from ._show_versions import *
 from ._types import *
 from .goldstein import goldstein
+from .interpolation import interpolate
diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
new file mode 100644
index 000000000..267ed2eed
--- /dev/null
+++ b/src/dolphin/interpolation.py
@@ -0,0 +1,149 @@
+import numba
+import numpy as np
+
+from dolphin._log import get_log
+
+from .similarity import get_circle_idxs
+
+logger = get_log(__name__)
+
+
+def interpolate(
+    ifg: np.ndarray,
+    weights: np.ndarray,
+    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.
+
+    Parameters
+    ----------
+    ifg : np.ndarray, 2D complex array
+        wrapped interferogram to interpolate
+    weights : 2D boolean 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.
+    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
+    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}"
+        )
+
+    # Ensure weights are between 0 and 1
+    if np.any(weights.astype(np.float32) > 1):
+        logger.warning("weights array has values greater than 1. Clipping to 1.")
+    if np.any(weights.astype(np.float32) < 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)
+
+    interpolated_ifg = np.zeros((nrow, ncol), dtype=np.complex64)
+
+    indices = np.array(get_circle_idxs(max_radius, min_radius=min_radius))
+
+    _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):
+        # convert linear idx to row, col
+        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
+
+            # TODO : why use the "counter - 1" here to normalize?
+            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))
diff --git a/src/dolphin/unwrap/_unwrap.py b/src/dolphin/unwrap/_unwrap.py
index 036178bd1..b96791db5 100644
--- a/src/dolphin/unwrap/_unwrap.py
+++ b/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
@@ -49,6 +49,8 @@ def run(
     overwrite: bool = False,
     run_goldstein: bool = False,
     alpha: float = 0.5,
+    run_interpolation: bool = False,
+    max_radius: int = 51,
 ) -> tuple[list[Path], list[Path]]:
     """Run snaphu on all interferograms in a directory.
 
@@ -103,6 +105,10 @@ 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
 
     Returns
     -------
@@ -167,6 +173,8 @@ def run(
                 scratchdir=scratchdir,
                 run_goldstein=run_goldstein,
                 alpha=alpha,
+                run_interpolation=run_interpolation,
+                max_radius=max_radius,
             )
             for ifg_file, out_file, cor_file in zip(in_files, out_files, cor_filenames)
         ]
@@ -200,6 +208,8 @@ def unwrap(
     scratchdir: Optional[Filename] = None,
     run_goldstein: bool = False,
     alpha: float = 0.5,
+    run_interpolation: bool = False,
+    max_radius: int = 51,
 ) -> tuple[Path, Path]:
     """Unwrap a single interferogram using snaphu, isce3, or tophu.
 
@@ -258,6 +268,10 @@ 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
 
     Returns
     -------
@@ -266,6 +280,9 @@ def unwrap(
     conncomp_path : Path
         Path to output connected component label file.
 
+    Attention: If both interpolation and goldstein filter are true,
+    only goldstein filter will apply and interpolation will be skipped
+
     """
     if isinstance(downsample_factor, int):
         downsample_factor = (downsample_factor, downsample_factor)
@@ -286,6 +303,9 @@ def unwrap(
             output_filename=combined_mask_file,
         )
 
+    unwrapper_ifg_filename = Path(ifg_filename)
+    unwrapper_unw_filename = Path(unw_filename)
+
     if run_goldstein:
         suf = Path(unw_filename).suffix
         if suf == ".tif":
@@ -304,10 +324,10 @@ def unwrap(
 
         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,
@@ -315,9 +335,41 @@ def unwrap(
         )
         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)
+
+    elif 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)
+
+        # temporarily storing the intermediate interpolated rasters in the scratch dir.
+        interp_ifg_filename = (
+            Path(scratchdir or ".")
+            / Path(ifg_filename).with_suffix(".interp" + suf).name
+        )
+        scratch_unw_filename = Path(unw_filename).with_suffix(".interp.unw" + suf)
+
+        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"Interpolating {ifg_filename} -> {interp_ifg_filename}")
+        modified_ifg = interpolate(
+            ifg=ifg, weights=coherent_pixel_mask, 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 = scratch_unw_filename
 
     if unwrap_method == UnwrapMethod.SNAPHU:
         from ._snaphu_py import unwrap_snaphu_py
@@ -370,15 +422,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)
 
-        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,
diff --git a/src/dolphin/workflows/config/_common.py b/src/dolphin/workflows/config/_common.py
index 119fbb664..6b951abec 100644
--- a/src/dolphin/workflows/config/_common.py
+++ b/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(
@@ -210,13 +214,16 @@ class UnwrapOptions(BaseModel, extra="forbid"):
         "smooth",
         description="Statistical cost mode method for SNAPHU.",
     )
-
     alpha: float = Field(
         0.5,
         description=(
             "(for Goldstein filtering) Power parameter for Goldstein algorithm."
         ),
     )
+    max_radius: int = Field(
+        51,
+        description=("(for interpolation) maximum radius to find scatterers."),
+    )
 
     @field_validator("ntiles", "downsample_factor", mode="before")
     @classmethod
diff --git a/src/dolphin/workflows/unwrapping.py b/src/dolphin/workflows/unwrapping.py
index 9b767bf63..0f34ab9c5 100644
--- a/src/dolphin/workflows/unwrapping.py
+++ b/src/dolphin/workflows/unwrapping.py
@@ -89,6 +89,8 @@ 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,
     )
 
     if add_overviews:
diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index a7140932b..1db7cc8c8 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -130,6 +130,21 @@ def test_goldstein(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
         assert unw_path.exists()
         assert conncomp_path.exists()
 
+    @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
+    def test_interpolation(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
+        # test other init_method
+        unw_filename = tmp_path / "unwrapped.unw.tif"
+        unw_path, conncomp_path = dolphin.unwrap.unwrap(
+            ifg_filename=list_of_gtiff_ifgs[0],
+            corr_filename=corr_raster,
+            unw_filename=unw_filename,
+            nlooks=1,
+            unwrap_method=method,
+            run_interpolation=True,
+        )
+        assert unw_path.exists()
+        assert conncomp_path.exists()
+
 
 class TestUnwrapRun:
     def test_run_gtiff(self, list_of_gtiff_ifgs, corr_raster):

From 1012e7c1248e37b5625c8b4cfd6977bd73805b7c Mon Sep 17 00:00:00 2001
From: Scott Staniewicz <scott.j.staniewicz@jpl.nasa.gov>
Date: Tue, 12 Mar 2024 21:07:20 -0400
Subject: [PATCH 02/15] add reference, convert float weights once

---
 docs/references.bib          | 15 ++++++++++++++-
 src/dolphin/interpolation.py | 26 +++++++++++++-------------
 2 files changed, 27 insertions(+), 14 deletions(-)

diff --git a/docs/references.bib b/docs/references.bib
index c7557bcd5..67f70d1c8 100644
--- a/docs/references.bib
+++ b/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},
diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index 267ed2eed..ad4cc9833 100644
--- a/src/dolphin/interpolation.py
+++ b/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,
     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
     ----------
@@ -58,12 +64,6 @@ def interpolate(
         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
 
@@ -79,13 +79,13 @@ def interpolate(
             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)
 

From 97246065bcfd1e56e6a4c5911b47360e17ea7918 Mon Sep 17 00:00:00 2001
From: Scott Staniewicz <scott.j.staniewicz@jpl.nasa.gov>
Date: Tue, 12 Mar 2024 21:10:49 -0400
Subject: [PATCH 03/15] remove the extra all check

---
 src/dolphin/interpolation.py | 12 ------------
 1 file changed, 12 deletions(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index ad4cc9833..269ff3814 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -67,18 +67,6 @@ def interpolate(
     """
     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_float > 1):

From 1102da17a187f8b015d4ef61e97da9d14ea7f067 Mon Sep 17 00:00:00 2001
From: Scott Staniewicz <scott.j.staniewicz@jpl.nasa.gov>
Date: Tue, 12 Mar 2024 21:11:11 -0400
Subject: [PATCH 04/15] Remove pymp reference

---
 src/dolphin/interpolation.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index 269ff3814..3086a9869 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -98,7 +98,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]

From 1ef74c5df742d9ac4b4efc3cc5025cea07d9cab5 Mon Sep 17 00:00:00 2001
From: Scott Staniewicz <scott.j.staniewicz@jpl.nasa.gov>
Date: Tue, 12 Mar 2024 21:13:00 -0400
Subject: [PATCH 05/15] remove todo

---
 src/dolphin/interpolation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index 3086a9869..9ea073507 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -128,7 +128,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]

From ab26f93794b84ce1fb416d5a6c2adaf10b9ecdae Mon Sep 17 00:00:00 2001
From: Scott Staniewicz <scott.j.staniewicz@jpl.nasa.gov>
Date: Tue, 12 Mar 2024 21:54:10 -0400
Subject: [PATCH 06/15] Pass through the correlation threshold for
 interpolation

The last implementation had a default of 0.0
Since every weight was >=0, no pixels were interpolated.
---
 src/dolphin/interpolation.py            | 16 +++++++---------
 src/dolphin/unwrap/_unwrap.py           | 22 +++++++++++++++++++---
 src/dolphin/workflows/config/_common.py |  9 +++++++++
 src/dolphin/workflows/unwrapping.py     |  1 +
 4 files changed, 36 insertions(+), 12 deletions(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index 9ea073507..4162b6fe6 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -12,11 +12,11 @@
 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,
-    weight_cutoff: float = 0.0,
 ) -> np.ndarray:
     """Interpolate a complex interferogram based on pixel weights.
 
@@ -36,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
@@ -49,14 +55,6 @@ 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
     -------
diff --git a/src/dolphin/unwrap/_unwrap.py b/src/dolphin/unwrap/_unwrap.py
index b96791db5..02074a4cf 100644
--- a/src/dolphin/unwrap/_unwrap.py
+++ b/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(
diff --git a/src/dolphin/workflows/config/_common.py b/src/dolphin/workflows/config/_common.py
index 6b951abec..aded917dc 100644
--- a/src/dolphin/workflows/config/_common.py
+++ b/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
diff --git a/src/dolphin/workflows/unwrapping.py b/src/dolphin/workflows/unwrapping.py
index 0f34ab9c5..babd08dee 100644
--- a/src/dolphin/workflows/unwrapping.py
+++ b/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:

From 46aea756b02b7f4c6419f88bef716b878b53a515 Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Thu, 14 Mar 2024 16:53:27 -0700
Subject: [PATCH 07/15] support to have both interpolation and filtering

---
 src/dolphin/interpolation.py  |  4 +++-
 src/dolphin/unwrap/_unwrap.py | 38 ++++++++++++++++++++---------------
 2 files changed, 25 insertions(+), 17 deletions(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index 4162b6fe6..e7cf7eadb 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -75,7 +75,9 @@ def interpolate(
 
     interpolated_ifg = np.zeros((nrow, ncol), dtype=np.complex64)
 
-    indices = np.array(get_circle_idxs(max_radius, min_radius=min_radius))
+    indices = np.array(
+        get_circle_idxs(max_radius, min_radius=min_radius, sort_output=False)
+    )
 
     _interp_loop(
         ifg,
diff --git a/src/dolphin/unwrap/_unwrap.py b/src/dolphin/unwrap/_unwrap.py
index 02074a4cf..9b3d59b90 100644
--- a/src/dolphin/unwrap/_unwrap.py
+++ b/src/dolphin/unwrap/_unwrap.py
@@ -291,9 +291,6 @@ def unwrap(
     conncomp_path : Path
         Path to output connected component label file.
 
-    Attention: If both interpolation and goldstein filter are true,
-    only goldstein filter will apply and interpolation will be skipped
-
     """
     if isinstance(downsample_factor, int):
         downsample_factor = (downsample_factor, downsample_factor)
@@ -316,6 +313,7 @@ def unwrap(
 
     unwrapper_ifg_filename = Path(ifg_filename)
     unwrapper_unw_filename = Path(unw_filename)
+    name_change = "."
 
     if run_goldstein:
         suf = Path(unw_filename).suffix
@@ -326,12 +324,15 @@ 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}")
@@ -345,9 +346,9 @@ def unwrap(
             options=opts,
         )
         unwrapper_ifg_filename = filt_ifg_filename
-        unwrapper_unw_filename = scratch_unw_filename
+        unwrapper_unw_filename = filt_unw_filename
 
-    elif run_interpolation:
+    if run_interpolation:
         suf = Path(ifg_filename).suffix
         if suf == ".tif":
             driver = "GTiff"
@@ -356,20 +357,25 @@ def unwrap(
             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 ".")
-            / Path(ifg_filename).with_suffix(".interp" + suf).name
+        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)
         )
-        scratch_unw_filename = Path(unw_filename).with_suffix(".interp.unw" + suf)
 
-        ifg = io.load_gdal(ifg_filename)
+        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 {ifg_filename} -> {interp_ifg_filename}")
+        logger.info(f"Interpolating {pre_interp_ifg_filename} -> {interp_ifg_filename}")
         modified_ifg = interpolate(
             ifg=ifg,
             weights=coherent_pixel_mask,
@@ -385,7 +391,7 @@ def unwrap(
             options=opts,
         )
         unwrapper_ifg_filename = interp_ifg_filename
-        unwrapper_unw_filename = scratch_unw_filename
+        unwrapper_unw_filename = interp_unw_filename
 
     if unwrap_method == UnwrapMethod.SNAPHU:
         from ._snaphu_py import unwrap_snaphu_py
@@ -445,7 +451,7 @@ def unwrap(
             "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(modified_ifg))
 

From 1f09151d1b11e1887bfd7fcf105aa4f50fb836b2 Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Mon, 18 Mar 2024 13:38:59 -0700
Subject: [PATCH 08/15] change the values of all zero sample interferogram

---
 tests/conftest.py | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/tests/conftest.py b/tests/conftest.py
index 203a26afe..ac8e4fe9d 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -457,11 +457,21 @@ def dem_file(tmp_path, slc_stack):
 @pytest.fixture()
 def list_of_gtiff_ifgs(tmp_path, raster_100_by_200):
     ifg_list = []
+    constants = [10, 15, 5]
     for i in range(3):
         # Create a copy of the raster in the same directory
         f = tmp_path / f"ifg_{i}.int.tif"
+        array = np.exp(
+            1j
+            * np.arange(
+                -constants[i] * np.pi,
+                constants[i] * np.pi,
+                2 * constants[i] * np.pi / 20000,
+            )
+        ).reshape(100, 200)
+        array[20:70, 50:150] *= constants[i] + 1j * constants[i]
         write_arr(
-            arr=np.ones((100, 200), dtype=np.complex64),
+            arr=array,
             output_name=f,
             like_filename=raster_100_by_200,
             driver="GTiff",

From 04454ea2ef909c54da4911566a27a11ab1aadeec Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Mon, 18 Mar 2024 13:39:45 -0700
Subject: [PATCH 09/15] add test for interpolation loop

---
 tests/test_unwrap.py | 36 +++++++++++++++++++++++++++++++++++-
 1 file changed, 35 insertions(+), 1 deletion(-)

diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index 1db7cc8c8..569f07b06 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -1,3 +1,4 @@
+import math
 import os
 from pathlib import Path
 
@@ -27,8 +28,10 @@ def corr_raster(raster_100_by_200):
     # Make a correlation raster of all 1s in the same directory as the raster
     d = Path(raster_100_by_200).parent
     corr_raster = d / "corr_raster.cor.tif"
+    array = np.ones((100, 200), dtype=np.float32)
+    array[:, 100:102] = 0
     io.write_arr(
-        arr=np.ones((100, 200), dtype=np.float32),
+        arr=array,
         output_name=corr_raster,
         like_filename=raster_100_by_200,
         driver="GTiff",
@@ -119,6 +122,8 @@ def test_unwrap_snaphu_nodata(self, tmp_path, list_of_gtiff_ifgs, corr_raster):
     def test_goldstein(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
         # test other init_method
         unw_filename = tmp_path / "unwrapped.unw.tif"
+        scratch_dir = tmp_path / "scratch"
+        scratch_dir.mkdir()
         unw_path, conncomp_path = dolphin.unwrap.unwrap(
             ifg_filename=list_of_gtiff_ifgs[0],
             corr_filename=corr_raster,
@@ -126,14 +131,42 @@ def test_goldstein(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
             nlooks=1,
             unwrap_method=method,
             run_goldstein=True,
+            scratchdir=scratch_dir,
         )
         assert unw_path.exists()
         assert conncomp_path.exists()
 
+    def test_interp_loop(self, list_of_gtiff_ifgs, corr_raster):
+        ifg = io.load_gdal(list_of_gtiff_ifgs[0])
+        corr = io.load_gdal(corr_raster)
+        interpolated_ifg = np.zeros((100, 200), dtype=np.complex64)
+        indices = np.array(
+            dolphin.similarity.get_circle_idxs(
+                max_radius=51, min_radius=0, sort_output=False
+            )
+        )
+        dolphin.interpolation._interp_loop(
+            ifg,
+            corr,
+            weight_cutoff=0.5,
+            num_neighbors=20,
+            alpha=0.75,
+            indices=indices,
+            interpolated_ifg=interpolated_ifg,
+        )
+        assert math.isclose(
+            np.angle(interpolated_ifg[0, 100]), 0.73608005, rel_tol=0.0000001
+        )
+        assert math.isclose(
+            np.angle(interpolated_ifg[70, 101]), 0.56214905, rel_tol=0.0000001
+        )
+
     @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
     def test_interpolation(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
         # test other init_method
         unw_filename = tmp_path / "unwrapped.unw.tif"
+        scratch_dir = tmp_path / "scratch"
+        scratch_dir.mkdir()
         unw_path, conncomp_path = dolphin.unwrap.unwrap(
             ifg_filename=list_of_gtiff_ifgs[0],
             corr_filename=corr_raster,
@@ -141,6 +174,7 @@ def test_interpolation(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
             nlooks=1,
             unwrap_method=method,
             run_interpolation=True,
+            scratchdir=scratch_dir,
         )
         assert unw_path.exists()
         assert conncomp_path.exists()

From 3eac61a8251908adc570306770d59cd14d420b5e Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Mon, 18 Mar 2024 14:27:13 -0700
Subject: [PATCH 10/15] adjust the equality tolerance

---
 tests/test_unwrap.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index 569f07b06..7859e5e8f 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -155,10 +155,10 @@ def test_interp_loop(self, list_of_gtiff_ifgs, corr_raster):
             interpolated_ifg=interpolated_ifg,
         )
         assert math.isclose(
-            np.angle(interpolated_ifg[0, 100]), 0.73608005, rel_tol=0.0000001
+            np.angle(interpolated_ifg[0, 100]), 0.736080, rel_tol=1e-06
         )
         assert math.isclose(
-            np.angle(interpolated_ifg[70, 101]), 0.56214905, rel_tol=0.0000001
+            np.angle(interpolated_ifg[70, 101]), 0.562149, rel_tol=1e-06
         )
 
     @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])

From 41acf599ea3fd8030b0a0da368ccc5a084430d85 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 18 Mar 2024 21:27:48 +0000
Subject: [PATCH 11/15] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 tests/test_unwrap.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index 7859e5e8f..091b33e26 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -154,9 +154,7 @@ def test_interp_loop(self, list_of_gtiff_ifgs, corr_raster):
             indices=indices,
             interpolated_ifg=interpolated_ifg,
         )
-        assert math.isclose(
-            np.angle(interpolated_ifg[0, 100]), 0.736080, rel_tol=1e-06
-        )
+        assert math.isclose(np.angle(interpolated_ifg[0, 100]), 0.736080, rel_tol=1e-06)
         assert math.isclose(
             np.angle(interpolated_ifg[70, 101]), 0.562149, rel_tol=1e-06
         )

From 4efc8a086e2fba24a5f874ccaee0b284f60228ba Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Mon, 18 Mar 2024 14:36:59 -0700
Subject: [PATCH 12/15] change docstring to match datatype

---
 src/dolphin/interpolation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index e7cf7eadb..85f5ec36d 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -29,7 +29,7 @@ def interpolate(
     ----------
     ifg : np.ndarray, 2D complex array
         wrapped interferogram to interpolate
-    weights : 2D boolean array
+    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

From b242ca0a0d3a6fd78acdfee17b83f65b844e6fa6 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 18 Mar 2024 21:38:17 +0000
Subject: [PATCH 13/15] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/dolphin/interpolation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/dolphin/interpolation.py b/src/dolphin/interpolation.py
index 85f5ec36d..43854f7b3 100644
--- a/src/dolphin/interpolation.py
+++ b/src/dolphin/interpolation.py
@@ -29,7 +29,7 @@ def interpolate(
     ----------
     ifg : np.ndarray, 2D complex array
         wrapped interferogram to interpolate
-    weights : 2D float array 
+    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

From ee7795a511d546ae11f2f2e4652c64cbe660304c Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Mon, 18 Mar 2024 16:19:42 -0700
Subject: [PATCH 14/15] simpler interferogram for unit test

---
 tests/conftest.py    | 9 ++++-----
 tests/test_unwrap.py | 8 +++++---
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/tests/conftest.py b/tests/conftest.py
index ac8e4fe9d..c20f092a9 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -457,19 +457,18 @@ def dem_file(tmp_path, slc_stack):
 @pytest.fixture()
 def list_of_gtiff_ifgs(tmp_path, raster_100_by_200):
     ifg_list = []
-    constants = [10, 15, 5]
+    constants = [1, 2, 3]
     for i in range(3):
         # Create a copy of the raster in the same directory
         f = tmp_path / f"ifg_{i}.int.tif"
         array = np.exp(
             1j
             * np.arange(
-                -constants[i] * np.pi,
-                constants[i] * np.pi,
-                2 * constants[i] * np.pi / 20000,
+                -constants[i],
+                constants[i],
+                2 * constants[i] / 20000,
             )
         ).reshape(100, 200)
-        array[20:70, 50:150] *= constants[i] + 1j * constants[i]
         write_arr(
             arr=array,
             output_name=f,
diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index 091b33e26..97bee67f0 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -29,7 +29,7 @@ def corr_raster(raster_100_by_200):
     d = Path(raster_100_by_200).parent
     corr_raster = d / "corr_raster.cor.tif"
     array = np.ones((100, 200), dtype=np.float32)
-    array[:, 100:102] = 0
+    array[70, 100] = 0
     io.write_arr(
         arr=array,
         output_name=corr_raster,
@@ -154,9 +154,11 @@ def test_interp_loop(self, list_of_gtiff_ifgs, corr_raster):
             indices=indices,
             interpolated_ifg=interpolated_ifg,
         )
-        assert math.isclose(np.angle(interpolated_ifg[0, 100]), 0.736080, rel_tol=1e-06)
+        neighbors_average_value = (
+            np.sum(np.angle(ifg[69:72, 99:102])) - np.angle(ifg[70, 100])
+        ) / 8
         assert math.isclose(
-            np.angle(interpolated_ifg[70, 101]), 0.562149, rel_tol=1e-06
+            np.angle(interpolated_ifg[70, 100]), neighbors_average_value, rel_tol=1e-02
         )
 
     @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])

From e92240e7fb4aa771edd8474abd6ab94013ad1d45 Mon Sep 17 00:00:00 2001
From: mirzaees <smirzaee@aurora.jpl.nasa.gov>
Date: Tue, 19 Mar 2024 09:34:54 -0700
Subject: [PATCH 15/15] change to more transparent interpolation check

---
 tests/conftest.py    | 13 +++++--------
 tests/test_unwrap.py | 30 ++++++++++++++++++++----------
 2 files changed, 25 insertions(+), 18 deletions(-)

diff --git a/tests/conftest.py b/tests/conftest.py
index c20f092a9..f816ccfda 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -456,19 +456,16 @@ def dem_file(tmp_path, slc_stack):
 # For unwrapping/overviews
 @pytest.fixture()
 def list_of_gtiff_ifgs(tmp_path, raster_100_by_200):
+    x, y = np.meshgrid(np.arange(200), np.arange(100))
+    # simulate a simple phase ramp
+    phase = 0.003 * x + 0.002 * y
     ifg_list = []
+    # simulate three interferograms
     constants = [1, 2, 3]
     for i in range(3):
         # Create a copy of the raster in the same directory
         f = tmp_path / f"ifg_{i}.int.tif"
-        array = np.exp(
-            1j
-            * np.arange(
-                -constants[i],
-                constants[i],
-                2 * constants[i] / 20000,
-            )
-        ).reshape(100, 200)
+        array = np.exp(1j * (phase + constants[i]))
         write_arr(
             arr=array,
             output_name=f,
diff --git a/tests/test_unwrap.py b/tests/test_unwrap.py
index 97bee67f0..e7d756d72 100644
--- a/tests/test_unwrap.py
+++ b/tests/test_unwrap.py
@@ -1,4 +1,3 @@
-import math
 import os
 from pathlib import Path
 
@@ -136,15 +135,26 @@ def test_goldstein(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):
         assert unw_path.exists()
         assert conncomp_path.exists()
 
-    def test_interp_loop(self, list_of_gtiff_ifgs, corr_raster):
-        ifg = io.load_gdal(list_of_gtiff_ifgs[0])
-        corr = io.load_gdal(corr_raster)
-        interpolated_ifg = np.zeros((100, 200), dtype=np.complex64)
+    def test_interp_loop(self):
+        x, y = np.meshgrid(np.arange(200), np.arange(100))
+        # simulate a simple phase ramp
+        phase = 0.003 * x + 0.002 * y
+        # interferogram with the simulated phase ramp and with a constant amplitude
+        ifg = np.exp(1j * phase)
+        corr = np.ones(ifg.shape)
+        # mask out the ifg/corr at a given pixel for example at pixel 50,40
+        # index of the pixel of interest
+        x_idx = 50
+        y_idx = 40
+        corr[y_idx, x_idx] = 0
+        # generate indices for the pixels to be used in interpolation
         indices = np.array(
             dolphin.similarity.get_circle_idxs(
                 max_radius=51, min_radius=0, sort_output=False
             )
         )
+        # interpolate pixels with zero in the corr and write to interpolated_ifg
+        interpolated_ifg = np.zeros((100, 200), dtype=np.complex64)
         dolphin.interpolation._interp_loop(
             ifg,
             corr,
@@ -154,12 +164,12 @@ def test_interp_loop(self, list_of_gtiff_ifgs, corr_raster):
             indices=indices,
             interpolated_ifg=interpolated_ifg,
         )
-        neighbors_average_value = (
-            np.sum(np.angle(ifg[69:72, 99:102])) - np.angle(ifg[70, 100])
-        ) / 8
-        assert math.isclose(
-            np.angle(interpolated_ifg[70, 100]), neighbors_average_value, rel_tol=1e-02
+        # expected phase based on the model above used for simulation
+        expected_phase = 0.003 * x_idx + 0.002 * y_idx
+        phase_error = np.angle(
+            interpolated_ifg[y_idx, x_idx] * np.exp(-1j * expected_phase)
         )
+        assert np.allclose(phase_error, 0.0, atol=1e-3)
 
     @pytest.mark.parametrize("method", [UnwrapMethod.SNAPHU, UnwrapMethod.PHASS])
     def test_interpolation(self, tmp_path, list_of_gtiff_ifgs, corr_raster, method):