From 00b0bbcee9f2521211bd53c889f70369af99e8e1 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Thu, 2 Jan 2025 14:04:05 -0500 Subject: [PATCH] clean up commented, unused --- src/dolphin/filtering.py | 58 ++++++++++------------------------------ 1 file changed, 14 insertions(+), 44 deletions(-) diff --git a/src/dolphin/filtering.py b/src/dolphin/filtering.py index 0fb66292..a1b33f42 100644 --- a/src/dolphin/filtering.py +++ b/src/dolphin/filtering.py @@ -9,7 +9,6 @@ import numpy as np from numpy.typing import ArrayLike, NDArray -from osgeo import gdal from osgeo_utils.gdal_fillnodata import gdal_fillnodata from scipy import fft, ndimage @@ -76,14 +75,10 @@ def filter_long_wavelength( # Work on a copy of the displacement field to avoid modifying the input displacement = np.nan_to_num(unwrapped_phase) - good_pixel_mask = np.logical_not(bad_pixel_mask) - # Take either nan or 0 pixels in `unwrapped_phase` to be nodata nodata_mask = displacement == 0 in_bounds_pixels = np.logical_not(nodata_mask) - total_valid_mask = in_bounds_pixels & good_pixel_mask - # Convert bad pixels to NaN for GDAL fillnodata displacement[bad_pixel_mask] = np.nan @@ -104,12 +99,13 @@ def filter_long_wavelength( temp_dst = tmp_path / "filled.tif" # Save the array to a temporary GeoTIFF - driver = gdal.GetDriverByName("GTiff") - dset = driver.Create(str(temp_src), cols, rows, 1, gdal.GDT_Float32) - band = dset.GetRasterBand(1) - band.WriteArray(displacement) - band.SetNoDataValue(np.nan) - dset = None # Close the dset + io.write_arr( + arr=displacement, + nodata=np.nan, + output_name=temp_src, + shape=(rows, cols), + dtype=displacement.dtype, + ) # Fill nodata using GDAL gdal_fillnodata( @@ -123,44 +119,18 @@ def filter_long_wavelength( filled_data = io.load_gdal(temp_dst, masked=True) - padded = filled_data - - # We need to find the new pixels which are valid after filling - # These are going to be out of bounds in the original image, - # but we don't want to ignore them during the low pass filter. - # The boundary pixels are the ones which will allow the edge values - # to be smoothed out without. - new_valid = np.isnan(displacement) & ~np.isnan(filled_data) - - bad_pixel_mask_ext = bad_pixel_mask.copy() - bad_pixel_mask_ext[new_valid] = False - - # # Remove the plane, setting to 0 where we had no data for the plane fit: - # unw_ifg_interp = np.where(total_valid_mask, displacement, 0) - - # # 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( - # unw_ifg_interp, ((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) + lowpass_filtered = fft.fft2(filled_data, workers=workers) + lowpass_filtered = ndimage.fourier_gaussian(lowpass_filtered, 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] - lowpass_filtered = result + lowpass_filtered = fft.ifft2(lowpass_filtered, workers=workers).real.astype( + unwrapped_phase.dtype + ) filtered_ifg = filled_data - lowpass_filtered * in_bounds_pixels if fill_value is not None: + good_pixel_mask = np.logical_not(bad_pixel_mask) + total_valid_mask = in_bounds_pixels & good_pixel_mask return np.where(total_valid_mask, filtered_ifg, fill_value) else: return filtered_ifg