From 1bd5ed6b59af18b842ac82ce3783ed5c94d1a63f Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Wed, 28 Oct 2026 10:51:03 -0400 Subject: [PATCH] Mask phase estimates without covariance support --- src/dolphin/phase_link/_core.py | 9 +++- src/dolphin/phase_link/_ps_filling.py | 15 ++++-- tests/test_phase_link_core.py | 73 +++++++++++++++++++++++++++ 3 files changed, 92 insertions(+), 5 deletions(-) diff --git a/src/dolphin/phase_link/_core.py b/src/dolphin/phase_link/_core.py index 0f3939855..d28f8332d 100644 --- a/src/dolphin/phase_link/_core.py +++ b/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) diff --git a/src/dolphin/phase_link/_ps_filling.py b/src/dolphin/phase_link/_ps_filling.py index 5b86c9cef..b369225ba 100644 --- a/src/dolphin/phase_link/_ps_filling.py +++ b/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: diff --git a/tests/test_phase_link_core.py b/tests/test_phase_link_core.py index 89530808e..6ecd9a586 100644 --- a/tests/test_phase_link_core.py +++ b/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)])