Mask phase estimates without covariance support

src/dolphin/phase_link/_core.py

@@ -41,7 +41,7 @@ class PhaseLinkOutput(NamedTuple):
     """Output of the MLE solver."""
 
     cpx_phase: np.ndarray
-    """Estimated linked phase."""
+    """Estimated linked phase; NaN where no estimate or PS fill is available."""
 
     temp_coh: np.ndarray
     """Temporal coherence of the optimization.
@@ -360,6 +360,13 @@ def run_cpl(
     # Get the temporal coherence
     temp_coh = metrics.estimate_temp_coh(cpx_phase, C_arrays)
 
+    # An unsupported acquisition has a zero diagonal. Without this mask the
+    # solver can return zero phase and unit temporal coherence for an empty window.
+    no_support = jnp.any(jnp.diagonal(C_arrays, axis1=-2, axis2=-1) == 0, axis=-1)
+    cpx_phase = jnp.where(no_support[..., None], jnp.nan, cpx_phase)
+    crlb_std_dev = jnp.where(no_support[..., None], jnp.nan, crlb_std_dev)
+    temp_coh = jnp.where(no_support, jnp.nan, temp_coh)
+
     # Reshape the (rows, cols, nslcs) output to be same as input stack
     cpx_phase_reshaped = jnp.moveaxis(cpx_phase, -1, 0)
     crlb_std_dev_reshaped = jnp.moveaxis(crlb_std_dev, -1, 0)

src/dolphin/phase_link/_ps_filling.py

@@ -31,6 +31,7 @@ def fill_ps_pixels(
     ----------
     cpx_phase : ndarray, shape = (nslc, rows, cols)
         The complex valued-MLE estimate of the phase.
+        PS pixels retain its magnitude, or use unit magnitude where it is NaN.
     temp_coh : ndarray, shape = (rows, cols)
         The temporal coherence of the estimate.
     slc_stack : ArrayLike
@@ -75,7 +76,7 @@ def fill_ps_pixels(
         ref = np.exp(-1j * np.angle(slc_stack[reference_idx][slc_r_idxs, slc_c_idxs]))
         for i in range(len(slc_stack)):
             slc_phase = np.angle(slc_stack[i][slc_r_idxs, slc_c_idxs])
-            cur_amp = np.abs(cpx_phase[i][ps_mask_looked])
+            cur_amp = _amp_or_one(cpx_phase[i][ps_mask_looked])
             new_value = np.exp(1j * slc_phase) * ref
             cpx_phase[i][ps_mask_looked] = cur_amp * new_value
 
@@ -88,14 +89,20 @@ def fill_ps_pixels(
         ps_phases = np.angle(avg_ps[:, ps_mask_looked])
 
         # Set the angle only, don't change magnitude
-        cpx_phase[:, ps_mask_looked] = np.abs(cpx_phase[:, ps_mask_looked]) * np.exp(
-            1j * ps_phases
-        )
+        cpx_phase[:, ps_mask_looked] = _amp_or_one(
+            cpx_phase[:, ps_mask_looked]
+        ) * np.exp(1j * ps_phases)
 
     # Force PS pixels to have high temporal coherence
     temp_coh[ps_mask_looked] = 1
 
 
+def _amp_or_one(cpx: np.ndarray) -> np.ndarray:
+    """Return unit magnitude where the estimate is NaN."""
+    amp = np.abs(cpx)
+    return np.where(np.isnan(amp), 1.0, amp)
+
+
 def _get_avg_ps(
     slc_stack: np.ndarray, ps_mask: np.ndarray, strides: Strides
 ) -> np.ndarray:

tests/test_phase_link_core.py

@@ -4,6 +4,7 @@
 import numpy.testing as npt
 import pytest
 
+from dolphin import shp
 from dolphin._types import HalfWindow, Strides
 from dolphin.phase_link import _core, covariance, simulate
 from dolphin.shp._gaussian import estimate_neighbors as gaussian_estimate_neighbors
@@ -142,6 +143,78 @@ def test_run_phase_linking_with_shift(slc_samples):
     assert np.all(~np.isnan(pl_out.temp_coh[1:-1, :]))
 
 
+@pytest.mark.parametrize("with_ps", [False, True])
+@pytest.mark.parametrize("use_max_ps", [False, True])
+def test_run_phase_linking_no_shp_neighbors_at_edge(slc_samples, with_ps, use_max_ps):
+    """Mark an unsupported burst edge as nodata unless a PS fills it."""
+    slc_stack = slc_samples.copy().reshape(NUM_ACQ, 11, 11)
+    strides = Strides(3, 3)
+    # The first 3x3 looked cell keeps only its corner pixel: not all nodata, so the
+    # looked nodata mask keeps it, but its center pixel (1, 1) is nodata
+    nodata_mask = np.zeros((11, 11), dtype=bool)
+    nodata_mask[:3, :3] = True
+    nodata_mask[0, 0] = False
+    slc_stack[:, nodata_mask] = np.nan
+    ps_mask = np.zeros((11, 11), dtype=bool)
+    ps_mask[0, 0] = with_ps
+
+    amp = np.abs(slc_stack)
+    neighbor_arrays = shp.estimate_neighbors(
+        halfwin_rowcol=(5, 5),
+        alpha=0.001,
+        strides=strides,
+        mean=np.mean(amp, axis=0),
+        var=np.var(amp, axis=0),
+        nslc=NUM_ACQ,
+    )
+    assert not neighbor_arrays[0, 0].any()
+
+    pl_out = _core.run_phase_linking(
+        slc_stack,
+        half_window=HalfWindow(5, 5),
+        strides=strides,
+        ps_mask=ps_mask,
+        use_max_ps=use_max_ps,
+        neighbor_arrays=neighbor_arrays,
+    )
+    if with_ps:
+        ps_phase = slc_stack[:, 0, 0] * slc_stack[0, 0, 0].conj()
+        npt.assert_array_almost_equal(
+            np.angle(ps_phase), np.angle(pl_out.cpx_phase[:, 0, 0])
+        )
+        npt.assert_allclose(np.abs(pl_out.cpx_phase[:, 0, 0]), 1, atol=1e-6)
+        assert pl_out.temp_coh[0, 0] == 1
+    else:
+        assert np.isnan(pl_out.cpx_phase[:, 0, 0]).all()
+        assert np.isnan(pl_out.temp_coh[0, 0])
+        assert np.isnan(pl_out.crlb_std_dev[:, 0, 0]).all()
+    # The fully covered cells are unaffected
+    assert np.isfinite(pl_out.cpx_phase[:, 1:, 1:]).all()
+    assert np.isfinite(pl_out.temp_coh[1:, 1:]).all()
+
+
+def test_run_cpl_missing_acquisition_and_valid_zero_phase():
+    """Keep supported zero phase while rejecting a zero-power acquisition."""
+    slc_stack = np.ones((4, 5, 5), dtype=np.complex64)
+    slc_stack[1, 2, 2] = 0
+    neighbors = np.zeros((5, 5, 3, 3), dtype=bool)
+    neighbors[:, :, 1, 1] = True
+
+    result = _core.run_cpl(
+        slc_stack,
+        half_window=HalfWindow(1, 1),
+        strides=Strides(1, 1),
+        neighbor_arrays=neighbors,
+        compute_crlb=False,
+    )
+
+    assert np.isnan(result.cpx_phase[:, 2, 2]).all()
+    assert np.isnan(result.temp_coh[2, 2])
+    assert np.isnan(result.crlb_std_dev[:, 2, 2]).all()
+    npt.assert_allclose(result.cpx_phase[:, 1, 1], 1, atol=1e-6)
+    npt.assert_allclose(result.temp_coh[1, 1], 1, atol=1e-6)
+
+
 @pytest.mark.parametrize("strides", [1, 2, 3, 4, 5])
 @pytest.mark.parametrize("half_window", [5, 11])
 @pytest.mark.parametrize("shape", [(10, 11), (15, 20), (50, 75)])