Merge pull request #69 from jhkennedy/s1-filter-widtha

Use a different filter width for Sentinel-1
nasa-jpl · Nov 1, 2022 · ade077f · ade077f
2 parents 3a2a17a + 038326f
commit ade077f
Show file tree

Hide file tree

Showing 3 changed files with 45 additions and 18 deletions.
diff --git a/geo_autoRIFT/autoRIFT/autoRIFT.py b/geo_autoRIFT/autoRIFT/autoRIFT.py
@@ -64,17 +64,17 @@ def _wallis_filter(image, filter_width):
 
 
 def _wallis_filter_fill(image, filter_width, std_cutoff):
-    image = np.ma.masked_values(image, 0.0)
+    invalid_data = np.isclose(image, 0.0)
     buff = np.sqrt(2 * ((filter_width - 1) / 2) ** 2) + 0.01
 
     # find edges of image, this makes missing scan lines valid and will
     # later be filled with random white noise
-    potential_data = distance_transform_edt(image.mask) < 30
-    missing_data = potential_data & image.mask
-    missing_data = distance_transform_edt(~missing_data) <= buff
+    potential_data = cv2.distanceTransform(invalid_data.astype(np.uint8), distanceType=cv2.DIST_L2, maskSize=cv2.DIST_MASK_PRECISE) < 30
+    missing_data = potential_data & invalid_data
+    missing_data = cv2.distanceTransform((~missing_data).astype(np.uint8), distanceType=cv2.DIST_L2, maskSize=cv2.DIST_MASK_PRECISE) <= buff
 
     # trying to frame out the image
-    valid_domain = ~image.mask | missing_data
+    valid_domain = ~invalid_data | missing_data
     zero_mask = ~valid_domain
 
     kernel = np.ones((filter_width, filter_width), dtype=np.float32)
@@ -84,13 +84,13 @@ def _wallis_filter_fill(image, filter_width, std_cutoff):
     std = _preprocess_filt_std(image, kernel)
 
     low_std = std < std_cutoff
-    low_std = distance_transform_edt(~low_std) <= buff
+    low_std = cv2.distanceTransform((~low_std).astype(np.uint8), distanceType=cv2.DIST_L2, maskSize=cv2.DIST_MASK_PRECISE) <= buff
     missing_data = (missing_data | low_std) & valid_domain
 
     image = shifted / std
 
     valid_data = valid_domain & ~missing_data
-    image.mask |= ~valid_data
+    invalid_data |= ~valid_data
 
     # wallis filter normalizes the imagery to have a mean=0 and std=1;
     # fill with random values from a normal distribution with same mean and std

diff --git a/testautoRIFT.py b/testautoRIFT.py
@@ -130,7 +130,8 @@ def loadProductOptical(file_m, file_s):
 
 
 def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, optflag,
-                nodata, mpflag, geogrid_run_info=None, preprocessing_methods=('hps', 'hps')):
+                nodata, mpflag, geogrid_run_info=None, preprocessing_methods=('hps', 'hps'),
+                preprocessing_filter_width=5):
     '''
     Wire and run geogrid.
     '''
@@ -145,7 +146,10 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS
 
 
     obj = autoRIFT()
-#    obj.configure()
+    # obj.configure()
+
+    obj.WallisFilterWidth = preprocessing_filter_width
+    print(f'Setting Wallis Filter Width to {preprocessing_filter_width}')
 
 #    ##########     uncomment if starting from preprocessed images
 #    I1 = I1.astype(np.uint8)
@@ -281,6 +285,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS
     ######## preprocessing
     t1 = time.time()
     print("Pre-process Start!!!")
+    print(f"Using Wallis Filter Width: {obj.WallisFilterWidth}")
 #    obj.zeroMask = 1
 
     # TODO: Allow different filters to be applied images independently
@@ -518,6 +523,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search
         m_name = os.path.basename(indir_m)
         s_name = os.path.basename(indir_s)
 
+        # FIXME: Filter width is a magic variable here and not exposed well.
+        preprocessing_filter_width = 5
+        if nc_sensor == 'S1':
+            preprocessing_filter_width = 21
+
+        print(f'Preprocessing filter width {preprocessing_filter_width}')
+
         preprocessing_methods = ['hps', 'hps']
         for ii, name in enumerate((m_name, s_name)):
             if len(re.findall("L[EO]07_", name)) > 0:
@@ -529,9 +541,11 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search
 
         print(f'Using preprocessing methods {preprocessing_methods}')
 
-        Dx, Dy, InterpMask, ChipSizeX, GridSpacingX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask = runAutorift(
-            data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0,
-            noDataMask, optical_flag, nodata, mpflag, geogrid_run_info=geogrid_run_info, preprocessing_methods=preprocessing_methods,
+        Dx, Dy, InterpMask, ChipSizeX, GridSpacingX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask = \
+            runAutorift(
+                data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0,
+                noDataMask, optical_flag, nodata, mpflag, geogrid_run_info=geogrid_run_info,
+                preprocessing_methods=preprocessing_methods, preprocessing_filter_width=preprocessing_filter_width,
         )
         if nc_sensor is not None:
             import netcdf_output as no
@@ -618,7 +632,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search
 
 
         if offset2vx is not None:
-        
+
             ds = gdal.Open(scale_factor)
             band = ds.GetRasterBand(1)
             scale_factor_1 = band.ReadAsArray()

diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py
@@ -130,7 +130,8 @@ def loadProductOptical(file_m, file_s):
 
 
 def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, optflag,
-                nodata, mpflag, geogrid_run_info=None, preprocessing_methods=('hps', 'hps')):
+                nodata, mpflag, geogrid_run_info=None, preprocessing_methods=('hps', 'hps'),
+                preprocessing_filter_width=5):
     '''
     Wire and run geogrid.
     '''
@@ -146,6 +147,8 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS
     obj = autoRIFT_ISCE()
     obj.configure()
 
+    obj.WallisFilterWidth = preprocessing_filter_width
+
 #    ##########     uncomment if starting from preprocessed images
 #    I1 = I1.astype(np.uint8)
 #    I2 = I2.astype(np.uint8)
@@ -280,6 +283,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS
     ######## preprocessing
     t1 = time.time()
     print("Pre-process Start!!!")
+    print(f"Using Wallis Filter Width: {obj.WallisFilterWidth}")
 #    obj.zeroMask = 1
 
     # TODO: Allow different filters to be applied images independently
@@ -517,6 +521,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search
         m_name = os.path.basename(indir_m)
         s_name = os.path.basename(indir_s)
 
+        # FIXME: Filter width is a magic variable here and not exposed well.
+        preprocessing_filter_width = 5
+        if nc_sensor == 'S1':
+            preprocessing_filter_width = 21
+
+        print(f'Preprocessing filter width {preprocessing_filter_width}')
+
         preprocessing_methods = ['hps', 'hps']
         for ii, name in enumerate((m_name, s_name)):
             if len(re.findall("L[EO]07_", name)) > 0:
@@ -528,9 +539,11 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search
 
         print(f'Using preprocessing methods {preprocessing_methods}')
 
-        Dx, Dy, InterpMask, ChipSizeX, GridSpacingX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask = runAutorift(
-            data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0,
-            noDataMask, optical_flag, nodata, mpflag, geogrid_run_info=geogrid_run_info, preprocessing_methods=preprocessing_methods,
+        Dx, Dy, InterpMask, ChipSizeX, GridSpacingX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask = \
+            runAutorift(
+                data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0,
+                noDataMask, optical_flag, nodata, mpflag, geogrid_run_info=geogrid_run_info,
+                preprocessing_methods=preprocessing_methods, preprocessing_filter_width=preprocessing_filter_width,
         )
         if nc_sensor is not None:
             import netcdf_output as no
@@ -617,7 +630,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search
 
 
         if offset2vx is not None:
-        
+
             ds = gdal.Open(scale_factor)
             band = ds.GetRasterBand(1)
             scale_factor_1 = band.ReadAsArray()