From 41fd1416c0f0917fd878d916ac1388402fd6b503 Mon Sep 17 00:00:00 2001 From: Jungkyo Jung Date: Thu, 30 Apr 2026 01:30:17 +0000 Subject: [PATCH] add gradient mask --- .../isce3/atmosphere/ionosphere_filter.py | 217 ++++++++++++- .../isce3/atmosphere/main_band_estimation.py | 112 +++++++ .../isce3/atmosphere/split_band_estimation.py | 129 ++++++++ python/packages/nisar/workflows/ionosphere.py | 289 ++++++++++++++---- share/nisar/defaults/insar.yaml | 14 + share/nisar/schemas/insar.yaml | 13 +- 6 files changed, 715 insertions(+), 59 deletions(-) diff --git a/python/packages/isce3/atmosphere/ionosphere_filter.py b/python/packages/isce3/atmosphere/ionosphere_filter.py index 62da47386..468539b31 100644 --- a/python/packages/isce3/atmosphere/ionosphere_filter.py +++ b/python/packages/isce3/atmosphere/ionosphere_filter.py @@ -11,7 +11,10 @@ label, find_objects, gaussian_filter, - generic_filter) + generic_filter, + maximum_filter, + uniform_filter, + ) import isce3 from isce3.core.block_param_generator import block_param_generator @@ -1178,3 +1181,215 @@ def remove_small_components(image, min_cluster_pixels): cleaned_image[slice_tuple][current_object] = np.nan return cleaned_image + + + +def build_mask( + phase_image, + *, + coherence_image=None, + conncomp_image=None, + subswath_mask_image=None, + coh_threshold=None, + gradient_threshold=None, + gradient_window_size=25, + gradient_method='mean', + gradient_mask_percentile=None, + water_mask_image=None, + mask_type=None, + fill_value=0.0, +): + """ + Apply a mask to one native-grid phase image. + """ + if phase_image is None: + return None + + if mask_type is None: + mask_type = [] + + mask = np.isfinite(phase_image) & (phase_image != 0) + + if "coherence" in mask_type and coherence_image is not None: + mask &= np.isfinite(coherence_image) + mask &= (coherence_image >= coh_threshold) + + if "connected_components" in mask_type and conncomp_image is not None: + mask &= np.isfinite(conncomp_image) + mask &= (conncomp_image > 0) + + if "subswath_mask" in mask_type and subswath_mask_image is not None: + mask &= np.isfinite(subswath_mask_image) + mask &= (subswath_mask_image > 0) + + if "water" in mask_type and water_mask_image is not None: + mask &= water_mask_image + + if "gradient" in mask_type: + _, _, grad_mask = detect_high_phase_gradient_local( + phase_image, + gradient_threshold, + window_size=gradient_window_size, + mask=mask, + method=gradient_method, + percentile=gradient_mask_percentile) + mask &= grad_mask + + out = phase_image.copy() + out[~mask] = fill_value + return out, mask + + +def detect_high_phase_gradient_local( + unwrapped_phase, + threshold, + window_size=5, + mask=None, + method="mean", + percentile=90, + return_gradient_components=False, +): + """ + Detect high-gradient areas from an unwrapped phase image using + a local windowed gradient score. + Parameters + ---------- + unwrapped_phase : np.ndarray + 2D unwrapped phase array in radians. + threshold : float + Threshold applied to the local gradient score. + Unit: radians/pixel. + window_size : int, optional + Size of the moving window used to summarize gradient magnitude. + Must be >= 1. Typical values: 3, 5, 7. + mask : np.ndarray or None, optional + Boolean valid-data mask with same shape as input. + True means valid pixel. If None, finite values are treated as valid. + method : str, optional + Local summary method: + - "mean": local mean of gradient magnitude + - "max": local max of gradient magnitude + - "percentile": local percentile of gradient magnitude + percentile : float, optional + Percentile used when method="percentile". + return_gradient_components : bool, optional + If True, also return gx and gy. + Returns + ------- + grad_mag : np.ndarray + Pixelwise gradient magnitude in radians/pixel. + local_grad_score : np.ndarray + Local windowed summary of grad_mag in radians/pixel. + high_gradient_mask : np.ndarray + Boolean mask where local_grad_score > threshold. + gx, gy : np.ndarray, optional + Gradient components if return_gradient_components=True. + Notes + ----- + For unwrapped phase, gradients are computed using central differences + for interior pixels and first differences at the boundaries. + """ + phase = np.asarray(unwrapped_phase, dtype=np.float64) + + if phase.ndim != 2: + raise ValueError("unwrapped_phase must be a 2D array") + + if window_size < 1 or int(window_size) != window_size: + raise ValueError("window_size must be a positive integer") + + if method not in {"mean", "max", "percentile"}: + raise ValueError("method must be 'mean', 'max', or 'percentile'") + + if mask is None: + valid = np.isfinite(phase) + else: + valid = np.asarray(mask, dtype=bool) & np.isfinite(phase) + + # Require valid neighbors for safe gradient estimation + # We will compute gradients first on a filled array, then invalidate + # locations influenced by invalid pixels. + phase_filled = phase.copy() + phase_filled[~valid] = 0.0 + + # Gradient in y (rows) and x (cols) + gy = np.full_like(phase, np.nan, dtype=np.float64) + gx = np.full_like(phase, np.nan, dtype=np.float64) + + # Central differences for interior + gy[1:-1, :] = (phase_filled[2:, :] - phase_filled[:-2, :]) / 2.0 + gx[:, 1:-1] = (phase_filled[:, 2:] - phase_filled[:, :-2]) / 2.0 + + # Forward/backward difference at borders + gy[0, :] = phase_filled[1, :] - phase_filled[0, :] + gy[-1, :] = phase_filled[-1, :] - phase_filled[-2, :] + + gx[:, 0] = phase_filled[:, 1] - phase_filled[:, 0] + gx[:, -1] = phase_filled[:, -1] - phase_filled[:, -2] + + # Build validity masks for gradients + valid_gx = np.zeros_like(valid, dtype=bool) + valid_gy = np.zeros_like(valid, dtype=bool) + + # gx validity + valid_gx[:, 1:-1] = valid[:, 2:] & valid[:, :-2] & valid[:, 1:-1] + valid_gx[:, 0] = valid[:, 0] & valid[:, 1] + valid_gx[:, -1] = valid[:, -1] & valid[:, -2] + + # gy validity + valid_gy[1:-1, :] = valid[2:, :] & valid[:-2, :] & valid[1:-1, :] + valid_gy[0, :] = valid[0, :] & valid[1, :] + valid_gy[-1, :] = valid[-1, :] & valid[-2, :] + + gx[~valid_gx] = np.nan + gy[~valid_gy] = np.nan + + # Gradient magnitude + grad_mag = np.sqrt(gx**2 + gy**2) + + # Valid where both phase and grad_mag are valid + valid_grad = np.isfinite(grad_mag) & valid + + # Local summary of gradient magnitude + if method == "mean": + data = np.where(valid_grad, grad_mag, 0.0) + weights = valid_grad.astype(np.float64) + + local_sum = uniform_filter(data, size=window_size, + mode="constant", cval=0.0) * (window_size ** 2) + local_count = uniform_filter(weights, size=window_size, mode="constant", cval=0.0) * (window_size ** 2) + + with np.errstate(invalid="ignore", divide="ignore"): + local_grad_score = local_sum / local_count + + local_grad_score[local_count == 0] = np.nan + + elif method == "max": + data = np.where(valid_grad, grad_mag, -np.inf) + local_grad_score = maximum_filter(data, size=window_size, mode="constant", cval=-np.inf) + local_grad_score[~valid] = np.nan + local_grad_score[~np.isfinite(local_grad_score)] = np.nan + + elif method == "percentile": + def nanpercentile_func(values): + vals = values[np.isfinite(values)] + if vals.size == 0: + return np.nan + return np.percentile(vals, percentile) + + data = np.where(valid_grad, grad_mag, np.nan) + local_grad_score = generic_filter( + data, + function=nanpercentile_func, + size=window_size, + mode="constant", + cval=np.nan, + ) + + # Final mask + high_gradient_mask = np.isfinite(local_grad_score) & (local_grad_score < threshold) + high_gradient_mask[~valid] = False + + if return_gradient_components: + return grad_mag, local_grad_score, high_gradient_mask, gx, gy + + return grad_mag, local_grad_score, high_gradient_mask diff --git a/python/packages/isce3/atmosphere/main_band_estimation.py b/python/packages/isce3/atmosphere/main_band_estimation.py index b8c58967d..365fdb582 100644 --- a/python/packages/isce3/atmosphere/main_band_estimation.py +++ b/python/packages/isce3/atmosphere/main_band_estimation.py @@ -3,6 +3,8 @@ import numpy as np from .ionosphere_estimation import IonosphereEstimation +from isce3.atmosphere.ionosphere_filter import ( + detect_high_phase_gradient_local) from isce3.signal.interpolate_by_range import decimate_freq_a_array from isce3.unwrap.preprocess import interpret_subswath_mask @@ -367,6 +369,116 @@ def get_valid_area( return mask_array + def get_gradient_mask( + self, + main_array=None, + side_array=None, + diff_ms_array=None, + low_band_array=None, + high_band_array=None, + diff_low_high_band_array=None, + slant_main=None, + slant_side=None, + mask=None, + window=None, + threshold_first=None, + threshold_second=None, + method='mean', + percentile=None): + """ + Generate a boolean mask identifying pixels with consistently high phase + gradient based on one or two input interferometric arrays. + This function detects high phase gradient regions using + `detect_high_phase_gradient_local` and combines two gradient masks + (first and second) using a logical AND operation. The resulting mask + highlights pixels that simultaneously exceed gradient thresholds in + both inputs. Optionally, frequency A/B alignment (decimation) is applied + when side-band data is involved. + The function supports two modes: + 1. Using `main_array` and `diff_low_high_band_array` + 2. Using `low_band_array` and `high_band_array` + Small isolated pixels are removed from the final mask. + Parameters + ---------- + main_array : numpy.ndarray, optional + Interferometric phase (or coherence) array for the main band. + Used when `diff_low_high_band_array` is provided. + side_array : numpy.ndarray, optional + Interferometric array for the side band. If provided, low/high band + arrays are decimated to match the geometry of the main band. + diff_ms_array : numpy.ndarray, optional + Difference array between main and side bands (currently unused). + low_band_array : numpy.ndarray, optional + Interferometric array for the low-frequency sub-band. + high_band_array : numpy.ndarray, optional + Interferometric array for the high-frequency sub-band. + diff_low_high_band_array : numpy.ndarray, optional + Difference array between high and low sub-band interferograms. + If provided, this is used instead of separate low/high arrays. + slant_main : numpy.ndarray, optional + Slant range coordinates for the main (frequency A) band. + slant_side : numpy.ndarray, optional + Slant range coordinates for the side (frequency B) band. + window : int, optional + Window size used for local gradient estimation. + threshold_first : float + Threshold for the first gradient mask (e.g., main or low band). + threshold_second : float + Threshold for the second gradient mask (e.g., diff or high band). + method : str, optional + Method used for gradient aggregation within the local window. + Options depend on `detect_high_phase_gradient_local` + (e.g., 'mean', 'max', 'percentile). + percentile : float, optional + Percentile value used when method requires percentile-based + thresholding. + Returns + ------- + mask_array : numpy.ndarray (bool) + 2D boolean mask where: + - True (1): pixels with high phase gradient in both inputs + - False (0): pixels not meeting the criteria + The mask is post-processed to remove isolated single pixels. + """ + # decimate coherence or connected components + # when side array is also used. + if side_array is not None or diff_ms_array is not None: + main_array = decimate_freq_a_array( + slant_main, + slant_side, + main_array) + + _, _, grad_mask_first = detect_high_phase_gradient_local( + main_array, + threshold_first, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + + if diff_ms_array is not None: + _, _, grad_mask_second = detect_high_phase_gradient_local( + diff_ms_array, + threshold_second, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + elif side_array is not None: + _, _, grad_mask_second = detect_high_phase_gradient_local( + side_array, + threshold_second, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + + mask_array = grad_mask_first & grad_mask_second + mask_array = self.remove_single_pixels(mask_array) + + mask_array = np.array(mask_array, dtype=bool) + return mask_array + def get_subswath_mask_array( self, main_array=None, diff --git a/python/packages/isce3/atmosphere/split_band_estimation.py b/python/packages/isce3/atmosphere/split_band_estimation.py index 2f4cdee6b..0d5e0b20e 100644 --- a/python/packages/isce3/atmosphere/split_band_estimation.py +++ b/python/packages/isce3/atmosphere/split_band_estimation.py @@ -3,6 +3,8 @@ import numpy as np from .ionosphere_estimation import IonosphereEstimation +from isce3.atmosphere.ionosphere_filter import ( + detect_high_phase_gradient_local) from isce3.signal.interpolate_by_range import decimate_freq_a_array from isce3.unwrap.preprocess import interpret_subswath_mask @@ -457,6 +459,133 @@ def get_valid_area( return mask_array + def get_gradient_mask( + self, + main_array=None, + side_array=None, + diff_ms_array=None, + low_band_array=None, + high_band_array=None, + diff_low_high_band_array=None, + mask=None, + slant_main=None, + slant_side=None, + window=None, + threshold_first=None, + threshold_second=None, + method='mean', + percentile=None): + """ + Generate a boolean mask identifying pixels with consistently high phase + gradient based on one or two input interferometric arrays. + This function detects high phase gradient regions using + `detect_high_phase_gradient_local` and combines two gradient masks + (first and second) using a logical AND operation. The resulting mask + highlights pixels that simultaneously exceed gradient thresholds in + both inputs. Optionally, frequency A/B alignment (decimation) is applied + when side-band data is involved. + The function supports two modes: + 1. Using `main_array` and `diff_low_high_band_array` + 2. Using `low_band_array` and `high_band_array` + Small isolated pixels are removed from the final mask. + Parameters + ---------- + main_array : numpy.ndarray, optional + Interferometric phase (or coherence) array for the main band. + Used when `diff_low_high_band_array` is provided. + side_array : numpy.ndarray, optional + Interferometric array for the side band. If provided, low/high band + arrays are decimated to match the geometry of the main band. + diff_ms_array : numpy.ndarray, optional + Difference array between main and side bands (currently unused). + low_band_array : numpy.ndarray, optional + Interferometric array for the low-frequency sub-band. + high_band_array : numpy.ndarray, optional + Interferometric array for the high-frequency sub-band. + diff_low_high_band_array : numpy.ndarray, optional + Difference array between high and low sub-band interferograms. + If provided, this is used instead of separate low/high arrays. + slant_main : numpy.ndarray, optional + Slant range coordinates for the main (frequency A) band. + slant_side : numpy.ndarray, optional + Slant range coordinates for the side (frequency B) band. + window : int, optional + Window size used for local gradient estimation. + threshold_first : float + Threshold for the first gradient mask (e.g., main or low band). + threshold_second : float + Threshold for the second gradient mask (e.g., diff or high band). + method : str, optional + Method used for gradient aggregation within the local window. + Options depend on `detect_high_phase_gradient_local` + (e.g., 'mean', 'max', 'percentile). + percentile : float, optional + Percentile value used when method requires percentile-based + thresholding. + Returns + ------- + mask_array : numpy.ndarray (bool) + 2D boolean mask where: + - True (1): pixels with high phase gradient in both inputs + - False (0): pixels not meeting the criteria + The mask is post-processed to remove isolated single pixels. + """ + # decimate coherence or connected components + # when side array is also used. + if side_array is not None: + if slant_main is None: + slant_main = self.slant_main + if slant_side is None: + slant_side = self.slant_side + + if low_band_array is not None: + low_band_array = decimate_freq_a_array( + slant_main, + slant_side, + low_band_array) + if high_band_array is not None: + high_band_array = decimate_freq_a_array( + slant_main, + slant_side, + high_band_array) + + if diff_low_high_band_array is not None: + _, _, grad_mask_first = detect_high_phase_gradient_local( + main_array, + threshold_first, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + + _, _, grad_mask_second = detect_high_phase_gradient_local( + diff_low_high_band_array, + threshold_second, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + + elif low_band_array is not None and high_band_array is not None: + _, _, grad_mask_first = detect_high_phase_gradient_local( + low_band_array, + threshold_first, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + + _, _, grad_mask_second = detect_high_phase_gradient_local( + high_band_array, + threshold_second, + window_size=window, + mask=mask, + method=method, + percentile=percentile) + mask_array = grad_mask_first & grad_mask_second + mask_array = self.remove_single_pixels(mask_array) + return mask_array + def estimate_iono_std( self, main_coh=None, diff --git a/python/packages/nisar/workflows/ionosphere.py b/python/packages/nisar/workflows/ionosphere.py index 5be4a9647..b82a72ac5 100644 --- a/python/packages/nisar/workflows/ionosphere.py +++ b/python/packages/nisar/workflows/ionosphere.py @@ -10,7 +10,10 @@ import journal import numpy as np from isce3.atmosphere.ionosphere_filter import ( - IonosphereFilter, read_block_array, unwrapping_correction_with_filter, + IonosphereFilter, + build_mask, + read_block_array, + unwrapping_correction_with_filter, write_array) from isce3.atmosphere.main_band_estimation import ( MainDiffMsBandIonosphereEstimation, MainSideBandIonosphereEstimation) @@ -366,6 +369,85 @@ def copy_iono_datasets(iono_insar_cfg, compute_stats_real_hdf5_dataset(dst_h5[dst_iono_path]) +def apply_bridge_with_mask( + phase_image, + *, + shared_params, + coherence_image=None, + conncomp_image=None, + subswath_mask_image=None, + water_mask_image=None, + gradient_threshold=None, + zero_to_nan=False, +): + """ + Apply pre-bridge masking and bridge algorithm using shared parameters. + + Parameters + ---------- + phase_image : numpy.ndarray + Input phase image. + shared_params : dict + Dictionary containing common mask and bridge parameters. + coherence_image : numpy.ndarray, optional + Coherence image used for masking. + conncomp_image : numpy.ndarray, optional + Connected component image used for masking. + subswath_mask_image : numpy.ndarray, optional + Subswath mask image used for masking. + water_mask_image : numpy.ndarray, optional + Water mask image used for masking. + gradient_threshold : float, optional + Per-image gradient threshold. If None, use the value from shared_params. + zero_to_nan : bool, optional + Convert zeros in bridged result to NaN. + + Returns + ------- + bridged_image : numpy.ndarray + Bridged output image. + mask : numpy.ndarray + Boolean mask used before bridging. + """ + if phase_image is None: + return None, None + + local_gradient_threshold = ( + mask_params['gradient_threshold'] + if gradient_threshold is None + else gradient_threshold + ) + mask_params = shared_params['mask'] + masked_image, mask = build_mask( + phase_image, + coherence_image=coherence_image, + conncomp_image=conncomp_image, + subswath_mask_image=subswath_mask_image, + water_mask_image=water_mask_image, + coh_threshold=mask_params['coh_threshold'], + gradient_threshold=local_gradient_threshold, + gradient_window_size=mask_params['gradient_window_size'], + gradient_method=mask_params['gradient_method'], + gradient_mask_percentile=mask_params['gradient_mask_percentile'], + mask_type=mask_params['mask_type'], + fill_value=mask_params.get('fill_value', 0.0), + ) + bridge_params = shared_params['bridge'] + bridged_image = bridge_unwrapped_phase( + masked_image, + radius=bridge_params['radius'], + min_num_pixel=bridge_params['min_num_pixel'], + erosion_size=bridge_params['erosion_size'], + ramp_type=bridge_params['ramp_type'], + deramp_max_num_sample=bridge_params['deramp_max_num_sample'], + ) + + if zero_to_nan: + bridged_image[bridged_image == 0] = np.nan + + return bridged_image, mask + + def compute_differential_phase( phase_first, phase_second, @@ -1072,7 +1154,6 @@ def run(cfg: dict, runw_hdf5: str): # pull parameters for dispersive filter filter_bool = filter_cfg['enabled'] mask_type = filter_cfg['filter_mask_type'] - filter_coh_thresh = filter_cfg['filter_coherence_threshold'] kernel_range_size = filter_cfg['kernel_range'] kernel_azimuth_size = filter_cfg['kernel_azimuth'] kernel_sigma_range = filter_cfg['sigma_range'] @@ -1090,6 +1171,15 @@ def run(cfg: dict, runw_hdf5: str): min_cluster_pixels = filter_cfg['min_cluster_pixels'] unwrap_correction_bool = filter_cfg['unwrap_correction'] + # Mask parameter + mask_type = filter_cfg['filter_mask_type'] + filter_coh_thresh = filter_cfg['filter_coherence_threshold'] + gradient_mask_method = filter_cfg['gradient_mask_method'] + gradient_mask_window = filter_cfg['gradient_mask_window'] + gradient_mask_threshold_first = filter_cfg['gradient_mask_threshold_first'] + gradient_mask_threshold_second = filter_cfg['gradient_mask_threshold_second'] + gradient_mask_percentile = filter_cfg['gradient_mask_percentile'] + # bridge algorithm options bridge_algorithm_bool = filter_cfg['bridge_algorithm_enabled'] bridge_minimum_samples = filter_cfg['bridge_minimum_samples'] @@ -1098,6 +1188,24 @@ def run(cfg: dict, runw_hdf5: str): bridge_deramp_type = filter_cfg['bridge_ramp_type'] bridge_ramp_maximum_pixel = filter_cfg['bridge_ramp_maximum_pixel'] + shared_bridge_mask_params = { + 'mask': { + 'coh_threshold': filter_coh_thresh, + 'gradient_window_size': gradient_mask_window, + 'gradient_method': gradient_mask_method, + 'gradient_mask_percentile': gradient_mask_percentile, + 'mask_type': mask_type, + 'fill_value': 0.0, + }, + 'bridge': { + 'radius': bridge_radius, + 'min_num_pixel': bridge_minimum_samples, + 'erosion_size': bridge_erosion_size, + 'ramp_type': bridge_deramp_type, + 'deramp_max_num_sample': bridge_ramp_maximum_pixel, + } + } + rg_looks = cfg['processing']['crossmul']['range_looks'] az_looks = cfg['processing']['crossmul']['azimuth_looks'] unwrap_rg_looks = cfg['processing']['phase_unwrap']['range_looks'] @@ -1121,7 +1229,7 @@ def run(cfg: dict, runw_hdf5: str): # Run InSAR for sub-band SLCs (split-main-bands) or # for main and side bands for iono_freq_pols (main-side-bands) - insar_ionosphere_pair(iono_insar_cfg, runw_hdf5) + # insar_ionosphere_pair(iono_insar_cfg, runw_hdf5) t_all = time.time() # Define methods to use subband or sideband @@ -1216,6 +1324,19 @@ def run(cfg: dict, runw_hdf5: str): mad_scale_factor=filling_outlier_mad_scale_factor, outputdir=os.path.join(iono_path, iono_method)) + water_mask_b_blk = None + water_mask_a_blk = None + if "water" in mask_type and filter_bool: + water_mask_path = cfg[ + "dynamic_ancillary_file_group"]["water_mask_file"] + + water_distance_a = project_map_to_radar(cfg, water_mask_path, 'A') + water_mask_a = (water_distance_a == 0) + + if iono_method in iono_method_sideband: + water_distance_b = project_map_to_radar(cfg, water_mask_path, 'B') + water_mask_b = (water_distance_b == 0) + # pull parameters for polarizations pol_list_a = list(iono_freq_pols['A']) if iono_method in iono_method_sideband: @@ -1389,6 +1510,10 @@ def run(cfg: dict, runw_hdf5: str): [block_rows_data, cols_main], dtype=int) + if "water" in mask_type: + water_mask_a_blk = None if water_mask_a is None else \ + water_mask_a[row_start:row_start + block_rows_data, :] + if iono_method == 'main_diff_low_high_subband': main_image = np.empty([block_rows_data, cols_main], dtype=float) @@ -1443,40 +1568,6 @@ def run(cfg: dict, runw_hdf5: str): subswath_mask_image, np.s_[row_start:row_start + block_rows_data, :]) - if bridge_algorithm_bool: - sub_high_image = bridge_unwrapped_phase( - sub_high_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) - sub_low_image = bridge_unwrapped_phase( - sub_low_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) - - if unwrap_correction_bool: - sub_high_image = unwrapping_correction_with_filter( - sub_high_image, - kernel_width=kernel_range_size, - kernel_length=kernel_azimuth_size, - sig_kernel_x=kernel_sigma_range, - sig_kernel_y=kernel_sigma_azimuth, - iterations=filter_iterations, - filter_method='convolution') - sub_low_image = unwrapping_correction_with_filter( - sub_low_image, - kernel_width=kernel_range_size, - kernel_length=kernel_azimuth_size, - sig_kernel_x=kernel_sigma_range, - sig_kernel_y=kernel_sigma_azimuth, - iterations=filter_iterations, - filter_method='convolution') - if iono_method == 'main_diff_low_high_subband': with HDF5OptimizedReader( name=runw_path_insar, mode='r', @@ -1518,6 +1609,54 @@ def run(cfg: dict, runw_hdf5: str): row_start:row_start + block_rows_data, :]) + if "water" in mask_type: + water_mask_a_blk = \ + water_mask_a[row_start:row_start + block_rows_data, :] + + if bridge_algorithm_bool: + if iono_method == "split_main_band": + sub_low_image, _ = apply_bridge_with_mask( + sub_low_image, + shared_params=shared_bridge_mask_params, + coherence_image=sub_low_coh_image, + conncomp_image=sub_low_conn_image, + subswath_mask_image=subswath_mask_image, + water_mask_image=water_mask_a_blk, + gradient_threshold=gradient_mask_threshold_first, + ) + + sub_high_image, _ = apply_bridge_with_mask( + sub_high_image, + shared_params=shared_bridge_mask_params, + coherence_image=sub_high_coh_image, + conncomp_image=sub_high_conn_image, + subswath_mask_image=subswath_mask_image, + water_mask_image=water_mask_a_blk, + gradient_threshold=gradient_mask_threshold_second, + ) + elif iono_method == "main_diff_low_high_subband": + main_image, _ = apply_bridge_with_mask( + main_image, + shared_params=shared_bridge_mask_params, + coherence_image=main_coh_image, + conncomp_image=main_conn_image, + subswath_mask_image=subswath_mask_image, + water_mask_image=water_mask_a_blk, + gradient_threshold=gradient_mask_threshold_first, + zero_to_nan=True, + ) + + diff_subband_image, _ = apply_bridge_with_mask( + diff_subband_image, + shared_params=shared_bridge_mask_params, + coherence_image=diff_coh_image, + conncomp_image=diff_subband_conn_image, + subswath_mask_image=subswath_mask_image, + water_mask_image=water_mask_a_blk, + gradient_threshold=gradient_mask_threshold_second, + zero_to_nan=True, + ) + if iono_method in iono_method_sideband: main_image = np.empty([block_rows_data, cols_main], @@ -1555,6 +1694,15 @@ def run(cfg: dict, runw_hdf5: str): [block_rows_data, cols_side], dtype=int) + if "water" in mask_type: + water_mask_a_blk = None if water_mask_a is None \ + else water_mask_a[row_start:row_start + + block_rows_data, :] + + water_mask_b_blk = None if water_mask_b is None \ + else water_mask_b[row_start:row_start + + block_rows_data, :] + with HDF5OptimizedReader( name=runw_freq_a_str, mode='r', libver='latest', swmr=True) as src_main_h5, \ @@ -1626,32 +1774,43 @@ def run(cfg: dict, runw_hdf5: str): diff_ms_conn_image, np.s_[row_start:row_start + block_rows_data, :]) + if "water" in mask_type: + water_mask_a_blk = \ + water_mask_a[row_start:row_start + block_rows_data, :] + water_mask_b_blk = \ + water_mask_b[row_start:row_start + block_rows_data, :] + if bridge_algorithm_bool: - main_image = bridge_unwrapped_phase( + main_image, _ = apply_bridge_with_mask( main_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) + shared_params=shared_bridge_mask_params, + coherence_image=main_coh_image, + conncomp_image=main_conn_image, + subswath_mask_image=subswath_mask_main_image, + water_mask_image=water_mask_a_blk, + gradient_threshold=gradient_mask_threshold_first, + ) if iono_method == 'main_side_band': - side_image = bridge_unwrapped_phase( + side_image, _ = apply_bridge_with_mask( side_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) + shared_params=shared_bridge_mask_params, + coherence_image=side_coh_image, + conncomp_image=side_conn_image, + subswath_mask_image=subswath_mask_side_image, + water_mask_image=water_mask_b_blk, + gradient_threshold=gradient_mask_threshold_second, + ) elif iono_method == 'main_diff_ms_band': - diff_ms_image = bridge_unwrapped_phase( + diff_ms_image, _ = apply_bridge_with_mask( diff_ms_image, - radius=bridge_radius, - min_num_pixel=bridge_minimum_samples, - erosion_size=bridge_erosion_size, - ramp_type=bridge_deramp_type, - deramp_max_num_sample=bridge_ramp_maximum_pixel) - + shared_params=shared_bridge_mask_params, + coherence_image=diff_ms_coh_image, + conncomp_image=diff_ms_conn_image, + subswath_mask_image=subswath_mask_side_image, + water_mask_image=water_mask_b_blk, + gradient_threshold=gradient_mask_threshold_second, + ) if unwrap_correction_bool: main_image = unwrapping_correction_with_filter( main_image, @@ -1830,6 +1989,24 @@ def run(cfg: dict, runw_hdf5: str): slant_side=side_slant) mask_array &= mask_subswath + if "gradient" in mask_type: + mask_gradient = iono_phase_obj.get_gradient_mask( + main_array=main_image, + side_array=side_image, + diff_ms_array=diff_ms_image, + low_band_array=sub_low_image, + high_band_array=sub_high_image, + diff_low_high_band_array=diff_subband_image, + slant_main=main_slant, + slant_side=side_slant, + mask=mask_array, + window=gradient_mask_window, + threshold_first=gradient_mask_threshold_first, + threshold_second=gradient_mask_threshold_second, + method=gradient_mask_method, + percentile=gradient_mask_percentile) + mask_array = np.logical_and(mask_array, mask_gradient) + if "water" in mask_type: # Extract preprocessing dictionary and open arrays # water_mask_file is expected to have distance from the diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml index b85fb632a..86661050c 100644 --- a/share/nisar/defaults/insar.yaml +++ b/share/nisar/defaults/insar.yaml @@ -287,6 +287,20 @@ runconfig: median_filter_threshold: 0.6 # Coherence thresholds 0<= thresholds <=1 filter_coherence_threshold: 0.5 + # Gradient estimation window + gradient_mask_window: 25 + # Gradient estimation method + gradient_mask_method: 'mean' + # Gradient estimation threshold + # split_main_band: first -> low subband, second -> high subband + # main_diff_low_high_subband: first-> main band, second -> diff + # main_side_band: first -> main band, second -> side band + # main_diff_ms_band: first -> main_band, second-> diff + gradient_mask_threshold_first: 1 + gradient_mask_threshold_second: 1 + # Gredient method percentile + gradient_mask_percentile: 90 + # Filling method for invalid regions # (nearest, smoothed) filling_method: nearest diff --git a/share/nisar/schemas/insar.yaml b/share/nisar/schemas/insar.yaml index c3999eb75..6a92eb9a5 100644 --- a/share/nisar/schemas/insar.yaml +++ b/share/nisar/schemas/insar.yaml @@ -789,13 +789,22 @@ dispersive_filter_options: # When multiple validation sources are enabled, take the intersection of # their invalid masks to form a combined mask. # Masked areas are then filled using spatial interpolation. - filter_mask_type: list('water', 'median_filter', 'coherence', 'connected_components', 'subswath_mask', required=False) - # Median filter size for mask generation + filter_mask_type: list(enum('water', 'median_filter', 'coherence', 'connected_components', 'subswath_mask', 'gradient'), required=False) # Median filter size for mask generation median_filter_size: num(min=1, required=False) # Median filter threshold for mask generation median_filter_threshold: num(min=0, max=1, required=False) # coherence threshold for mask generation filter_coherence_threshold: num(min=0, max=1, required=False) + # Gradient estimation window + gradient_mask_window: num(min=1, required=False) + # Gradient estimation method + gradient_mask_method: enum('mean', 'percentile', 'max', required=False) + # Gradient estimation threshold + gradient_mask_threshold_first: num(min=0, required=False) + gradient_mask_threshold_second: num(min=0, required=False) + # Gredient method percentile + gradient_mask_percentile: num(min=0, required=False) + # Filling method for invalid regions # nearest : apply nearest interpolation for whole invalid regions # smoothed : apply lienar interpolation and nearest extrapolations