From da95af1971b8ea0069b12182f026a51e92d87be9 Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Tue, 19 Sep 2023 20:40:36 +0000
Subject: [PATCH 01/15] sm r3.4

---
 tools/imagesets/oracle8conda/distrib_nisar/Dockerfile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile b/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile
index c34cb42af..45ecd19c0 100644
--- a/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile
+++ b/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile
@@ -12,7 +12,7 @@ RUN cd /opt \
  && git clone https://$GIT_OAUTH_TOKEN@github-fn.jpl.nasa.gov/NISAR-ADT/SoilMoisture \
  && git clone https://$GIT_OAUTH_TOKEN@github-fn.jpl.nasa.gov/NISAR-ADT/QualityAssurance \
  && cd /opt/QualityAssurance && git checkout v4.0.0 && rm -rf .git \
- && cd /opt/SoilMoisture && git checkout f62fe7b47001aea2195f3c8e88d5f7d3a30e71a7 && rm -rf .git
+ && cd /opt/SoilMoisture && git checkout 93a364c05de2819fce0df704126320cfb5face68 && rm -rf .git
 
 FROM $distrib_img
 

From b6cfbbb18c2c162e4f6c0e08996d3f9531d1295f Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Thu, 23 May 2024 21:23:36 +0000
Subject: [PATCH 02/15] change the SM commit id for R4.0.2

---
 tools/imagesets/oracle8conda/distrib_nisar/Dockerfile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile b/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile
index 27c27c5c2..d118c6252 100644
--- a/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile
+++ b/tools/imagesets/oracle8conda/distrib_nisar/Dockerfile
@@ -12,7 +12,7 @@ RUN cd /opt \
  && git clone https://$GIT_OAUTH_TOKEN@github-fn.jpl.nasa.gov/NISAR-ADT/SoilMoisture \
  && git clone https://$GIT_OAUTH_TOKEN@github-fn.jpl.nasa.gov/NISAR-ADT/QualityAssurance \
  && cd /opt/QualityAssurance && git checkout v9.0.0 && rm -rf .git \
- && cd /opt/SoilMoisture && git checkout d4391e350dd8781a79b2f736ec84bc05967b38d9 && rm -rf .git
+ && cd /opt/SoilMoisture && git checkout 9b437761673d97004b46e22a74ee67cc1b26e280 && rm -rf .git
 
 FROM $distrib_img
 

From e50145ff6b4680607e51380c0f8589e39cc9744a Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 17:19:55 +0000
Subject: [PATCH 03/15] remove redundant imports

---
 python/packages/nisar/workflows/rubbersheet.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index e70d029eb..2cac147ee 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -728,11 +728,6 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
     return offset_az, offset_rg
 
 
-import journal
-import numpy as np
-from scipy import interpolate
-
-
 def _interpolate_offsets(offset, interp_method):
     '''
     Replace NaN in offset with interpolated values

From 43f4d1fce088b0a75f5c0a9b7f82c57e5863db64 Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 17:58:41 +0000
Subject: [PATCH 04/15] add the idw

---
 .../packages/nisar/workflows/rubbersheet.py   | 112 +++++++++++++++++-
 1 file changed, 109 insertions(+), 3 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 2cac147ee..77dd15c91 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -23,7 +23,7 @@
 from nisar.workflows.yaml_argparse import YamlArgparse
 from osgeo import gdal
 from scipy import interpolate, ndimage, signal
-
+from scipy.spatial import cKDTree
 
 def run_rubbersheet_with_polyfit(cfg: dict, output_hdf5: str = None):
     '''
@@ -336,8 +336,12 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                     # If there are residual NaNs, use interpolation to fill residual holes
                     nan_count = np.count_nonzero(np.isnan(offset))
                     if nan_count > 0:
-                        offsets[k] = _interpolate_offsets(offset,
-                                                          rubbersheet_params['interpolation_method'])
+                        if rubbersheet_params['interpolation_method'] == 'idw':
+                            offsets[k] = _interpolate_offsets_by_idw(offset)
+                            offsets[k] =  _interpolate_offsets(offset[k],'linear')
+                        else:
+                            offsets[k] = _interpolate_offsets(offset,
+                                                              rubbersheet_params['interpolation_method'])
                     # If required, filter offsets
                     offsets[k] = _filter_offsets(offsets[k], rubbersheet_params)
                     # Save offsets on disk for resampling
@@ -728,6 +732,108 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
     return offset_az, offset_rg
 
 
+def _interpolate_offsets_by_idw(
+    Z: np.ndarray,
+    power: float = 2.0,
+    number: int | None = 100,
+    radius: float | None = 200.0,
+    eps: float = 1e-12):
+    '''
+    Performs IDW interpolation on a 2D grid containing NaN values.
+    Missing pixels are filled using inverse distance weighting (IDW)
+    based on nearby valid samples. The interpolation can optionally
+    limit the number of nearest neighbors and the maximum search
+    radius to improve efficiency and control smoothing.
+
+    Parameters
+    ---------
+    Z: np.ndarray
+        2D array of shape (H, W) containing NaN values to be filled
+    power: float
+        Power parameter p controlling the distance weighting
+    number: int or None
+        Number of nearest neighbors used for interpolation.
+        If None, all valid points are used (default: 100)
+    radius: float or None
+        Maximum search radius in pixels. If None, no radius
+        limitation is applied (default: 200)
+    eps: float
+        Small value added to distance to avoid division by zero
+
+    Returns
+    -------
+    filled: np.ndarray
+        2D array with the same shape as Z where NaN values
+        have been filled using IDW interpolation
+    '''
+
+    valid = np.isfinite(Z)
+    yy, xx = np.nonzero(valid)
+    if yy.size == 0:
+        raise ValueError("No valid points found to interpolate from.")
+
+    points = np.column_stack([xx, yy]).astype(float)   # (x, y)
+    values = Z[yy, xx].astype(float)
+
+    # Query points (NaNs)
+    qmask = ~valid
+    yq, xq = np.nonzero(qmask)
+    if yq.size == 0:
+        return Z.copy()
+
+    qpoints = np.column_stack([xq, yq]).astype(float)
+
+    # KDTree for neighbors
+    tree = cKDTree(points)
+
+    # If number is None, use all valid points
+    k = points.shape[0] if number is None else number
+
+    # Query neighbors
+    if radius is None:
+        dists, idxs = tree.query(qpoints, k=min(k, points.shape[0]))
+        # dists: (M, k)
+        # idxs : (M, k)
+        use_mask = np.isfinite(dists)
+    else:
+        dists, idxs = tree.query(qpoints, k=min(k, points.shape[0]), distance_upper_bound=radius)
+        # For neighbors not found within radius, idxs == points.shape[0], dists == inf
+        use_mask = np.isfinite(dists) & (idxs < points.shape[0])
+
+    # Filled values
+    filled_vals = np.full(qpoints.shape[0], np.nan, dtype=float)
+
+    # Fill the value with nearest neighbor if distance is too small.
+    zero_hit = np.any((dists <= eps) & use_mask, axis=1)
+    if np.any(zero_hit):
+        row = np.where(zero_hit)[0]
+        col = np.argmax(((dists <= eps) & use_mask)[row], axis=1)
+        filled_vals[row] = values[idxs[row, col]]
+
+    # For the rest, compute IDW
+    rest = ~zero_hit
+    if np.any(rest):
+        d = dists[rest]
+        ii = idxs[rest]
+        m = use_mask[rest]
+
+        # weights = 1 / d^p, ignore invalid neighbors
+        w = np.zeros_like(d, dtype=float)
+        w[m] = 1.0 / np.maximum(d[m], eps) ** power
+
+        v = np.zeros_like(d, dtype=float)
+        v[m] = values[ii[m]]
+
+        ws = w.sum(axis=1)
+        # If no neighbors (ws==0), leave as NaN (no extrapolation)
+        ok = ws > 0
+        filled_vals[np.where(rest)[0][ok]] = (w[ok] * v[ok]).sum(axis=1) / ws[ok]
+
+    out = Z.copy()
+    out[yq, xq] = filled_vals
+
+    return out
+
 def _interpolate_offsets(offset, interp_method):
     '''
     Replace NaN in offset with interpolated values

From ea3acd5ef190e0ed7311a0c37f4232126c5095eb Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 18:14:56 +0000
Subject: [PATCH 05/15] update the schema

---
 python/packages/nisar/workflows/rubbersheet.py |  6 +++++-
 share/nisar/defaults/insar.yaml                | 13 +++++++++++--
 share/nisar/schemas/insar.yaml                 | 11 ++++++++++-
 3 files changed, 26 insertions(+), 4 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 77dd15c91..d32054290 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -337,7 +337,11 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                     nan_count = np.count_nonzero(np.isnan(offset))
                     if nan_count > 0:
                         if rubbersheet_params['interpolation_method'] == 'idw':
-                            offsets[k] = _interpolate_offsets_by_idw(offset)
+                            # Get the IDW parameters
+                            power, number, radius = ( rubbersheet_params['idw_interpolation'].get(k)
+                                                     for k in ('power', 'number', 'radius'))
+                            offsets[k] = _interpolate_offsets_by_idw(offset, power, number, radius)
+                            # Additional linear interpolation is required to fill the extra outliers beyond the radius
                             offsets[k] =  _interpolate_offsets(offset[k],'linear')
                         else:
                             offsets[k] = _interpolate_offsets(offset,
diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml
index 36f3cbc4f..8505fbe7f 100644
--- a/share/nisar/defaults/insar.yaml
+++ b/share/nisar/defaults/insar.yaml
@@ -306,7 +306,7 @@ runconfig:
                     # When `min_cluster_pixels` is 0, this pre-removal step is
                     # skipped entirely.
                     min_cluster_pixels: 2
-                    
+
                     bridge_algorithm_enabled: True
                     # Maximum radius used when bridging disconnected
                     # regions.
@@ -545,6 +545,15 @@ runconfig:
                 # Interpolation method. Interpolation is used to fill residual outlier
                 # holes if present
                 interpolation_method: linear
+                # IDW interpolation method parameters
+                idw_interpolation:
+                    # Power parameter controlling the distance weighting
+                    power: 2
+                    # Number of nearest neighbors used for interpolation
+                    number: 100
+                    # Maximum search radius in pixels
+                    radius: 200
+
                 # It is good practice to filter the offsets prior to interferogram
                 # formation to reduce noise. We expose: no_filter: do not filter the offsets/
                 # degrade offsets resolution. median_filter, boxcar_filter (moving average),
@@ -780,7 +789,7 @@ runconfig:
                     # During deramping, uniform sampling is applied
                     # when number of samples is larger than maximum pixel
                     bridge_ramp_maximum_pixel: 1e6
-    
+
             correction_luts:
                 # Boolean flag to activate/deactivate model-based solid earth tide corrections
                 solid_earth_tides_enabled: True
diff --git a/share/nisar/schemas/insar.yaml b/share/nisar/schemas/insar.yaml
index ad1766f14..c1f089a45 100644
--- a/share/nisar/schemas/insar.yaml
+++ b/share/nisar/schemas/insar.yaml
@@ -548,7 +548,8 @@ rubbersheet_options:
 
     # Interpolation method. Interpolation is used to fill residual
     # outlier holes (if present)
-    interpolation_method: enum('no_interpolation', 'nearest', 'linear', 'cubic', required=False)
+    interpolation_method: enum('no_interpolation', 'nearest', 'linear', 'cubic', 'idw', required=False)
+    idw_interpolation: include('idw_interpolation_options', required=False)
 
     # It is good practice to filter the offsets prior to interferogram
     # formation to reduce noise. We expose: no_filter: do not filter the offsets/
@@ -595,6 +596,14 @@ polyfitting_options:
     # Number of samples  along the azimuth in pixels for polyfitting
     samples_along_azimuth: num(min=0, required=False)
 
+idw_interpolation_options:
+    # Power parameter controlling the distance weighting
+    power: int(min=0, required=False)
+    # Number of nearest neighbors used for interpolation
+    number: int(min=0, required=False)
+    # Maximum search radius in pixels
+    radius: num(min=0, required=False)
+
 crossmul_options:
     # Path to HDF5 file or directory with coregistered SLCs
     # If directory then the following directory tree is required:

From dd786c7b57828077bc7be8753bb28eef02afc614 Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 18:29:58 +0000
Subject: [PATCH 06/15] add the connected components function

---
 .../packages/nisar/workflows/rubbersheet.py   | 87 +++++++++++++++++++
 1 file changed, 87 insertions(+)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index d32054290..2d1eb5660 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -735,6 +735,93 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
 
     return offset_az, offset_rg
 
+def _build_connected_components(
+    arr: np.ndarray,
+    connectivity: int = 4,
+    return_sizes: bool = True):
+    '''
+    Build connected components from a 2D floating-point array.
+
+    Parameters
+    ----------
+    arr : (H, W) float array
+        Input array.
+    connectivity : int
+        Pixel connectivity for labeling (4 or 8).
+    return_sizes : bool
+        If True, return the pixel counts for each component.
+
+    Returns
+    -------
+    labeled : (H, W) int array
+        Labeled component map (0 = background, 1..N = components).
+    num_features : int
+        Number of connected components.
+    valid_mask : (H, W) bool array
+        Boolean mask used to determine foreground pixels.
+    sizes : 1D int array, optional
+        Pixel count for each component label.
+    '''
+    if arr.ndim != 2:
+        raise ValueError(f"Expected a 2D array, got shape {arr.shape}")
+
+    if connectivity not in (4, 8):
+        raise ValueError("connectivity must be 4 or 8")
+
+    valid_mask = ~np.isnan(arr)
+
+    # Choose connectivity structure for 2D
+    if connectivity == 4:
+        structure = np.array([[0, 1, 0],
+                              [1, 1, 1],
+                              [0, 1, 0]], dtype=np.int8)
+    else:  # 8-connectivity
+        structure = np.ones((3, 3), dtype=np.int8)
+
+    labeled, num_features = ndimage.label(valid_mask, structure=structure)
+
+    if return_sizes:
+        sizes = np.bincount(labeled.ravel())  # sizes[0] is background count
+        return labeled, num_features, valid_mask, sizes
+
+    return labeled, num_features, valid_mask
+
+def _remove_small_components(arr: np.ndarray,
+                            labeled: np.ndarray,
+                            sizes: np.ndarray,
+                            min_size: int = 21,
+                            fill_value: float = np.nan):
+    '''
+    Remove small connected components from an array.
+
+    Components with pixel count smaller than min_size are set to fill_value.
+
+    Parameters
+    ----------
+    arr : (H, W) float array
+        Input array to be modified.
+    labeled : (H, W) int array
+        Connected component labels (0 = background), typically from ndimage.label.
+    sizes : 1D int array
+        Pixel counts for each label (e.g., np.bincount(labeled.ravel())).
+    min_size : int
+        Minimum number of pixels required to keep a component.
+    fill_value : float
+        Value used to replace pixels in removed components.
+
+    Returns
+    -------
+    arr : (H, W) float array
+        Array with small components removed.
+    '''
+
+    # labels to remove (exclude background label 0)
+    small_labels = np.where((sizes < min_size) & (np.arange(len(sizes)) != 0))[0]
+    mask_small = np.isin(labeled, small_labels)
+
+    out = arr.copy()
+    out[mask_small] = fill_value
+    return out, small_labels
 
 def _interpolate_offsets_by_idw(
     Z: np.ndarray,

From 2b6b9b63fb5a852cc57980eefa9a67c5bea56ab1 Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 18:30:54 +0000
Subject: [PATCH 07/15] update the default value

---
 share/nisar/defaults/insar.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml
index 8505fbe7f..38d833afa 100644
--- a/share/nisar/defaults/insar.yaml
+++ b/share/nisar/defaults/insar.yaml
@@ -544,7 +544,7 @@ runconfig:
                     kernel_size: 3
                 # Interpolation method. Interpolation is used to fill residual outlier
                 # holes if present
-                interpolation_method: linear
+                interpolation_method: idw
                 # IDW interpolation method parameters
                 idw_interpolation:
                     # Power parameter controlling the distance weighting

From ac487fb6f0a495d3468723c921732e5ace358a2c Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 18:36:25 +0000
Subject: [PATCH 08/15] run isort

---
 python/packages/nisar/workflows/rubbersheet.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 2d1eb5660..976e321ec 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -25,6 +25,7 @@
 from scipy import interpolate, ndimage, signal
 from scipy.spatial import cKDTree
 
+
 def run_rubbersheet_with_polyfit(cfg: dict, output_hdf5: str = None):
     '''
     Run rubbersheet via polyfitting method

From 5318933cd81b86a08123972fac12e8724102f60a Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 20:58:48 +0000
Subject: [PATCH 09/15] fix bugs

---
 python/packages/nisar/workflows/rubbersheet.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 976e321ec..eca689228 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -14,7 +14,7 @@
 from isce3.math import offsets_polyfit
 from nisar.products.insar.product_paths import RIFGGroupsPaths
 from nisar.products.readers import SLC
-from nisar.workflows import prepare_insar_hdf5
+from nisar.workflows import h5_prep
 from nisar.workflows.compute_stats import compute_stats_real_hdf5_dataset
 from nisar.workflows.helpers import (get_cfg_freq_pols,
                                      get_ground_track_velocity_product,
@@ -343,7 +343,7 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                                                      for k in ('power', 'number', 'radius'))
                             offsets[k] = _interpolate_offsets_by_idw(offset, power, number, radius)
                             # Additional linear interpolation is required to fill the extra outliers beyond the radius
-                            offsets[k] =  _interpolate_offsets(offset[k],'linear')
+                            offsets[k] =  _interpolate_offsets(offsets[k],'linear')
                         else:
                             offsets[k] = _interpolate_offsets(offset,
                                                               rubbersheet_params['interpolation_method'])
@@ -1044,5 +1044,5 @@ def run(cfg: dict, output_hdf5: str = None):
 
     # Prepare RIFG. Culled offsets will be
     # allocated in RIFG product
-    out_paths = prepare_insar_hdf5.run(rubbersheet_runconfig.cfg)
+    _, out_paths = h5_prep.get_products_and_paths(rubbersheet_runconfig.cfg)
     run(rubbersheet_runconfig.cfg, out_paths['RIFG'])

From 44c2902b8041c0500de327a9440d6a247f5a7807 Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 21:06:55 +0000
Subject: [PATCH 10/15] update the crossmul main

---
 python/packages/nisar/workflows/crossmul.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/packages/nisar/workflows/crossmul.py b/python/packages/nisar/workflows/crossmul.py
index 4b79fb65a..62d8ce2d4 100644
--- a/python/packages/nisar/workflows/crossmul.py
+++ b/python/packages/nisar/workflows/crossmul.py
@@ -14,7 +14,7 @@
 
 from isce3.io import HDF5OptimizedReader
 from nisar.products.readers import SLC
-from nisar.workflows import prepare_insar_hdf5
+from nisar.workflows import h5_prep
 from nisar.workflows.compute_stats import compute_stats_real_data
 from nisar.workflows.crossmul_runconfig import CrossmulRunConfig
 from nisar.workflows.helpers import (copy_raster,
@@ -220,6 +220,6 @@ def stats_offsets(h5_ds, freq, pol):
     # get a runconfig dict from command line args
     crossmul_runconfig = CrossmulRunConfig(args, resample_type)
     # prepare RIFG HDF5
-    out_paths = prepare_insar_hdf5.run(crossmul_runconfig.cfg)
+    _, out_paths = h5_prep.get_products_and_paths(crossmul_runconfig.cfg)
     # run crossmul
-    run(crossmul_runconfig.cfg, out_paths['RIFG'], resample_type)
+    run(crossmul_runconfig.cfg, out_paths['RIFG'], resample_type, dump_on_disk=True)

From f1631544bddc6d56606961e5a345d010591b0dbb Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Wed, 18 Mar 2026 22:22:30 +0000
Subject: [PATCH 11/15] add subswath mask

---
 .../packages/nisar/workflows/rubbersheet.py   | 30 ++++++++++++++-----
 1 file changed, 23 insertions(+), 7 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index eca689228..6cb427fef 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -12,6 +12,7 @@
 import numpy as np
 from isce3.io import HDF5OptimizedReader
 from isce3.math import offsets_polyfit
+from isce3.unwrap.preprocess import interpret_subswath_mask
 from nisar.products.insar.product_paths import RIFGGroupsPaths
 from nisar.products.readers import SLC
 from nisar.workflows import h5_prep
@@ -292,6 +293,12 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                                                                            zero_doppler_time,
                                                                            dem_file,
                                                                            rubbersheet_dir)
+            # Apply the subswath mask for the interpolation
+            pixel_offsets_mask_path = f'{pixel_offsets_path}/mask'
+            subswath_mask = dst_h5[pixel_offsets_mask_path][()]
+            ref_valid, sec_valid, _ = interpret_subswath_mask(subswath_mask)
+            valid_mask = ref_valid & sec_valid
+
             for pol in pol_list:
                 # Create input and output directories for pol under processing
                 pol_group_path = f'{pixel_offsets_path}/{pol}'
@@ -306,7 +313,7 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                     # Identify outliers
                     offset_az_culled, offset_rg_culled = identify_outliers(
                         str(dense_offsets_dir),
-                        rubbersheet_params)
+                        rubbersheet_params, valid_mask)
                     # Fill outliers holes
                     offset_az = fill_outliers_holes(offset_az_culled,
                                                     rubbersheet_params)
@@ -326,7 +333,9 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                                   key.startswith('layer')]
                     # Apply offset blending
                     offset_az, offset_rg = _offset_blending(off_product_dir,
-                                                            rubbersheet_params, layer_keys)
+                                                            rubbersheet_params,
+                                                            layer_keys,
+                                                            valid_mask)
 
                     # Get correlation peak path for the first offset layer
                     corr_peak_path = str(f'{off_prod_dir}/{layer_keys[0]}/correlation_peak')
@@ -443,7 +452,7 @@ def _write_to_disk(outpath, array, format='ENVI',
     ds.GetRasterBand(1).WriteArray(array)
     ds.FlushCache()
 
-def identify_outliers(offsets_dir, rubbersheet_params):
+def identify_outliers(offsets_dir, rubbersheet_params, mask = None):
     '''
     Identify outliers in the offset fields.
     Outliers are identified by a thresholding
@@ -457,6 +466,8 @@ def identify_outliers(offsets_dir, rubbersheet_params):
         where pixel offsets are located
     rubbersheet_params: cfg
         Dictionary containing rubbersheet parameters
+    mask: np.ndarray
+        A mask indicating the interested region (default: None)
 
     Returns
     -------
@@ -477,6 +488,9 @@ def identify_outliers(offsets_dir, rubbersheet_params):
     elif metric == 'median_filter':
         offset_az = _open_raster(f'{offsets_dir}/dense_offsets', 1)
         offset_rg = _open_raster(f'{offsets_dir}/dense_offsets', 2)
+        if mask is not None:
+            offset_az[~mask] = np.nan
+            offset_rg[~mask] = np.nan
         mask_data = compute_mad_mask(offset_az, window_az, window_rg, threshold) | \
                     compute_mad_mask(offset_rg, window_az, window_rg, threshold)
     elif metric == 'covariance':
@@ -680,7 +694,7 @@ def _fill_nan_with_mean(arr_in, arr_out, neighborhood_size):
     return filled_arr
 
 
-def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
+def _offset_blending(off_product_dir, rubbersheet_params, layer_keys, mask = None):
     '''
     Blends offsets layers at different resolution. Implements a
     pyramidal filling algorithm using the offset layer at higher
@@ -698,6 +712,8 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
         Dictionary containing the user-defined rubbersheet options
     layer_keys: list
         List of layers within the offset product
+    mask: np.ndarray
+        A mask indicting the interested region, default: None
 
     Returns
     -------
@@ -709,7 +725,7 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
 
     # Filter outliers from layer one
     offset_az, offset_rg = identify_outliers(str(off_product_dir / layer_keys[0]),
-                                             rubbersheet_params)
+                                             rubbersheet_params, mask)
 
     # Replace the NaN locations in layer1 with the mean of pixels in layers
     # at lower resolution computed in a neighborhood centered at the NaN location
@@ -720,12 +736,12 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys):
 
         if nan_count_az > 0:
             offset_az_culled, _ = identify_outliers(str(off_product_dir / layer_key),
-                                                    rubbersheet_params)
+                                                    rubbersheet_params, mask)
             offset_az = _fill_nan_with_mean(offset_az, offset_az_culled, filter_size)
 
         if nan_count_rg > 0:
             _, offset_rg_culled = identify_outliers(str(off_product_dir / layer_key),
-                                                    rubbersheet_params)
+                                                    rubbersheet_params, mask)
             offset_rg = _fill_nan_with_mean(offset_rg, offset_rg_culled, filter_size)
 
     # Fill remaining holes by iteratively filling the output offset layer

From f301f716db10cdb2b1168eb0b859a34890f3b97f Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Thu, 19 Mar 2026 18:28:04 +0000
Subject: [PATCH 12/15] add the subswath mask to the runconfig

---
 python/packages/nisar/workflows/rubbersheet.py | 12 +++++++-----
 share/nisar/defaults/insar.yaml                |  4 ++++
 share/nisar/schemas/insar.yaml                 |  5 +++++
 3 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 6cb427fef..59228a8ff 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -293,11 +293,13 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                                                                            zero_doppler_time,
                                                                            dem_file,
                                                                            rubbersheet_dir)
-            # Apply the subswath mask for the interpolation
-            pixel_offsets_mask_path = f'{pixel_offsets_path}/mask'
-            subswath_mask = dst_h5[pixel_offsets_mask_path][()]
-            ref_valid, sec_valid, _ = interpret_subswath_mask(subswath_mask)
-            valid_mask = ref_valid & sec_valid
+            valid_mask = None
+            # Apply the subswath mask to the outlier detection
+            if rubbersheet_params['subswath_mask_apply_enabled']:
+                pixel_offsets_mask_path = f'{pixel_offsets_path}/mask'
+                subswath_mask = dst_h5[pixel_offsets_mask_path][()]
+                ref_valid, sec_valid, _ = interpret_subswath_mask(subswath_mask)
+                valid_mask = ref_valid & sec_valid
 
             for pol in pol_list:
                 # Create input and output directories for pol under processing
diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml
index 38d833afa..0757c5035 100644
--- a/share/nisar/defaults/insar.yaml
+++ b/share/nisar/defaults/insar.yaml
@@ -545,6 +545,10 @@ runconfig:
                 # Interpolation method. Interpolation is used to fill residual outlier
                 # holes if present
                 interpolation_method: idw
+
+                # Flag to enable the use of subswath mask for the outlier detection
+                subswath_mask_apply_enabled: True
+
                 # IDW interpolation method parameters
                 idw_interpolation:
                     # Power parameter controlling the distance weighting
diff --git a/share/nisar/schemas/insar.yaml b/share/nisar/schemas/insar.yaml
index c1f089a45..5e2dcf95c 100644
--- a/share/nisar/schemas/insar.yaml
+++ b/share/nisar/schemas/insar.yaml
@@ -549,8 +549,13 @@ rubbersheet_options:
     # Interpolation method. Interpolation is used to fill residual
     # outlier holes (if present)
     interpolation_method: enum('no_interpolation', 'nearest', 'linear', 'cubic', 'idw', required=False)
+
+    # Inverse distance weighting (IDW) interpolation method
     idw_interpolation: include('idw_interpolation_options', required=False)
 
+    # Flag to enable the use of the subswath mask for the outlier detection
+    subswath_mask_apply_enabled: bool(required=False)
+
     # It is good practice to filter the offsets prior to interferogram
     # formation to reduce noise. We expose: no_filter: do not filter the offsets/
     # degrade offsets resolution. median_filter, boxcar_filter (moving average),

From 972331033cce5324eeade71db25c19b3b453001e Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Tue, 7 Apr 2026 17:53:27 +0000
Subject: [PATCH 13/15] remove unused functions

---
 .../packages/nisar/workflows/rubbersheet.py   | 88 -------------------
 1 file changed, 88 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 59228a8ff..56792eab8 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -754,94 +754,6 @@ def _offset_blending(off_product_dir, rubbersheet_params, layer_keys, mask = Non
 
     return offset_az, offset_rg
 
-def _build_connected_components(
-    arr: np.ndarray,
-    connectivity: int = 4,
-    return_sizes: bool = True):
-    '''
-    Build connected components from a 2D floating-point array.
-
-    Parameters
-    ----------
-    arr : (H, W) float array
-        Input array.
-    connectivity : int
-        Pixel connectivity for labeling (4 or 8).
-    return_sizes : bool
-        If True, return the pixel counts for each component.
-
-    Returns
-    -------
-    labeled : (H, W) int array
-        Labeled component map (0 = background, 1..N = components).
-    num_features : int
-        Number of connected components.
-    valid_mask : (H, W) bool array
-        Boolean mask used to determine foreground pixels.
-    sizes : 1D int array, optional
-        Pixel count for each component label.
-    '''
-    if arr.ndim != 2:
-        raise ValueError(f"Expected a 2D array, got shape {arr.shape}")
-
-    if connectivity not in (4, 8):
-        raise ValueError("connectivity must be 4 or 8")
-
-    valid_mask = ~np.isnan(arr)
-
-    # Choose connectivity structure for 2D
-    if connectivity == 4:
-        structure = np.array([[0, 1, 0],
-                              [1, 1, 1],
-                              [0, 1, 0]], dtype=np.int8)
-    else:  # 8-connectivity
-        structure = np.ones((3, 3), dtype=np.int8)
-
-    labeled, num_features = ndimage.label(valid_mask, structure=structure)
-
-    if return_sizes:
-        sizes = np.bincount(labeled.ravel())  # sizes[0] is background count
-        return labeled, num_features, valid_mask, sizes
-
-    return labeled, num_features, valid_mask
-
-def _remove_small_components(arr: np.ndarray,
-                            labeled: np.ndarray,
-                            sizes: np.ndarray,
-                            min_size: int = 21,
-                            fill_value: float = np.nan):
-    '''
-    Remove small connected components from an array.
-
-    Components with pixel count smaller than min_size are set to fill_value.
-
-    Parameters
-    ----------
-    arr : (H, W) float array
-        Input array to be modified.
-    labeled : (H, W) int array
-        Connected component labels (0 = background), typically from ndimage.label.
-    sizes : 1D int array
-        Pixel counts for each label (e.g., np.bincount(labeled.ravel())).
-    min_size : int
-        Minimum number of pixels required to keep a component.
-    fill_value : float
-        Value used to replace pixels in removed components.
-
-    Returns
-    -------
-    arr : (H, W) float array
-        Array with small components removed.
-    '''
-
-    # labels to remove (exclude background label 0)
-    small_labels = np.where((sizes < min_size) & (np.arange(len(sizes)) != 0))[0]
-    mask_small = np.isin(labeled, small_labels)
-
-    out = arr.copy()
-    out[mask_small] = fill_value
-    return out, small_labels
-
 def _interpolate_offsets_by_idw(
     Z: np.ndarray,
     power: float = 2.0,

From f350df397126ce8cbf4766df494a9ab89cbb4221 Mon Sep 17 00:00:00 2001
From: xhuang-jpl <118782850+xhuang-jpl@users.noreply.github.com>
Date: Thu, 9 Apr 2026 14:59:02 -0700
Subject: [PATCH 14/15] Apply suggestions from code review

Co-authored-by: xhuang-jpl <118782850+xhuang-jpl@users.noreply.github.com>
---
 python/packages/nisar/workflows/rubbersheet.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 56792eab8..64aef7c04 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -350,8 +350,8 @@ def run_rubbersheet_with_interpolation(cfg: dict, output_hdf5: str = None):
                     if nan_count > 0:
                         if rubbersheet_params['interpolation_method'] == 'idw':
                             # Get the IDW parameters
-                            power, number, radius = ( rubbersheet_params['idw_interpolation'].get(k)
-                                                     for k in ('power', 'number', 'radius'))
+                            power, number, radius = ( rubbersheet_params['idw_interpolation'].get(i)
+                                                     for i in ('power', 'number', 'radius'))
                             offsets[k] = _interpolate_offsets_by_idw(offset, power, number, radius)
                             # Additional linear interpolation is required to fill the extra outliers beyond the radius
                             offsets[k] =  _interpolate_offsets(offsets[k],'linear')

From 6a895d5e6409718d5acd7333eea48f87e5852947 Mon Sep 17 00:00:00 2001
From: Xiaodong Huang <xhuang@nisar-adt-dev-3.jpl.nasa.gov>
Date: Thu, 9 Apr 2026 23:12:32 +0000
Subject: [PATCH 15/15] update all comments

---
 .../packages/nisar/workflows/rubbersheet.py   | 71 +++++++------------
 share/nisar/defaults/insar.yaml               |  2 +-
 share/nisar/schemas/insar.yaml                | 11 +--
 3 files changed, 32 insertions(+), 52 deletions(-)

diff --git a/python/packages/nisar/workflows/rubbersheet.py b/python/packages/nisar/workflows/rubbersheet.py
index 64aef7c04..b74ec2bb4 100644
--- a/python/packages/nisar/workflows/rubbersheet.py
+++ b/python/packages/nisar/workflows/rubbersheet.py
@@ -758,12 +758,11 @@ def _interpolate_offsets_by_idw(
     Z: np.ndarray,
     power: float = 2.0,
     number: int | None = 100,
-    radius: float | None = 200.0,
-    eps: float = 1e-12):
+    radius: float | None = 200.0):
     '''
-    Performs IDW interpolation on a 2D grid containing NaN values.
-    Missing pixels are filled using inverse distance weighting (IDW)
-    based on nearby valid samples. The interpolation can optionally
+    Performs Inverse Distance Weighting (IDW) interpolation on a 2D grid containing NaN values.
+    Missing pixels are filled using IDW interpolation based on nearby
+    valid samples. The interpolation can optionally
     limit the number of nearest neighbors and the maximum search
     radius to improve efficiency and control smoothing.
 
@@ -794,10 +793,9 @@ def _interpolate_offsets_by_idw(
     if yy.size == 0:
         raise ValueError("No valid points found to interpolate from.")
 
-    points = np.column_stack([xx, yy]).astype(float)   # (x, y)
+    points = np.column_stack([xx, yy]).astype(float)
     values = Z[yy, xx].astype(float)
 
-    # Query points (NaNs)
     qmask = ~valid
     yq, xq = np.nonzero(qmask)
     if yq.size == 0:
@@ -805,51 +803,30 @@ def _interpolate_offsets_by_idw(
 
     qpoints = np.column_stack([xq, yq]).astype(float)
 
-    # KDTree for neighbors
     tree = cKDTree(points)
+    k = points.shape[0] if number is None else min(number, points.shape[0])
 
-    # If number is None, use all valid points
-    k = points.shape[0] if number is None else number
+    dists, idxs = tree.query(
+        qpoints, k=k, workers=-1,
+        **({'distance_upper_bound': radius} if radius is not None else {}))
 
-    # Query neighbors
-    if radius is None:
-        dists, idxs = tree.query(qpoints, k=min(k, points.shape[0]))
-        # dists: (M, k)
-        # idxs : (M, k)
-        use_mask = np.isfinite(dists)
-    else:
-        dists, idxs = tree.query(qpoints, k=min(k, points.shape[0]), distance_upper_bound=radius)
-        # For neighbors not found within radius, idxs == points.shape[0], dists == inf
-        use_mask = np.isfinite(dists) & (idxs < points.shape[0])
+    # Ensure 2D — tree.query returns 1D when k=1
+    dists = np.reshape(dists, (qpoints.shape[0], -1))
+    idxs  = np.reshape(idxs,  (qpoints.shape[0], -1))
 
-    # Filled values
-    filled_vals = np.full(qpoints.shape[0], np.nan, dtype=float)
+    use_mask = np.isfinite(dists) & (idxs < points.shape[0])
+
+    # Clip out-of-bounds sentinel indices before indexing into values
+    safe_idxs = np.clip(idxs, 0, points.shape[0] - 1)
 
-    # Fill the value with nearest neighbor if distance is too small.
-    zero_hit = np.any((dists <= eps) & use_mask, axis=1)
-    if np.any(zero_hit):
-        row = np.where(zero_hit)[0]
-        col = np.argmax(((dists <= eps) & use_mask)[row], axis=1)
-        filled_vals[row] = values[idxs[row, col]]
-
-    # For the rest, compute IDW
-    rest = ~zero_hit
-    if np.any(rest):
-        d = dists[rest]
-        ii = idxs[rest]
-        m = use_mask[rest]
-
-        # weights = 1 / d^p, ignore invalid neighbors
-        w = np.zeros_like(d, dtype=float)
-        w[m] = 1.0 / np.maximum(d[m], eps) ** power
-
-        v = np.zeros_like(d, dtype=float)
-        v[m] = values[ii[m]]
-
-        ws = w.sum(axis=1)
-        # If no neighbors (ws==0), leave as NaN (no extrapolation)
-        ok = ws > 0
-        filled_vals[np.where(rest)[0][ok]] = (w[ok] * v[ok]).sum(axis=1) / ws[ok]
+    w = np.where(use_mask, 1.0 / dists**power, 0.0)
+    v = np.where(use_mask, values[safe_idxs], 0.0)
+
+    ws = w.sum(axis=1)
+    ok = ws > 0
+
+    filled_vals = np.full(qpoints.shape[0], np.nan, dtype=float)
+    filled_vals[ok] = (w[ok] * v[ok]).sum(axis=1) / ws[ok]
 
     out = Z.copy()
     out[yq, xq] = filled_vals
diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml
index 0757c5035..f8a36bb6f 100644
--- a/share/nisar/defaults/insar.yaml
+++ b/share/nisar/defaults/insar.yaml
@@ -544,7 +544,7 @@ runconfig:
                     kernel_size: 3
                 # Interpolation method. Interpolation is used to fill residual outlier
                 # holes if present
-                interpolation_method: idw
+                interpolation_method: linear
 
                 # Flag to enable the use of subswath mask for the outlier detection
                 subswath_mask_apply_enabled: True
diff --git a/share/nisar/schemas/insar.yaml b/share/nisar/schemas/insar.yaml
index 5e2dcf95c..9c4e42382 100644
--- a/share/nisar/schemas/insar.yaml
+++ b/share/nisar/schemas/insar.yaml
@@ -602,12 +602,15 @@ polyfitting_options:
     samples_along_azimuth: num(min=0, required=False)
 
 idw_interpolation_options:
-    # Power parameter controlling the distance weighting
+    # Power parameter controlling the distance weighting;
+    # If the power is 0, then it is equally weighted
     power: int(min=0, required=False)
     # Number of nearest neighbors used for interpolation
-    number: int(min=0, required=False)
-    # Maximum search radius in pixels
-    radius: num(min=0, required=False)
+    number: int(min=1, required=False)
+    # Maximum search radius in pixels;
+    # If unspecified, all valid pixels are used,
+    # which can significantly increase runtime and is not recommended
+    radius: num(min=1, required=False)
 
 crossmul_options:
     # Path to HDF5 file or directory with coregistered SLCs