From b1aa869259d61e65f51f4100f4a490a2ada24e41 Mon Sep 17 00:00:00 2001 From: Seyeon Jeon Date: Tue, 4 Jun 2024 17:17:26 -0700 Subject: [PATCH 01/23] Add filtering module for unwrapped interferogram --- src/dolphin/filtering.py | 116 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 src/dolphin/filtering.py diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py new file mode 100644 index 00000000..677198ea --- /dev/null +++ b/src/dolphin/filtering.py @@ -0,0 +1,116 @@ +import numpy as np +from numpy.typing import ArrayLike +from scipy import ndimage + + +def filtering( + unw_ifg: ArrayLike, + corr: ArrayLike, + mask_cutoff: float = 0.5, + wavelength_cutoff: float = 50*1e3, + dx: float = 30, + +) -> np.ndarray: + """Filter out signals with spatial wavelength longer than a threshold. + + Parameters + ---------- + unw_ifg : np.ndarray, 2D complex array + unwrapped interferogram to interpolate + corr : 2D float array + Array of interferometric correlation from 0 to 1. + mask_cutoff: float + Threshold to use on `corr` so that pixels where + `corr[i, j] > mask_cutoff` are True and the rest are False. + The default is 0.5. + wavelength_cutoff: float + Spatial wavelength threshold to filter unw_ifg. + Signals with wavelength longer than 'wavelength_curoff' in unw_ifg will be filtered out. + The default is 50*1e3 (m). + dx : float + Pixel spatial spacing. Assume same spacing for x, y axes. + The default is 30 (m). + + + Returns + ------- + filtered_ifg : 2D complex array + filtered interferogram that does not contain signals with spatial wavelength longer than a threshold. + + """ + nrow, ncol = corr.shape + + # Create Boolean mask for corr > mask_cutoff to be True and the rest to be False + mask = (corr>mask_cutoff).astype('bool') + + # Create Boolean mask for Zero-filled boundary area to be False and the rest to be True + mask_boundary = ~(corr==0).astype('bool') + + # Ramp plane fitting + Y = unw_ifg[mask] # get data of non NaN & masked pixels + Xdata = np.argwhere(mask) # get indices of non NaN & masked pixels + X = np.c_[np.ones((len(Xdata))), Xdata] + theta = np.dot(np.dot( np.linalg.pinv(np.dot(X.transpose(), X)), X.transpose()), Y) + X1_, X2_ = np.mgrid[:nrow, :ncol] + X_ = np.hstack(( np.reshape(X1_, (nrow*ncol, 1)) , np.reshape(X2_, (nrow*ncol, 1)) ) ) + X_ = np.hstack(( np.ones((nrow*ncol, 1)) , X_ )) + plane = np.reshape(np.dot(X_, theta), (nrow, ncol)) + + # Replace masked out pixels with the ramp plane + unw_ifg_interp = np.copy(unw_ifg) + unw_ifg_interp[~mask*mask_boundary] = plane[~mask*mask_boundary] + + # Copy the edge pixels for the boundary area before filling them by reflection + EV_fill = np.copy(unw_ifg_interp) + + NW = Xdata[np.argmin(Xdata[:,0])] # Get indices of upper left corner pixel + SE = Xdata[np.argmax(Xdata[:,0])] # Get indices of lower right corner pixel + SW = Xdata[np.argmin(Xdata[:,1])] # Get indices of lower left corner pixel + NE = Xdata[np.argmax(Xdata[:,1])] # Get indices of upper left corner pixel + + for k in range(NW[1], NE[1]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[0:NE[0]+1, k] == 0) # count zeros in North direction + EV_fill[0:n_zeros, k] = EV_fill[n_zeros,k] + for k in range(SW[1], SE[1]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[SW[0]+1:, k] == 0) # count zeros in South direction + EV_fill[-n_zeros:, k] = EV_fill[-n_zeros-1, k] + for k in range(NW[0], SW[0]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[k, 0:NW[1]+1] == 0) # count zeros in North direction + EV_fill[k, 0:n_zeros] = EV_fill[k, n_zeros] + for k in range(NE[0], SE[0]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[k, SE[1]+1:] == 0) # count zeros in North direction + EV_fill[k, -n_zeros:] = EV_fill[k, -n_zeros-1] + + # Fill the boundary area reflecting the pixel values + Reflect_fill = np.copy(EV_fill) + + for k in range(NW[1], NE[1]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[0:NE[0]+1, k] == 0) # count zeros in North direction + Reflect_fill[0:n_zeros, k] = np.flipud(EV_fill[n_zeros:n_zeros+n_zeros,k]) + for k in range(SW[1], SE[1]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[SW[0]+1:, k] == 0) # count zeros in South direction + Reflect_fill[-n_zeros:, k] = np.flipud(EV_fill[-n_zeros-n_zeros:-n_zeros, k]) + for k in range(NW[0], SW[0]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[k, 0:NW[1]+1] == 0) # count zeros in North direction + Reflect_fill[k, 0:n_zeros] = np.flipud(EV_fill[k, n_zeros:n_zeros+n_zeros]) + for k in range(NE[0], SE[0]+1): + n_zeros = np.count_nonzero(unw_ifg_interp[k, SE[1]+1:] == 0) # count zeros in North direction + Reflect_fill[k, -n_zeros:] = np.flipud(EV_fill[k, -n_zeros-n_zeros:-n_zeros]) + + Reflect_fill[0:NW[0], 0:NW[1]] = np.flipud(Reflect_fill[NW[0]:NW[0]+NW[0], 0:NW[1]]) # upper left corner area + Reflect_fill[0:NE[0], NE[1]+1:] = np.fliplr(Reflect_fill[0:NE[0], NE[1]+1-(ncol-NE[1]-1):NE[1]+1]) # upper right corner area + Reflect_fill[SW[0]+1:, 0:SW[1]] = np.fliplr(Reflect_fill[SW[0]+1:, SW[1]:SW[1]+SW[1]]) # lower left corner area + Reflect_fill[SE[0]+1:, SE[1]+1:] = np.flipud(Reflect_fill[SE[0]+1-(nrow-SE[0]-1):SE[0]+1, SE[1]+1:]) # lower right corner area + + # 2D filtering with Gaussian kernel + #wavelength_cutoff: float = 50*1e3, + #dx: float = 30, + cutoff_value = 0.5 # 0 < cutoff_value < 1 + sigma_f = 1/wavelength_cutoff/np.sqrt(np.log(1/cutoff_value)) # fc = sqrt(ln(1/cutoff_value))*sigma_f + sigma_x = 1/np.pi/2/sigma_f + sigma = sigma_x/dx + + lowpass_filtered = ndimage.gaussian_filter(Reflect_fill,sigma) + filtered_ifg = unw_ifg - lowpass_filtered*mask_boundary + + return filtered_ifg \ No newline at end of file From 99a431fed65f0d112f9078847a9408f888a1e1c6 Mon Sep 17 00:00:00 2001 From: Seyeon Jeon Date: Wed, 5 Jun 2024 14:23:37 -0700 Subject: [PATCH 02/23] spelling error --- src/dolphin/filtering.py | 150 ++++++++++++++++++++++++--------------- 1 file changed, 92 insertions(+), 58 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 677198ea..27f84466 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -7,9 +7,8 @@ def filtering( unw_ifg: ArrayLike, corr: ArrayLike, mask_cutoff: float = 0.5, - wavelength_cutoff: float = 50*1e3, + wavelength_cutoff: float = 50 * 1e3, dx: float = 30, - ) -> np.ndarray: """Filter out signals with spatial wavelength longer than a threshold. @@ -24,8 +23,9 @@ def filtering( `corr[i, j] > mask_cutoff` are True and the rest are False. The default is 0.5. wavelength_cutoff: float - Spatial wavelength threshold to filter unw_ifg. - Signals with wavelength longer than 'wavelength_curoff' in unw_ifg will be filtered out. + Spatial wavelength threshold to filter unw_ifg. + Signals with wavelength longer than 'wavelength_cutoff' in unw_ifg will be + filtered out. The default is 50*1e3 (m). dx : float Pixel spatial spacing. Assume same spacing for x, y axes. @@ -35,82 +35,116 @@ def filtering( Returns ------- filtered_ifg : 2D complex array - filtered interferogram that does not contain signals with spatial wavelength longer than a threshold. + filtered interferogram that does not contain signals with spatial wavelength + longer than a threshold. """ nrow, ncol = corr.shape # Create Boolean mask for corr > mask_cutoff to be True and the rest to be False - mask = (corr>mask_cutoff).astype('bool') + mask = (corr > mask_cutoff).astype("bool") - # Create Boolean mask for Zero-filled boundary area to be False and the rest to be True - mask_boundary = ~(corr==0).astype('bool') + # Create Boolean mask for Zero-filled boundary area to be False + # and the rest to be True + mask_boundary = ~(corr == 0).astype("bool") # Ramp plane fitting - Y = unw_ifg[mask] # get data of non NaN & masked pixels - Xdata = np.argwhere(mask) # get indices of non NaN & masked pixels + Y = unw_ifg[mask] # get data of non NaN & masked pixels + Xdata = np.argwhere(mask) # get indices of non NaN & masked pixels X = np.c_[np.ones((len(Xdata))), Xdata] - theta = np.dot(np.dot( np.linalg.pinv(np.dot(X.transpose(), X)), X.transpose()), Y) + theta = np.dot(np.dot(np.linalg.pinv(np.dot(X.transpose(), X)), X.transpose()), Y) X1_, X2_ = np.mgrid[:nrow, :ncol] - X_ = np.hstack(( np.reshape(X1_, (nrow*ncol, 1)) , np.reshape(X2_, (nrow*ncol, 1)) ) ) - X_ = np.hstack(( np.ones((nrow*ncol, 1)) , X_ )) + X_ = np.hstack( + (np.reshape(X1_, (nrow * ncol, 1)), np.reshape(X2_, (nrow * ncol, 1))) + ) + X_ = np.hstack((np.ones((nrow * ncol, 1)), X_)) plane = np.reshape(np.dot(X_, theta), (nrow, ncol)) # Replace masked out pixels with the ramp plane unw_ifg_interp = np.copy(unw_ifg) - unw_ifg_interp[~mask*mask_boundary] = plane[~mask*mask_boundary] + unw_ifg_interp[~mask * mask_boundary] = plane[~mask * mask_boundary] # Copy the edge pixels for the boundary area before filling them by reflection EV_fill = np.copy(unw_ifg_interp) - NW = Xdata[np.argmin(Xdata[:,0])] # Get indices of upper left corner pixel - SE = Xdata[np.argmax(Xdata[:,0])] # Get indices of lower right corner pixel - SW = Xdata[np.argmin(Xdata[:,1])] # Get indices of lower left corner pixel - NE = Xdata[np.argmax(Xdata[:,1])] # Get indices of upper left corner pixel - - for k in range(NW[1], NE[1]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[0:NE[0]+1, k] == 0) # count zeros in North direction - EV_fill[0:n_zeros, k] = EV_fill[n_zeros,k] - for k in range(SW[1], SE[1]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[SW[0]+1:, k] == 0) # count zeros in South direction - EV_fill[-n_zeros:, k] = EV_fill[-n_zeros-1, k] - for k in range(NW[0], SW[0]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[k, 0:NW[1]+1] == 0) # count zeros in North direction + NW = Xdata[np.argmin(Xdata[:, 0])] # Get indices of upper left corner pixel + SE = Xdata[np.argmax(Xdata[:, 0])] # Get indices of lower right corner pixel + SW = Xdata[np.argmin(Xdata[:, 1])] # Get indices of lower left corner pixel + NE = Xdata[np.argmax(Xdata[:, 1])] # Get indices of upper left corner pixel + + for k in range(NW[1], NE[1] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[0 : NE[0] + 1, k] == 0 + ) # count zeros in North direction + EV_fill[0:n_zeros, k] = EV_fill[n_zeros, k] + for k in range(SW[1], SE[1] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[SW[0] + 1 :, k] == 0 + ) # count zeros in South direction + EV_fill[-n_zeros:, k] = EV_fill[-n_zeros - 1, k] + for k in range(NW[0], SW[0] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[k, 0 : NW[1] + 1] == 0 + ) # count zeros in North direction EV_fill[k, 0:n_zeros] = EV_fill[k, n_zeros] - for k in range(NE[0], SE[0]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[k, SE[1]+1:] == 0) # count zeros in North direction - EV_fill[k, -n_zeros:] = EV_fill[k, -n_zeros-1] + for k in range(NE[0], SE[0] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[k, SE[1] + 1 :] == 0 + ) # count zeros in North direction + EV_fill[k, -n_zeros:] = EV_fill[k, -n_zeros - 1] # Fill the boundary area reflecting the pixel values Reflect_fill = np.copy(EV_fill) - for k in range(NW[1], NE[1]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[0:NE[0]+1, k] == 0) # count zeros in North direction - Reflect_fill[0:n_zeros, k] = np.flipud(EV_fill[n_zeros:n_zeros+n_zeros,k]) - for k in range(SW[1], SE[1]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[SW[0]+1:, k] == 0) # count zeros in South direction - Reflect_fill[-n_zeros:, k] = np.flipud(EV_fill[-n_zeros-n_zeros:-n_zeros, k]) - for k in range(NW[0], SW[0]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[k, 0:NW[1]+1] == 0) # count zeros in North direction - Reflect_fill[k, 0:n_zeros] = np.flipud(EV_fill[k, n_zeros:n_zeros+n_zeros]) - for k in range(NE[0], SE[0]+1): - n_zeros = np.count_nonzero(unw_ifg_interp[k, SE[1]+1:] == 0) # count zeros in North direction - Reflect_fill[k, -n_zeros:] = np.flipud(EV_fill[k, -n_zeros-n_zeros:-n_zeros]) - - Reflect_fill[0:NW[0], 0:NW[1]] = np.flipud(Reflect_fill[NW[0]:NW[0]+NW[0], 0:NW[1]]) # upper left corner area - Reflect_fill[0:NE[0], NE[1]+1:] = np.fliplr(Reflect_fill[0:NE[0], NE[1]+1-(ncol-NE[1]-1):NE[1]+1]) # upper right corner area - Reflect_fill[SW[0]+1:, 0:SW[1]] = np.fliplr(Reflect_fill[SW[0]+1:, SW[1]:SW[1]+SW[1]]) # lower left corner area - Reflect_fill[SE[0]+1:, SE[1]+1:] = np.flipud(Reflect_fill[SE[0]+1-(nrow-SE[0]-1):SE[0]+1, SE[1]+1:]) # lower right corner area + for k in range(NW[1], NE[1] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[0 : NE[0] + 1, k] == 0 + ) # count zeros in North direction + Reflect_fill[0:n_zeros, k] = np.flipud(EV_fill[n_zeros : n_zeros + n_zeros, k]) + for k in range(SW[1], SE[1] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[SW[0] + 1 :, k] == 0 + ) # count zeros in South direction + Reflect_fill[-n_zeros:, k] = np.flipud( + EV_fill[-n_zeros - n_zeros : -n_zeros, k] + ) + for k in range(NW[0], SW[0] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[k, 0 : NW[1] + 1] == 0 + ) # count zeros in North direction + Reflect_fill[k, 0:n_zeros] = np.flipud(EV_fill[k, n_zeros : n_zeros + n_zeros]) + for k in range(NE[0], SE[0] + 1): + n_zeros = np.count_nonzero( + unw_ifg_interp[k, SE[1] + 1 :] == 0 + ) # count zeros in North direction + Reflect_fill[k, -n_zeros:] = np.flipud( + EV_fill[k, -n_zeros - n_zeros : -n_zeros] + ) + + Reflect_fill[0 : NW[0], 0 : NW[1]] = np.flipud( + Reflect_fill[NW[0] : NW[0] + NW[0], 0 : NW[1]] + ) # upper left corner area + Reflect_fill[0 : NE[0], NE[1] + 1 :] = np.fliplr( + Reflect_fill[0 : NE[0], NE[1] + 1 - (ncol - NE[1] - 1) : NE[1] + 1] + ) # upper right corner area + Reflect_fill[SW[0] + 1 :, 0 : SW[1]] = np.fliplr( + Reflect_fill[SW[0] + 1 :, SW[1] : SW[1] + SW[1]] + ) # lower left corner area + Reflect_fill[SE[0] + 1 :, SE[1] + 1 :] = np.flipud( + Reflect_fill[SE[0] + 1 - (nrow - SE[0] - 1) : SE[0] + 1, SE[1] + 1 :] + ) # lower right corner area # 2D filtering with Gaussian kernel - #wavelength_cutoff: float = 50*1e3, - #dx: float = 30, - cutoff_value = 0.5 # 0 < cutoff_value < 1 - sigma_f = 1/wavelength_cutoff/np.sqrt(np.log(1/cutoff_value)) # fc = sqrt(ln(1/cutoff_value))*sigma_f - sigma_x = 1/np.pi/2/sigma_f - sigma = sigma_x/dx - - lowpass_filtered = ndimage.gaussian_filter(Reflect_fill,sigma) - filtered_ifg = unw_ifg - lowpass_filtered*mask_boundary - - return filtered_ifg \ No newline at end of file + # wavelength_cutoff: float = 50*1e3, + # dx: float = 30, + cutoff_value = 0.5 # 0 < cutoff_value < 1 + sigma_f = ( + 1 / wavelength_cutoff / np.sqrt(np.log(1 / cutoff_value)) + ) # fc = sqrt(ln(1/cutoff_value))*sigma_f + sigma_x = 1 / np.pi / 2 / sigma_f + sigma = sigma_x / dx + + lowpass_filtered = ndimage.gaussian_filter(Reflect_fill, sigma) + filtered_ifg = unw_ifg - lowpass_filtered * mask_boundary + + return filtered_ifg From 1cea785bfb14114beb185dabc8693dd7c5444e04 Mon Sep 17 00:00:00 2001 From: Seyeon Jeon Date: Wed, 5 Jun 2024 15:51:12 -0700 Subject: [PATCH 03/23] spelling error --- src/dolphin/filtering.py | 52 ++++++++++++++++++++++++++++++++-------- 1 file changed, 42 insertions(+), 10 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 27f84466..28b99cb3 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -49,16 +49,7 @@ def filtering( mask_boundary = ~(corr == 0).astype("bool") # Ramp plane fitting - Y = unw_ifg[mask] # get data of non NaN & masked pixels - Xdata = np.argwhere(mask) # get indices of non NaN & masked pixels - X = np.c_[np.ones((len(Xdata))), Xdata] - theta = np.dot(np.dot(np.linalg.pinv(np.dot(X.transpose(), X)), X.transpose()), Y) - X1_, X2_ = np.mgrid[:nrow, :ncol] - X_ = np.hstack( - (np.reshape(X1_, (nrow * ncol, 1)), np.reshape(X2_, (nrow * ncol, 1))) - ) - X_ = np.hstack((np.ones((nrow * ncol, 1)), X_)) - plane = np.reshape(np.dot(X_, theta), (nrow, ncol)) + plane = fit_ramp_plane(unw_ifg, mask) # Replace masked out pixels with the ramp plane unw_ifg_interp = np.copy(unw_ifg) @@ -67,6 +58,7 @@ def filtering( # Copy the edge pixels for the boundary area before filling them by reflection EV_fill = np.copy(unw_ifg_interp) + Xdata = np.argwhere(mask) # Get indices of non-NaN & masked pixels NW = Xdata[np.argmin(Xdata[:, 0])] # Get indices of upper left corner pixel SE = Xdata[np.argmax(Xdata[:, 0])] # Get indices of lower right corner pixel SW = Xdata[np.argmin(Xdata[:, 1])] # Get indices of lower left corner pixel @@ -148,3 +140,43 @@ def filtering( filtered_ifg = unw_ifg - lowpass_filtered * mask_boundary return filtered_ifg + + +def fit_ramp_plane(unw_ifg: ArrayLike, mask: ArrayLike) -> np.ndarray: + """Fit a ramp plane to the given data. + + Parameters + ---------- + unw_ifg : ArrayLike + 2D array where the unwrapped interferogram data is stored. + mask : ArrayLike + 2D boolean array indicating the valid (non-NaN) pixels. + + Returns + ------- + np.ndarray + 2D array of the fitted ramp plane. + + """ + # Extract data for non-NaN & masked pixels + Y = unw_ifg[mask] + Xdata = np.argwhere(mask) # Get indices of non-NaN & masked pixels + + # Include the intercept term (bias) in the model + X = np.c_[np.ones((len(Xdata))), Xdata] + + # Compute the parameter vector theta using the least squares solution + theta = np.linalg.pinv(X.T @ X) @ X.T @ Y + + # Prepare grid for the entire image + nrow, ncol = unw_ifg.shape + X1_, X2_ = np.mgrid[:nrow, :ncol] + X_ = np.hstack( + (np.reshape(X1_, (nrow * ncol, 1)), np.reshape(X2_, (nrow * ncol, 1))) + ) + X_ = np.hstack((np.ones((nrow * ncol, 1)), X_)) + + # Compute the fitted plane + plane = np.reshape(X_ @ theta, (nrow, ncol)) + + return plane From d438c470e81a27bb041c7a163a4771d661247b64 Mon Sep 17 00:00:00 2001 From: Seyeon Jeon Date: Wed, 5 Jun 2024 16:14:38 -0700 Subject: [PATCH 04/23] test_filtering.py added --- src/dolphin/filtering.py | 16 ++++++++++++++++ tests/test_filtering.py | 18 ++++++++++++++++++ 2 files changed, 34 insertions(+) create mode 100644 tests/test_filtering.py diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 28b99cb3..83343425 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -68,21 +68,29 @@ def filtering( n_zeros = np.count_nonzero( unw_ifg_interp[0 : NE[0] + 1, k] == 0 ) # count zeros in North direction + if n_zeros == 0: + continue EV_fill[0:n_zeros, k] = EV_fill[n_zeros, k] for k in range(SW[1], SE[1] + 1): n_zeros = np.count_nonzero( unw_ifg_interp[SW[0] + 1 :, k] == 0 ) # count zeros in South direction + if n_zeros == 0: + continue EV_fill[-n_zeros:, k] = EV_fill[-n_zeros - 1, k] for k in range(NW[0], SW[0] + 1): n_zeros = np.count_nonzero( unw_ifg_interp[k, 0 : NW[1] + 1] == 0 ) # count zeros in North direction + if n_zeros == 0: + continue EV_fill[k, 0:n_zeros] = EV_fill[k, n_zeros] for k in range(NE[0], SE[0] + 1): n_zeros = np.count_nonzero( unw_ifg_interp[k, SE[1] + 1 :] == 0 ) # count zeros in North direction + if n_zeros == 0: + continue EV_fill[k, -n_zeros:] = EV_fill[k, -n_zeros - 1] # Fill the boundary area reflecting the pixel values @@ -92,11 +100,15 @@ def filtering( n_zeros = np.count_nonzero( unw_ifg_interp[0 : NE[0] + 1, k] == 0 ) # count zeros in North direction + if n_zeros == 0: + continue Reflect_fill[0:n_zeros, k] = np.flipud(EV_fill[n_zeros : n_zeros + n_zeros, k]) for k in range(SW[1], SE[1] + 1): n_zeros = np.count_nonzero( unw_ifg_interp[SW[0] + 1 :, k] == 0 ) # count zeros in South direction + if n_zeros == 0: + continue Reflect_fill[-n_zeros:, k] = np.flipud( EV_fill[-n_zeros - n_zeros : -n_zeros, k] ) @@ -104,11 +116,15 @@ def filtering( n_zeros = np.count_nonzero( unw_ifg_interp[k, 0 : NW[1] + 1] == 0 ) # count zeros in North direction + if n_zeros == 0: + continue Reflect_fill[k, 0:n_zeros] = np.flipud(EV_fill[k, n_zeros : n_zeros + n_zeros]) for k in range(NE[0], SE[0] + 1): n_zeros = np.count_nonzero( unw_ifg_interp[k, SE[1] + 1 :] == 0 ) # count zeros in North direction + if n_zeros == 0: + continue Reflect_fill[k, -n_zeros:] = np.flipud( EV_fill[k, -n_zeros - n_zeros : -n_zeros] ) diff --git a/tests/test_filtering.py b/tests/test_filtering.py new file mode 100644 index 00000000..28d3a234 --- /dev/null +++ b/tests/test_filtering.py @@ -0,0 +1,18 @@ +import numpy as np + +from dolphin import filtering + + +def test_filtering_unwrapped_phase(): + # Check filtering with ramp phase + y, x = np.ogrid[-3:3:512j, -3:3:512j] + unw_ifg = np.pi * (x + y) + corr = np.ones(unw_ifg.shape, dtype=np.float32) + + # Filtering + filtered_ifg = filtering.filtering(unw_ifg, corr, dx=300) + np.testing.assert_allclose( + filtered_ifg[10:-10, 10:-10], + np.zeros(filtered_ifg[10:-10, 10:-10].shape), + atol=1.0, + ) From fbb885d8f39e1daf9c8f573d5987cf1729750773 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:47:27 -0700 Subject: [PATCH 05/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 83343425..192f3e6e 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -15,7 +15,7 @@ def filtering( Parameters ---------- unw_ifg : np.ndarray, 2D complex array - unwrapped interferogram to interpolate + Unwrapped interferogram phase to filter. corr : 2D float array Array of interferometric correlation from 0 to 1. mask_cutoff: float From b5392857e3f5c0b7a426c789ffe071608df7f4e0 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:48:06 -0700 Subject: [PATCH 06/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 192f3e6e..8b164047 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -16,7 +16,8 @@ def filtering( ---------- unw_ifg : np.ndarray, 2D complex array Unwrapped interferogram phase to filter. - corr : 2D float array + correlation : Arraylike, 2D + Array of interferometric correlation from 0 to 1. Array of interferometric correlation from 0 to 1. mask_cutoff: float Threshold to use on `corr` so that pixels where From fc7594bbb3620ba3c908bf183c5299c8933a79db Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:50:54 -0700 Subject: [PATCH 07/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 8b164047..89bb043c 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -20,8 +20,8 @@ def filtering( Array of interferometric correlation from 0 to 1. Array of interferometric correlation from 0 to 1. mask_cutoff: float - Threshold to use on `corr` so that pixels where - `corr[i, j] > mask_cutoff` are True and the rest are False. + Threshold to use on `correlation` so that pixels where + `correlation[i, j] > mask_cutoff` are used and the rest are ignored. The default is 0.5. wavelength_cutoff: float Spatial wavelength threshold to filter unw_ifg. From 160bb69780a76eb51bdfa049a6de7d2a4f43ce3c Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:51:12 -0700 Subject: [PATCH 08/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 89bb043c..0f62e56b 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -4,11 +4,11 @@ def filtering( - unw_ifg: ArrayLike, - corr: ArrayLike, + unwrapped_phase: ArrayLike, + correlation: ArrayLike, mask_cutoff: float = 0.5, wavelength_cutoff: float = 50 * 1e3, - dx: float = 30, + pixel_spacing: float = 30, ) -> np.ndarray: """Filter out signals with spatial wavelength longer than a threshold. From f720c6f28bb95a599a51d8194cd1349fa01215d9 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:52:04 -0700 Subject: [PATCH 09/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 0f62e56b..6224193a 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -14,7 +14,7 @@ def filtering( Parameters ---------- - unw_ifg : np.ndarray, 2D complex array + unwrapped_phase : np.ndarray, 2D complex array Unwrapped interferogram phase to filter. correlation : Arraylike, 2D Array of interferometric correlation from 0 to 1. From 3ee4134fd202adc747ec764920ddfe48b3e37735 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 15:52:27 -0700 Subject: [PATCH 10/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 6224193a..9347e6dd 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -24,7 +24,7 @@ def filtering( `correlation[i, j] > mask_cutoff` are used and the rest are ignored. The default is 0.5. wavelength_cutoff: float - Spatial wavelength threshold to filter unw_ifg. + Spatial wavelength threshold to filter the unwrapped phase. Signals with wavelength longer than 'wavelength_cutoff' in unw_ifg will be filtered out. The default is 50*1e3 (m). From bd57c6126a41cdcc491428cec31464ac0ddc41dc Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:17:23 -0700 Subject: [PATCH 11/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 9347e6dd..5da5b245 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -25,8 +25,7 @@ def filtering( The default is 0.5. wavelength_cutoff: float Spatial wavelength threshold to filter the unwrapped phase. - Signals with wavelength longer than 'wavelength_cutoff' in unw_ifg will be - filtered out. + Signals with wavelength longer than 'wavelength_cutoff' are filtered out. The default is 50*1e3 (m). dx : float Pixel spatial spacing. Assume same spacing for x, y axes. From 883e7fc2d0aa3572fe0c13077346e2ae8d3af60f Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:17:44 -0700 Subject: [PATCH 12/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 5da5b245..7be55d70 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -27,7 +27,7 @@ def filtering( Spatial wavelength threshold to filter the unwrapped phase. Signals with wavelength longer than 'wavelength_cutoff' are filtered out. The default is 50*1e3 (m). - dx : float + pixel_spacing : float Pixel spatial spacing. Assume same spacing for x, y axes. The default is 30 (m). From ada2911e68365dd0620fa43f68302187f62c897d Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:19:02 -0700 Subject: [PATCH 13/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 7be55d70..e144bb6b 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -39,10 +39,9 @@ def filtering( longer than a threshold. """ - nrow, ncol = corr.shape + nrow, ncol = correlation.shape - # Create Boolean mask for corr > mask_cutoff to be True and the rest to be False - mask = (corr > mask_cutoff).astype("bool") + mask = (correlation > mask_cutoff).astype("bool") # Create Boolean mask for Zero-filled boundary area to be False # and the rest to be True From cfcf3010a1bf4173d53bab052ae6e29ec036450b Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:19:12 -0700 Subject: [PATCH 14/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index e144bb6b..e42f26bd 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -45,7 +45,7 @@ def filtering( # Create Boolean mask for Zero-filled boundary area to be False # and the rest to be True - mask_boundary = ~(corr == 0).astype("bool") + mask_boundary = ~(correlation == 0).astype("bool") # Ramp plane fitting plane = fit_ramp_plane(unw_ifg, mask) From 48b59352b819e30b99c09cdd364e451f4c574b15 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:20:00 -0700 Subject: [PATCH 15/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index e42f26bd..3d98efe4 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -47,7 +47,6 @@ def filtering( # and the rest to be True mask_boundary = ~(correlation == 0).astype("bool") - # Ramp plane fitting plane = fit_ramp_plane(unw_ifg, mask) # Replace masked out pixels with the ramp plane From 10549fd0d1a31c55c3199e637ddf9a1a66463423 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:20:13 -0700 Subject: [PATCH 16/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 3d98efe4..16e35c60 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -47,10 +47,10 @@ def filtering( # and the rest to be True mask_boundary = ~(correlation == 0).astype("bool") - plane = fit_ramp_plane(unw_ifg, mask) + plane = fit_ramp_plane(unwrapped_phase, mask) # Replace masked out pixels with the ramp plane - unw_ifg_interp = np.copy(unw_ifg) + unw_ifg_interp = np.copy(unwrapped_phase) unw_ifg_interp[~mask * mask_boundary] = plane[~mask * mask_boundary] # Copy the edge pixels for the boundary area before filling them by reflection From a8ecdf0dd156458f16af8ad1af61e6f2241c9b82 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:20:38 -0700 Subject: [PATCH 17/23] Update src/dolphin/filtering.py Co-authored-by: Scott Staniewicz --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 16e35c60..066ea6c6 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -3,7 +3,7 @@ from scipy import ndimage -def filtering( +def filter_long_wavelength( unwrapped_phase: ArrayLike, correlation: ArrayLike, mask_cutoff: float = 0.5, From 77ed07832364be2ae6fe2c79e0bbe0486ce79564 Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:21:03 -0700 Subject: [PATCH 18/23] Update tests/test_filtering.py Co-authored-by: Scott Staniewicz --- tests/test_filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_filtering.py b/tests/test_filtering.py index 28d3a234..0ad244b3 100644 --- a/tests/test_filtering.py +++ b/tests/test_filtering.py @@ -3,7 +3,7 @@ from dolphin import filtering -def test_filtering_unwrapped_phase(): +def test_filter_long_wavelegnth(): # Check filtering with ramp phase y, x = np.ogrid[-3:3:512j, -3:3:512j] unw_ifg = np.pi * (x + y) From e39225df7349e266d8c65bb2396905069f5d005b Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 11 Jun 2024 16:22:43 -0700 Subject: [PATCH 19/23] Update tests/test_filtering.py Co-authored-by: Scott Staniewicz --- tests/test_filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_filtering.py b/tests/test_filtering.py index 0ad244b3..badffe89 100644 --- a/tests/test_filtering.py +++ b/tests/test_filtering.py @@ -10,7 +10,7 @@ def test_filter_long_wavelegnth(): corr = np.ones(unw_ifg.shape, dtype=np.float32) # Filtering - filtered_ifg = filtering.filtering(unw_ifg, corr, dx=300) + filtered_ifg = filtering.filter_long_wavelength(unw_ifg, corr, pixel_spacing=300) np.testing.assert_allclose( filtered_ifg[10:-10, 10:-10], np.zeros(filtered_ifg[10:-10, 10:-10].shape), From aeaf548b98f2fdea02474144bf72f4fc597b35f9 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 11 Jun 2024 20:07:18 -0400 Subject: [PATCH 20/23] Apply suggestions from code review Co-authored-by: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> --- src/dolphin/filtering.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 066ea6c6..a76c36b1 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -148,10 +148,10 @@ def filter_long_wavelength( 1 / wavelength_cutoff / np.sqrt(np.log(1 / cutoff_value)) ) # fc = sqrt(ln(1/cutoff_value))*sigma_f sigma_x = 1 / np.pi / 2 / sigma_f - sigma = sigma_x / dx + sigma = sigma_x / pixel_spacing lowpass_filtered = ndimage.gaussian_filter(Reflect_fill, sigma) - filtered_ifg = unw_ifg - lowpass_filtered * mask_boundary + filtered_ifg = unwrapped_phase - lowpass_filtered * mask_boundary return filtered_ifg From 013d7ec360b5116bf9d9f5304692e9ca4ee217fd Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 11 Jun 2024 20:08:27 -0400 Subject: [PATCH 21/23] Update src/dolphin/filtering.py --- src/dolphin/filtering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index a76c36b1..ecfd69e3 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -16,7 +16,7 @@ def filter_long_wavelength( ---------- unwrapped_phase : np.ndarray, 2D complex array Unwrapped interferogram phase to filter. - correlation : Arraylike, 2D + correlation : Arraylike, 2D Array of interferometric correlation from 0 to 1. Array of interferometric correlation from 0 to 1. mask_cutoff: float From 85701d0707bb71e4048e44c763b7bdcb59c9b6db Mon Sep 17 00:00:00 2001 From: seyeonjeon <133063223+seyeonjeon@users.noreply.github.com> Date: Tue, 29 Oct 2024 23:33:25 -0700 Subject: [PATCH 22/23] Update filtering.py reflect fill using np.pad --- src/dolphin/filtering.py | 348 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 348 insertions(+) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index e69de29b..56242082 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -0,0 +1,348 @@ +import multiprocessing as mp +from concurrent.futures import ProcessPoolExecutor +from itertools import repeat +from pathlib import Path +from typing import Sequence + +import numpy as np +from numpy.typing import ArrayLike, NDArray +from scipy import fft, ndimage + + +def filter_long_wavelength( + unwrapped_phase: ArrayLike, + bad_pixel_mask: ArrayLike, + wavelength_cutoff: float = 50 * 1e3, + pixel_spacing: float = 30, + workers: int = 1, +) -> np.ndarray: + """Filter out signals with spatial wavelength longer than a threshold. + + Parameters + ---------- + unwrapped_phase : ArrayLike + Unwrapped interferogram phase to filter. + bad_pixel_mask : ArrayLike + Boolean array with same shape as `unwrapped_phase` where `True` indicates a + pixel to ignore during ramp fitting + wavelength_cutoff : float + Spatial wavelength threshold to filter the unwrapped phase. + Signals with wavelength longer than 'wavelength_cutoff' are filtered out. + The default is 50*1e3 (m). + pixel_spacing : float + Pixel spatial spacing. Assume same spacing for x, y axes. + The default is 30 (m). + workers : int + Number of `fft` workers to use for `scipy.fft.fft2`. + Default is 1. + + Returns + ------- + filtered_ifg : 2D complex array + filtered interferogram that does not contain signals with spatial wavelength + longer than a threshold. + + Raises + ------ + ValueError + If wavelength_cutoff too large for image size/pixel spacing. + + """ + good_pixel_mask = ~bad_pixel_mask + + rows, cols = unwrapped_phase.shape + unw0 = np.nan_to_num(unwrapped_phase) + # Take either nan or 0 pixels in `unwrapped_phase` to be nodata + nodata_mask = unw0 == 0 + in_bounds_pixels = ~nodata_mask + + total_valid_mask = in_bounds_pixels & good_pixel_mask + + plane = fit_ramp_plane(unw0, total_valid_mask) + # Remove the plane, setting to 0 where we had no data for the plane fit: + unw_ifg_interp = np.where((~nodata_mask & good_pixel_mask), unw0, plane) + + # Fill the boundary area by reflecting pixel values + reflect_fill = _fill_boundary_area(unw_ifg_interp, in_bounds_pixels) + + # Find the filter `sigma` which gives the correct cutoff in meters + sigma = _compute_filter_sigma(wavelength_cutoff, pixel_spacing, cutoff_value=0.5) + + if sigma > unw0.shape[0] or sigma > unw0.shape[0]: + msg = f"{wavelength_cutoff = } too large for image." + msg += f"Shape = {(rows, cols)}, and {pixel_spacing = }" + raise ValueError(msg) + # Pad the array with edge values + # The padding extends further than the default "radius = 2*sigma + 1", + # which given specified in `gaussian_filter` + # https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html#scipy.ndimage.gaussian_filter + pad_rows = pad_cols = int(3 * sigma) + # See here for illustration of `mode="reflect"` + # https://scikit-image.org/docs/stable/auto_examples/transform/plot_edge_modes.html#interpolation-edge-modes + padded = np.pad( + reflect_fill, ((pad_rows, pad_rows), (pad_cols, pad_cols)), mode="reflect" + ) + + # Apply Gaussian filter + result = fft.fft2(padded, workers=workers) + result = ndimage.fourier_gaussian(result, sigma=sigma) + # Make sure to only take the real part (ifft returns complex) + result = fft.ifft2(result, workers=workers).real.astype(unwrapped_phase.dtype) + + # Crop back to original size + lowpass_filtered = result[pad_rows:-pad_rows, pad_cols:-pad_cols] + + filtered_ifg = unw_ifg_interp - lowpass_filtered * in_bounds_pixels + return np.where(in_bounds_pixels, filtered_ifg, 0) + + +def _fill_boundary_area(unw_ifg_interp: np.ndarray, mask: np.ndarray) -> np.ndarray: + """Fill the boundary area by reflecting pixel values.""" + # Rotate the data to align the frame + rt_unw_ifg = ndimage.rotate(unw_ifg_interp, 8.6, reshape=False) + + rt_pad = 700 # Apply padding for rotating back + # Fill boundary area using np.pad for each horizontal step and + reflect_filled = np.pad( + rt_unw_ifg[846:6059, :], + ((846, rt_unw_ifg.shape[0] - 6059), (0, 0)), + mode="reflect", + ) + reflect_filled = np.pad( + reflect_filled[618:6241, :], + ((618, rt_unw_ifg.shape[0] - 6241), (0, 0)), + mode="reflect", + ) + reflect_filled = np.pad( + reflect_filled[399:6447, :], + ((399 + rt_pad, rt_unw_ifg.shape[0] - 6447 + rt_pad), (0, 0)), + mode="reflect", + ) + # Fill vertical boundary area using np.pad + reflect_filled = np.pad( + reflect_filled[:, 591:8861], + ((0, 0), (591 + rt_pad, rt_unw_ifg.shape[1] - 8861 + rt_pad)), + mode="reflect", + ) + + # Rotate back to original angle + reflect_filled = ndimage.rotate(reflect_filled, -8.6, reshape=False) + reflect_filled = reflect_filled[rt_pad:-rt_pad, rt_pad:-rt_pad] + # Copy original in-bound pixels + reflect_filled[mask] = unw_ifg_interp[mask] + + return reflect_filled + + +def _compute_filter_sigma( + wavelength_cutoff: float, pixel_spacing: float, cutoff_value: float = 0.5 +) -> float: + sigma_f = 1 / wavelength_cutoff / np.sqrt(np.log(1 / cutoff_value)) + sigma_x = 1 / np.pi / 2 / sigma_f + sigma = sigma_x / pixel_spacing + return sigma + + +def fit_ramp_plane(unw_ifg: ArrayLike, mask: ArrayLike) -> np.ndarray: + """Fit a ramp plane to the given data. + + Parameters + ---------- + unw_ifg : ArrayLike + 2D array where the unwrapped interferogram data is stored. + mask : ArrayLike + 2D boolean array indicating the valid (non-NaN) pixels. + + Returns + ------- + np.ndarray + 2D array of the fitted ramp plane. + + """ + # Extract data for non-NaN & masked pixels + Y = unw_ifg[mask] + Xdata = np.argwhere(mask) # Get indices of non-NaN & masked pixels + + # Include the intercept term (bias) in the model + X = np.c_[np.ones((len(Xdata))), Xdata] + + # Compute the parameter vector theta using the least squares solution + theta = np.linalg.pinv(X.T @ X) @ X.T @ Y + + # Prepare grid for the entire image + nrow, ncol = unw_ifg.shape + X1_, X2_ = np.mgrid[:nrow, :ncol] + X_ = np.hstack( + (np.reshape(X1_, (nrow * ncol, 1)), np.reshape(X2_, (nrow * ncol, 1))) + ) + X_ = np.hstack((np.ones((nrow * ncol, 1)), X_)) + + # Compute the fitted plane + plane = np.reshape(X_ @ theta, (nrow, ncol)) + + return plane + + +def filter_rasters( + unw_filenames: list[Path], + cor_filenames: list[Path] | None = None, + conncomp_filenames: list[Path] | None = None, + temporal_coherence_filename: Path | None = None, + wavelength_cutoff: float = 50_000, + correlation_cutoff: float = 0.5, + output_dir: Path | None = None, + max_workers: int = 4, +) -> list[Path]: + """Filter a list of unwrapped interferogram files using a long-wavelength filter. + + Remove long-wavelength components from each unwrapped interferogram. + It can optionally use temporal coherence, correlation, and connected component + information for masking. + + Parameters + ---------- + unw_filenames : list[Path] + List of paths to unwrapped interferogram files to be filtered. + cor_filenames : list[Path] | None + List of paths to correlation files + Passing None skips filtering on correlation. + conncomp_filenames : list[Path] | None + List of paths to connected component files, filters any 0 labeled pixels. + Passing None skips filtering on connected component labels. + temporal_coherence_filename : Path | None + Path to the temporal coherence file for masking. + Passing None skips filtering on temporal coherence. + wavelength_cutoff : float, optional + Spatial wavelength cutoff (in meters) for the filter. Default is 50,000 meters. + correlation_cutoff : float, optional + Threshold of correlation (if passing `cor_filenames`) to use to ignore pixels + during filtering. + output_dir : Path | None, optional + Directory to save the filtered results. + If None, saves in the same location as inputs with .filt.tif extension. + max_workers : int, optional + Number of parallel images to process. Default is 4. + + Returns + ------- + list[Path] + Output filtered rasters. + + Notes + ----- + - If temporal_coherence_filename is provided, pixels with coherence < 0.5 are masked + + """ + from dolphin import io + + bad_pixel_mask = np.zeros( + io.get_raster_xysize(unw_filenames[0])[::-1], dtype="bool" + ) + if temporal_coherence_filename: + bad_pixel_mask = bad_pixel_mask | ( + io.load_gdal(temporal_coherence_filename) < 0.5 + ) + + if output_dir is None: + assert unw_filenames + output_dir = unw_filenames[0].parent + output_dir.mkdir(exist_ok=True) + ctx = mp.get_context("spawn") + + with ProcessPoolExecutor(max_workers, mp_context=ctx) as pool: + return list( + pool.map( + _filter_and_save, + unw_filenames, + cor_filenames or repeat(None), + conncomp_filenames or repeat(None), + repeat(output_dir), + repeat(wavelength_cutoff), + repeat(bad_pixel_mask), + repeat(correlation_cutoff), + ) + ) + + +def _filter_and_save( + unw_filename: Path, + cor_path: Path | None, + conncomp_path: Path | None, + output_dir: Path, + wavelength_cutoff: float, + bad_pixel_mask: NDArray[np.bool_], + correlation_cutoff: float = 0.5, +) -> Path: + """Filter one interferogram (wrapper for multiprocessing).""" + from dolphin import io + from dolphin._overviews import Resampling, create_image_overviews + + # Average for the pixel spacing for filtering + _, x_res, _, _, _, y_res = io.get_raster_gt(unw_filename) + pixel_spacing = (abs(x_res) + abs(y_res)) / 2 + + if cor_path is not None: + bad_pixel_mask |= io.load_gdal(cor_path) < correlation_cutoff + if conncomp_path is not None: + bad_pixel_mask |= io.load_gdal(conncomp_path, masked=True).astype(bool) == 0 + + unw = io.load_gdal(unw_filename) + filt_arr = filter_long_wavelength( + unwrapped_phase=unw, + wavelength_cutoff=wavelength_cutoff, + bad_pixel_mask=bad_pixel_mask, + pixel_spacing=pixel_spacing, + workers=1, + ) + io.round_mantissa(filt_arr, keep_bits=9) + output_name = output_dir / Path(unw_filename).with_suffix(".filt.tif").name + io.write_arr(arr=filt_arr, like_filename=unw_filename, output_name=output_name) + + create_image_overviews(output_name, resampling=Resampling.AVERAGE) + + return output_name + + +def gaussian_filter_nan( + image: ArrayLike, sigma: float | Sequence[float], mode="constant", **kwargs +) -> np.ndarray: + """Apply a gaussian filter to an image with NaNs (avoiding all nans). + + The scipy.ndimage `gaussian_filter` will make the output all NaNs if + any of the pixels in the input that touches the kernel is NaN + + Source: + https://stackoverflow.com/a/36307291 + + Parameters + ---------- + image : ndarray + Image with nans to filter + sigma : float + Size of filter kernel. passed into `gaussian_filter` + mode : str, default = "constant" + Boundary mode for `[scipy.ndimage.gaussian_filter][]` + **kwargs : Any + Passed into `[scipy.ndimage.gaussian_filter][]` + + Returns + ------- + ndarray + Filtered version of `image`. + + """ + from scipy.ndimage import gaussian_filter + + if np.sum(np.isnan(image)) == 0: + return gaussian_filter(image, sigma=sigma, mode=mode, **kwargs) + + V = image.copy() + nan_idxs = np.isnan(image) + V[nan_idxs] = 0 + V_filt = gaussian_filter(V, sigma, **kwargs) + + W = np.ones(image.shape) + W[nan_idxs] = 0 + W_filt = gaussian_filter(W, sigma, **kwargs) + + return V_filt / W_filt From de9276367f432f388dd154e2493be877706a8a83 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 22 Nov 2024 15:12:54 -0500 Subject: [PATCH 23/23] change CRLF to LF --- src/dolphin/filtering.py | 696 +++++++++++++++++++-------------------- 1 file changed, 348 insertions(+), 348 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 56242082..e23b19af 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -1,348 +1,348 @@ -import multiprocessing as mp -from concurrent.futures import ProcessPoolExecutor -from itertools import repeat -from pathlib import Path -from typing import Sequence - -import numpy as np -from numpy.typing import ArrayLike, NDArray -from scipy import fft, ndimage - - -def filter_long_wavelength( - unwrapped_phase: ArrayLike, - bad_pixel_mask: ArrayLike, - wavelength_cutoff: float = 50 * 1e3, - pixel_spacing: float = 30, - workers: int = 1, -) -> np.ndarray: - """Filter out signals with spatial wavelength longer than a threshold. - - Parameters - ---------- - unwrapped_phase : ArrayLike - Unwrapped interferogram phase to filter. - bad_pixel_mask : ArrayLike - Boolean array with same shape as `unwrapped_phase` where `True` indicates a - pixel to ignore during ramp fitting - wavelength_cutoff : float - Spatial wavelength threshold to filter the unwrapped phase. - Signals with wavelength longer than 'wavelength_cutoff' are filtered out. - The default is 50*1e3 (m). - pixel_spacing : float - Pixel spatial spacing. Assume same spacing for x, y axes. - The default is 30 (m). - workers : int - Number of `fft` workers to use for `scipy.fft.fft2`. - Default is 1. - - Returns - ------- - filtered_ifg : 2D complex array - filtered interferogram that does not contain signals with spatial wavelength - longer than a threshold. - - Raises - ------ - ValueError - If wavelength_cutoff too large for image size/pixel spacing. - - """ - good_pixel_mask = ~bad_pixel_mask - - rows, cols = unwrapped_phase.shape - unw0 = np.nan_to_num(unwrapped_phase) - # Take either nan or 0 pixels in `unwrapped_phase` to be nodata - nodata_mask = unw0 == 0 - in_bounds_pixels = ~nodata_mask - - total_valid_mask = in_bounds_pixels & good_pixel_mask - - plane = fit_ramp_plane(unw0, total_valid_mask) - # Remove the plane, setting to 0 where we had no data for the plane fit: - unw_ifg_interp = np.where((~nodata_mask & good_pixel_mask), unw0, plane) - - # Fill the boundary area by reflecting pixel values - reflect_fill = _fill_boundary_area(unw_ifg_interp, in_bounds_pixels) - - # Find the filter `sigma` which gives the correct cutoff in meters - sigma = _compute_filter_sigma(wavelength_cutoff, pixel_spacing, cutoff_value=0.5) - - if sigma > unw0.shape[0] or sigma > unw0.shape[0]: - msg = f"{wavelength_cutoff = } too large for image." - msg += f"Shape = {(rows, cols)}, and {pixel_spacing = }" - raise ValueError(msg) - # Pad the array with edge values - # The padding extends further than the default "radius = 2*sigma + 1", - # which given specified in `gaussian_filter` - # https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html#scipy.ndimage.gaussian_filter - pad_rows = pad_cols = int(3 * sigma) - # See here for illustration of `mode="reflect"` - # https://scikit-image.org/docs/stable/auto_examples/transform/plot_edge_modes.html#interpolation-edge-modes - padded = np.pad( - reflect_fill, ((pad_rows, pad_rows), (pad_cols, pad_cols)), mode="reflect" - ) - - # Apply Gaussian filter - result = fft.fft2(padded, workers=workers) - result = ndimage.fourier_gaussian(result, sigma=sigma) - # Make sure to only take the real part (ifft returns complex) - result = fft.ifft2(result, workers=workers).real.astype(unwrapped_phase.dtype) - - # Crop back to original size - lowpass_filtered = result[pad_rows:-pad_rows, pad_cols:-pad_cols] - - filtered_ifg = unw_ifg_interp - lowpass_filtered * in_bounds_pixels - return np.where(in_bounds_pixels, filtered_ifg, 0) - - -def _fill_boundary_area(unw_ifg_interp: np.ndarray, mask: np.ndarray) -> np.ndarray: - """Fill the boundary area by reflecting pixel values.""" - # Rotate the data to align the frame - rt_unw_ifg = ndimage.rotate(unw_ifg_interp, 8.6, reshape=False) - - rt_pad = 700 # Apply padding for rotating back - # Fill boundary area using np.pad for each horizontal step and - reflect_filled = np.pad( - rt_unw_ifg[846:6059, :], - ((846, rt_unw_ifg.shape[0] - 6059), (0, 0)), - mode="reflect", - ) - reflect_filled = np.pad( - reflect_filled[618:6241, :], - ((618, rt_unw_ifg.shape[0] - 6241), (0, 0)), - mode="reflect", - ) - reflect_filled = np.pad( - reflect_filled[399:6447, :], - ((399 + rt_pad, rt_unw_ifg.shape[0] - 6447 + rt_pad), (0, 0)), - mode="reflect", - ) - # Fill vertical boundary area using np.pad - reflect_filled = np.pad( - reflect_filled[:, 591:8861], - ((0, 0), (591 + rt_pad, rt_unw_ifg.shape[1] - 8861 + rt_pad)), - mode="reflect", - ) - - # Rotate back to original angle - reflect_filled = ndimage.rotate(reflect_filled, -8.6, reshape=False) - reflect_filled = reflect_filled[rt_pad:-rt_pad, rt_pad:-rt_pad] - # Copy original in-bound pixels - reflect_filled[mask] = unw_ifg_interp[mask] - - return reflect_filled - - -def _compute_filter_sigma( - wavelength_cutoff: float, pixel_spacing: float, cutoff_value: float = 0.5 -) -> float: - sigma_f = 1 / wavelength_cutoff / np.sqrt(np.log(1 / cutoff_value)) - sigma_x = 1 / np.pi / 2 / sigma_f - sigma = sigma_x / pixel_spacing - return sigma - - -def fit_ramp_plane(unw_ifg: ArrayLike, mask: ArrayLike) -> np.ndarray: - """Fit a ramp plane to the given data. - - Parameters - ---------- - unw_ifg : ArrayLike - 2D array where the unwrapped interferogram data is stored. - mask : ArrayLike - 2D boolean array indicating the valid (non-NaN) pixels. - - Returns - ------- - np.ndarray - 2D array of the fitted ramp plane. - - """ - # Extract data for non-NaN & masked pixels - Y = unw_ifg[mask] - Xdata = np.argwhere(mask) # Get indices of non-NaN & masked pixels - - # Include the intercept term (bias) in the model - X = np.c_[np.ones((len(Xdata))), Xdata] - - # Compute the parameter vector theta using the least squares solution - theta = np.linalg.pinv(X.T @ X) @ X.T @ Y - - # Prepare grid for the entire image - nrow, ncol = unw_ifg.shape - X1_, X2_ = np.mgrid[:nrow, :ncol] - X_ = np.hstack( - (np.reshape(X1_, (nrow * ncol, 1)), np.reshape(X2_, (nrow * ncol, 1))) - ) - X_ = np.hstack((np.ones((nrow * ncol, 1)), X_)) - - # Compute the fitted plane - plane = np.reshape(X_ @ theta, (nrow, ncol)) - - return plane - - -def filter_rasters( - unw_filenames: list[Path], - cor_filenames: list[Path] | None = None, - conncomp_filenames: list[Path] | None = None, - temporal_coherence_filename: Path | None = None, - wavelength_cutoff: float = 50_000, - correlation_cutoff: float = 0.5, - output_dir: Path | None = None, - max_workers: int = 4, -) -> list[Path]: - """Filter a list of unwrapped interferogram files using a long-wavelength filter. - - Remove long-wavelength components from each unwrapped interferogram. - It can optionally use temporal coherence, correlation, and connected component - information for masking. - - Parameters - ---------- - unw_filenames : list[Path] - List of paths to unwrapped interferogram files to be filtered. - cor_filenames : list[Path] | None - List of paths to correlation files - Passing None skips filtering on correlation. - conncomp_filenames : list[Path] | None - List of paths to connected component files, filters any 0 labeled pixels. - Passing None skips filtering on connected component labels. - temporal_coherence_filename : Path | None - Path to the temporal coherence file for masking. - Passing None skips filtering on temporal coherence. - wavelength_cutoff : float, optional - Spatial wavelength cutoff (in meters) for the filter. Default is 50,000 meters. - correlation_cutoff : float, optional - Threshold of correlation (if passing `cor_filenames`) to use to ignore pixels - during filtering. - output_dir : Path | None, optional - Directory to save the filtered results. - If None, saves in the same location as inputs with .filt.tif extension. - max_workers : int, optional - Number of parallel images to process. Default is 4. - - Returns - ------- - list[Path] - Output filtered rasters. - - Notes - ----- - - If temporal_coherence_filename is provided, pixels with coherence < 0.5 are masked - - """ - from dolphin import io - - bad_pixel_mask = np.zeros( - io.get_raster_xysize(unw_filenames[0])[::-1], dtype="bool" - ) - if temporal_coherence_filename: - bad_pixel_mask = bad_pixel_mask | ( - io.load_gdal(temporal_coherence_filename) < 0.5 - ) - - if output_dir is None: - assert unw_filenames - output_dir = unw_filenames[0].parent - output_dir.mkdir(exist_ok=True) - ctx = mp.get_context("spawn") - - with ProcessPoolExecutor(max_workers, mp_context=ctx) as pool: - return list( - pool.map( - _filter_and_save, - unw_filenames, - cor_filenames or repeat(None), - conncomp_filenames or repeat(None), - repeat(output_dir), - repeat(wavelength_cutoff), - repeat(bad_pixel_mask), - repeat(correlation_cutoff), - ) - ) - - -def _filter_and_save( - unw_filename: Path, - cor_path: Path | None, - conncomp_path: Path | None, - output_dir: Path, - wavelength_cutoff: float, - bad_pixel_mask: NDArray[np.bool_], - correlation_cutoff: float = 0.5, -) -> Path: - """Filter one interferogram (wrapper for multiprocessing).""" - from dolphin import io - from dolphin._overviews import Resampling, create_image_overviews - - # Average for the pixel spacing for filtering - _, x_res, _, _, _, y_res = io.get_raster_gt(unw_filename) - pixel_spacing = (abs(x_res) + abs(y_res)) / 2 - - if cor_path is not None: - bad_pixel_mask |= io.load_gdal(cor_path) < correlation_cutoff - if conncomp_path is not None: - bad_pixel_mask |= io.load_gdal(conncomp_path, masked=True).astype(bool) == 0 - - unw = io.load_gdal(unw_filename) - filt_arr = filter_long_wavelength( - unwrapped_phase=unw, - wavelength_cutoff=wavelength_cutoff, - bad_pixel_mask=bad_pixel_mask, - pixel_spacing=pixel_spacing, - workers=1, - ) - io.round_mantissa(filt_arr, keep_bits=9) - output_name = output_dir / Path(unw_filename).with_suffix(".filt.tif").name - io.write_arr(arr=filt_arr, like_filename=unw_filename, output_name=output_name) - - create_image_overviews(output_name, resampling=Resampling.AVERAGE) - - return output_name - - -def gaussian_filter_nan( - image: ArrayLike, sigma: float | Sequence[float], mode="constant", **kwargs -) -> np.ndarray: - """Apply a gaussian filter to an image with NaNs (avoiding all nans). - - The scipy.ndimage `gaussian_filter` will make the output all NaNs if - any of the pixels in the input that touches the kernel is NaN - - Source: - https://stackoverflow.com/a/36307291 - - Parameters - ---------- - image : ndarray - Image with nans to filter - sigma : float - Size of filter kernel. passed into `gaussian_filter` - mode : str, default = "constant" - Boundary mode for `[scipy.ndimage.gaussian_filter][]` - **kwargs : Any - Passed into `[scipy.ndimage.gaussian_filter][]` - - Returns - ------- - ndarray - Filtered version of `image`. - - """ - from scipy.ndimage import gaussian_filter - - if np.sum(np.isnan(image)) == 0: - return gaussian_filter(image, sigma=sigma, mode=mode, **kwargs) - - V = image.copy() - nan_idxs = np.isnan(image) - V[nan_idxs] = 0 - V_filt = gaussian_filter(V, sigma, **kwargs) - - W = np.ones(image.shape) - W[nan_idxs] = 0 - W_filt = gaussian_filter(W, sigma, **kwargs) - - return V_filt / W_filt +import multiprocessing as mp +from concurrent.futures import ProcessPoolExecutor +from itertools import repeat +from pathlib import Path +from typing import Sequence + +import numpy as np +from numpy.typing import ArrayLike, NDArray +from scipy import fft, ndimage + + +def filter_long_wavelength( + unwrapped_phase: ArrayLike, + bad_pixel_mask: ArrayLike, + wavelength_cutoff: float = 50 * 1e3, + pixel_spacing: float = 30, + workers: int = 1, +) -> np.ndarray: + """Filter out signals with spatial wavelength longer than a threshold. + + Parameters + ---------- + unwrapped_phase : ArrayLike + Unwrapped interferogram phase to filter. + bad_pixel_mask : ArrayLike + Boolean array with same shape as `unwrapped_phase` where `True` indicates a + pixel to ignore during ramp fitting + wavelength_cutoff : float + Spatial wavelength threshold to filter the unwrapped phase. + Signals with wavelength longer than 'wavelength_cutoff' are filtered out. + The default is 50*1e3 (m). + pixel_spacing : float + Pixel spatial spacing. Assume same spacing for x, y axes. + The default is 30 (m). + workers : int + Number of `fft` workers to use for `scipy.fft.fft2`. + Default is 1. + + Returns + ------- + filtered_ifg : 2D complex array + filtered interferogram that does not contain signals with spatial wavelength + longer than a threshold. + + Raises + ------ + ValueError + If wavelength_cutoff too large for image size/pixel spacing. + + """ + good_pixel_mask = ~bad_pixel_mask + + rows, cols = unwrapped_phase.shape + unw0 = np.nan_to_num(unwrapped_phase) + # Take either nan or 0 pixels in `unwrapped_phase` to be nodata + nodata_mask = unw0 == 0 + in_bounds_pixels = ~nodata_mask + + total_valid_mask = in_bounds_pixels & good_pixel_mask + + plane = fit_ramp_plane(unw0, total_valid_mask) + # Remove the plane, setting to 0 where we had no data for the plane fit: + unw_ifg_interp = np.where((~nodata_mask & good_pixel_mask), unw0, plane) + + # Fill the boundary area by reflecting pixel values + reflect_fill = _fill_boundary_area(unw_ifg_interp, in_bounds_pixels) + + # Find the filter `sigma` which gives the correct cutoff in meters + sigma = _compute_filter_sigma(wavelength_cutoff, pixel_spacing, cutoff_value=0.5) + + if sigma > unw0.shape[0] or sigma > unw0.shape[0]: + msg = f"{wavelength_cutoff = } too large for image." + msg += f"Shape = {(rows, cols)}, and {pixel_spacing = }" + raise ValueError(msg) + # Pad the array with edge values + # The padding extends further than the default "radius = 2*sigma + 1", + # which given specified in `gaussian_filter` + # https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html#scipy.ndimage.gaussian_filter + pad_rows = pad_cols = int(3 * sigma) + # See here for illustration of `mode="reflect"` + # https://scikit-image.org/docs/stable/auto_examples/transform/plot_edge_modes.html#interpolation-edge-modes + padded = np.pad( + reflect_fill, ((pad_rows, pad_rows), (pad_cols, pad_cols)), mode="reflect" + ) + + # Apply Gaussian filter + result = fft.fft2(padded, workers=workers) + result = ndimage.fourier_gaussian(result, sigma=sigma) + # Make sure to only take the real part (ifft returns complex) + result = fft.ifft2(result, workers=workers).real.astype(unwrapped_phase.dtype) + + # Crop back to original size + lowpass_filtered = result[pad_rows:-pad_rows, pad_cols:-pad_cols] + + filtered_ifg = unw_ifg_interp - lowpass_filtered * in_bounds_pixels + return np.where(in_bounds_pixels, filtered_ifg, 0) + + +def _fill_boundary_area(unw_ifg_interp: np.ndarray, mask: np.ndarray) -> np.ndarray: + """Fill the boundary area by reflecting pixel values.""" + # Rotate the data to align the frame + rt_unw_ifg = ndimage.rotate(unw_ifg_interp, 8.6, reshape=False) + + rt_pad = 700 # Apply padding for rotating back + # Fill boundary area using np.pad for each horizontal step and + reflect_filled = np.pad( + rt_unw_ifg[846:6059, :], + ((846, rt_unw_ifg.shape[0] - 6059), (0, 0)), + mode="reflect", + ) + reflect_filled = np.pad( + reflect_filled[618:6241, :], + ((618, rt_unw_ifg.shape[0] - 6241), (0, 0)), + mode="reflect", + ) + reflect_filled = np.pad( + reflect_filled[399:6447, :], + ((399 + rt_pad, rt_unw_ifg.shape[0] - 6447 + rt_pad), (0, 0)), + mode="reflect", + ) + # Fill vertical boundary area using np.pad + reflect_filled = np.pad( + reflect_filled[:, 591:8861], + ((0, 0), (591 + rt_pad, rt_unw_ifg.shape[1] - 8861 + rt_pad)), + mode="reflect", + ) + + # Rotate back to original angle + reflect_filled = ndimage.rotate(reflect_filled, -8.6, reshape=False) + reflect_filled = reflect_filled[rt_pad:-rt_pad, rt_pad:-rt_pad] + # Copy original in-bound pixels + reflect_filled[mask] = unw_ifg_interp[mask] + + return reflect_filled + + +def _compute_filter_sigma( + wavelength_cutoff: float, pixel_spacing: float, cutoff_value: float = 0.5 +) -> float: + sigma_f = 1 / wavelength_cutoff / np.sqrt(np.log(1 / cutoff_value)) + sigma_x = 1 / np.pi / 2 / sigma_f + sigma = sigma_x / pixel_spacing + return sigma + + +def fit_ramp_plane(unw_ifg: ArrayLike, mask: ArrayLike) -> np.ndarray: + """Fit a ramp plane to the given data. + + Parameters + ---------- + unw_ifg : ArrayLike + 2D array where the unwrapped interferogram data is stored. + mask : ArrayLike + 2D boolean array indicating the valid (non-NaN) pixels. + + Returns + ------- + np.ndarray + 2D array of the fitted ramp plane. + + """ + # Extract data for non-NaN & masked pixels + Y = unw_ifg[mask] + Xdata = np.argwhere(mask) # Get indices of non-NaN & masked pixels + + # Include the intercept term (bias) in the model + X = np.c_[np.ones((len(Xdata))), Xdata] + + # Compute the parameter vector theta using the least squares solution + theta = np.linalg.pinv(X.T @ X) @ X.T @ Y + + # Prepare grid for the entire image + nrow, ncol = unw_ifg.shape + X1_, X2_ = np.mgrid[:nrow, :ncol] + X_ = np.hstack( + (np.reshape(X1_, (nrow * ncol, 1)), np.reshape(X2_, (nrow * ncol, 1))) + ) + X_ = np.hstack((np.ones((nrow * ncol, 1)), X_)) + + # Compute the fitted plane + plane = np.reshape(X_ @ theta, (nrow, ncol)) + + return plane + + +def filter_rasters( + unw_filenames: list[Path], + cor_filenames: list[Path] | None = None, + conncomp_filenames: list[Path] | None = None, + temporal_coherence_filename: Path | None = None, + wavelength_cutoff: float = 50_000, + correlation_cutoff: float = 0.5, + output_dir: Path | None = None, + max_workers: int = 4, +) -> list[Path]: + """Filter a list of unwrapped interferogram files using a long-wavelength filter. + + Remove long-wavelength components from each unwrapped interferogram. + It can optionally use temporal coherence, correlation, and connected component + information for masking. + + Parameters + ---------- + unw_filenames : list[Path] + List of paths to unwrapped interferogram files to be filtered. + cor_filenames : list[Path] | None + List of paths to correlation files + Passing None skips filtering on correlation. + conncomp_filenames : list[Path] | None + List of paths to connected component files, filters any 0 labeled pixels. + Passing None skips filtering on connected component labels. + temporal_coherence_filename : Path | None + Path to the temporal coherence file for masking. + Passing None skips filtering on temporal coherence. + wavelength_cutoff : float, optional + Spatial wavelength cutoff (in meters) for the filter. Default is 50,000 meters. + correlation_cutoff : float, optional + Threshold of correlation (if passing `cor_filenames`) to use to ignore pixels + during filtering. + output_dir : Path | None, optional + Directory to save the filtered results. + If None, saves in the same location as inputs with .filt.tif extension. + max_workers : int, optional + Number of parallel images to process. Default is 4. + + Returns + ------- + list[Path] + Output filtered rasters. + + Notes + ----- + - If temporal_coherence_filename is provided, pixels with coherence < 0.5 are masked + + """ + from dolphin import io + + bad_pixel_mask = np.zeros( + io.get_raster_xysize(unw_filenames[0])[::-1], dtype="bool" + ) + if temporal_coherence_filename: + bad_pixel_mask = bad_pixel_mask | ( + io.load_gdal(temporal_coherence_filename) < 0.5 + ) + + if output_dir is None: + assert unw_filenames + output_dir = unw_filenames[0].parent + output_dir.mkdir(exist_ok=True) + ctx = mp.get_context("spawn") + + with ProcessPoolExecutor(max_workers, mp_context=ctx) as pool: + return list( + pool.map( + _filter_and_save, + unw_filenames, + cor_filenames or repeat(None), + conncomp_filenames or repeat(None), + repeat(output_dir), + repeat(wavelength_cutoff), + repeat(bad_pixel_mask), + repeat(correlation_cutoff), + ) + ) + + +def _filter_and_save( + unw_filename: Path, + cor_path: Path | None, + conncomp_path: Path | None, + output_dir: Path, + wavelength_cutoff: float, + bad_pixel_mask: NDArray[np.bool_], + correlation_cutoff: float = 0.5, +) -> Path: + """Filter one interferogram (wrapper for multiprocessing).""" + from dolphin import io + from dolphin._overviews import Resampling, create_image_overviews + + # Average for the pixel spacing for filtering + _, x_res, _, _, _, y_res = io.get_raster_gt(unw_filename) + pixel_spacing = (abs(x_res) + abs(y_res)) / 2 + + if cor_path is not None: + bad_pixel_mask |= io.load_gdal(cor_path) < correlation_cutoff + if conncomp_path is not None: + bad_pixel_mask |= io.load_gdal(conncomp_path, masked=True).astype(bool) == 0 + + unw = io.load_gdal(unw_filename) + filt_arr = filter_long_wavelength( + unwrapped_phase=unw, + wavelength_cutoff=wavelength_cutoff, + bad_pixel_mask=bad_pixel_mask, + pixel_spacing=pixel_spacing, + workers=1, + ) + io.round_mantissa(filt_arr, keep_bits=9) + output_name = output_dir / Path(unw_filename).with_suffix(".filt.tif").name + io.write_arr(arr=filt_arr, like_filename=unw_filename, output_name=output_name) + + create_image_overviews(output_name, resampling=Resampling.AVERAGE) + + return output_name + + +def gaussian_filter_nan( + image: ArrayLike, sigma: float | Sequence[float], mode="constant", **kwargs +) -> np.ndarray: + """Apply a gaussian filter to an image with NaNs (avoiding all nans). + + The scipy.ndimage `gaussian_filter` will make the output all NaNs if + any of the pixels in the input that touches the kernel is NaN + + Source: + https://stackoverflow.com/a/36307291 + + Parameters + ---------- + image : ndarray + Image with nans to filter + sigma : float + Size of filter kernel. passed into `gaussian_filter` + mode : str, default = "constant" + Boundary mode for `[scipy.ndimage.gaussian_filter][]` + **kwargs : Any + Passed into `[scipy.ndimage.gaussian_filter][]` + + Returns + ------- + ndarray + Filtered version of `image`. + + """ + from scipy.ndimage import gaussian_filter + + if np.sum(np.isnan(image)) == 0: + return gaussian_filter(image, sigma=sigma, mode=mode, **kwargs) + + V = image.copy() + nan_idxs = np.isnan(image) + V[nan_idxs] = 0 + V_filt = gaussian_filter(V, sigma, **kwargs) + + W = np.ones(image.shape) + W[nan_idxs] = 0 + W_filt = gaussian_filter(W, sigma, **kwargs) + + return V_filt / W_filt