diff --git a/docs/references.bib b/docs/references.bib index b4a01117..d5426750 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -165,6 +165,18 @@ @article{Fornaro2015CAESARApproachBased keywords = {3-D,4-D and multidimensional (Multi-D) SAR imaging,Covariance matrices,Covariance matrix decomposition,differential SAR tomography,differential synthetic aperture radar (SAR) interferometry (DInSAR),Interferometry,Monitoring,principal component analysis (PCA),SAR interferometry (InSAR),SAR tomography,Scattering,Spatial resolution,Synthetic aperture radar,Tomography} } +@article{Goldstein1998RadarInterferogramFiltering, + title = {Radar Interferogram Filtering for Geophysical Applications}, + author = {Goldstein, Richard M. and Werner, Charles L.}, + year = {1998}, + journal = {Geophysical Research Letters}, + volume = {25}, + number = {21}, + pages = {4035--4038}, + issn = {1944-8007}, + doi = {10.1029/1998GL900033} +} + @article{Guarnieri2008ExploitationTargetStatistics, title = {On the {{Exploitation}} of {{Target Statistics}} for {{SAR Interferometry Applications}}}, author = {Guarnieri, A. M. and Tebaldini, S.}, diff --git a/src/dolphin/goldstein.py b/src/dolphin/goldstein.py index 369846a0..521c4ec5 100644 --- a/src/dolphin/goldstein.py +++ b/src/dolphin/goldstein.py @@ -1,14 +1,16 @@ import numpy as np -from numpy.typing import ArrayLike +from numpy.typing import NDArray -def goldstein(phase: ArrayLike, alpha: float, psize: int = 32) -> np.ndarray: +def goldstein( + phase: NDArray[np.complex_] | NDArray[np.float_], alpha: float, psize: int = 32 +) -> np.ndarray: """Apply the Goldstein adaptive filter to the given data. Parameters ---------- phase : np.ndarray - 2D complex array containing the data to be filtered. + 2D array of floating point phase, or complex data, to be filtered. alpha : float Filtering parameter for Goldstein algorithm Must be between 0 (no filtering) and 1 (maximum filtering) @@ -20,18 +22,22 @@ def goldstein(phase: ArrayLike, alpha: float, psize: int = 32) -> np.ndarray: ------- 2D numpy array of filtered data. + References + ---------- + [@Goldstein1998RadarInterferogramFiltering] + """ - def apply_pspec(data: ArrayLike): + def apply_pspec(data: NDArray[np.complex_]) -> np.ndarray: # NaN is allowed value if alpha < 0: raise ValueError(f"alpha must be >= 0, got {alpha = }") - wgt = np.power(np.abs(data) ** 2, alpha / 2) - data = wgt * data + weight = np.power(np.abs(data) ** 2, alpha / 2) + data = weight * data return data - def make_wgt(nxp: int, nyp: int) -> np.ndarray: + def make_weight(nxp: int, nyp: int) -> np.ndarray: # Create arrays of horizontal and vertical weights wx = 1.0 - np.abs(np.arange(nxp // 2) - (nxp / 2.0 - 1.0)) / (nxp / 2.0 - 1.0) wy = 1.0 - np.abs(np.arange(nyp // 2) - (nyp / 2.0 - 1.0)) / (nyp / 2.0 - 1.0) @@ -39,16 +45,16 @@ def make_wgt(nxp: int, nyp: int) -> np.ndarray: # the top-left quadrant of the weight matrix quadrant = np.outer(wy, wx) # Create a full weight matrix by mirroring the quadrant along both axes - wgt = np.block( + weight = np.block( [ [quadrant, np.flip(quadrant, axis=1)], [np.flip(quadrant, axis=0), np.flip(np.flip(quadrant, axis=0), axis=1)], ] ) - return wgt + return weight def patch_goldstein_filter( - data: ArrayLike, wgt: ArrayLike, psize: int + data: NDArray[np.complex_], weight: NDArray[np.float_], psize: int ) -> np.ndarray: """Apply the filter to a single patch of data. @@ -56,7 +62,7 @@ def patch_goldstein_filter( ---------- data : np.ndarray 2D complex array containing the data to be filtered. - wgt : np.ndarray + weight : np.ndarray weight matrix for summing neighboring data psize : int edge length of square FFT area @@ -70,30 +76,39 @@ def patch_goldstein_filter( data = np.fft.fft2(data, s=(psize, psize)) data = apply_pspec(data) data = np.fft.ifft2(data, s=(psize, psize)) - return wgt * data + return weight * data - def apply_goldstein_filter(data: ArrayLike) -> np.ndarray: + def apply_goldstein_filter(data: NDArray[np.complex_]) -> np.ndarray: # Create an empty array for the output out = np.zeros(data.shape, dtype=np.complex64) - # ignore processing for empty chunks - if np.all(np.isnan(data)): + empty_mask = np.isnan(data) | (np.angle(data) == 0) + # ignore processing for a chunks + if np.all(empty_mask): return data # Create the weight matrix - wgt_matrix = make_wgt(psize, psize) + weight_matrix = make_weight(psize, psize) # Iterate over windows of the data for i in range(0, data.shape[0] - psize, psize // 2): for j in range(0, data.shape[1] - psize, psize // 2): - # Create proocessing windows + # Create processing windows data_window = data[i : i + psize, j : j + psize] - wgt_window = wgt_matrix[: data_window.shape[0], : data_window.shape[1]] + weight_window = weight_matrix[ + : data_window.shape[0], : data_window.shape[1] + ] # Apply the filter to the window - filtered_window = patch_goldstein_filter(data_window, wgt_window, psize) + filtered_window = patch_goldstein_filter( + data_window, weight_window, psize + ) # Add the result to the output array slice_i = slice(i, min(i + psize, out.shape[0])) slice_j = slice(j, min(j + psize, out.shape[1])) out[slice_i, slice_j] += filtered_window[ : slice_i.stop - slice_i.start, : slice_j.stop - slice_j.start ] + out[empty_mask] = 0 return out - return apply_goldstein_filter(phase) + if np.iscomplexobj(phase): + return apply_goldstein_filter(phase) + else: + return apply_goldstein_filter(np.exp(1j * phase))