Additonal block margin for GSLC and CUDA geocode (#833)

* update block margin
* revert changes brought by mistake from PR #783
* Update cxx/isce3/geocode/GeocodeCov.cpp
* added extra margin for GSLC block interp to address gaps in output block boundaries. changes mirror PR 830
* port extra margin fix to CUDA geocode. std::min/max replaced to account for unsigned size_t

Co-authored-by: Gustavo Shiroma <[email protected]>
Co-authored-by: Liang Yu <[email protected]>

Author: 3 people

cxx/isce3/cuda/geocode/Geocode.cu

@@ -24,6 +24,7 @@
 #include "Geocode.h"
 #include "MaskedMinMax.h"
 
+#include <cstdio>
 using isce3::core::Vec3;
 using isce3::io::Raster;
 
@@ -307,6 +308,16 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number)
     _az_last_line = std::max(static_cast<size_t>(0),
             static_cast<size_t>(std::ceil(rdr_y_max) - 1));
 
+    // Extra margin for interpolation to avoid gaps between blocks in output
+    const size_t interp_margin {5};
+    _range_first_pixel = interp_margin > _range_first_pixel ? 0 : _range_first_pixel - interp_margin;
+    _range_last_pixel = std::min(_rdr_geom.radarGrid().width() - 1,
+                             _range_last_pixel + interp_margin);
+
+    _az_first_line = interp_margin > _az_first_line ? 0 : _az_first_line - interp_margin;
+    _az_last_line = std::min(_rdr_geom.radarGrid().length() - 1,
+                             _az_last_line+interp_margin);
+
     // check if block entirely masked
     bool all_masked = std::isnan(rdr_y_min);
 

cxx/isce3/geocode/GeocodeCov.cpp

@@ -432,14 +432,6 @@ void Geocode<T>::geocodeInterp(
             out_geo_dem_array.fill(std::numeric_limits<float>::quiet_NaN());
         }
 
-        // First and last line of the data block in radar coordinates
-        int azimuthFirstLine = radar_grid.length() - 1;
-        int azimuthLastLine = 0;
-
-        // First and last pixel of the data block in radar coordinates
-        int rangeFirstPixel = radar_grid.width() - 1;
-        int rangeLastPixel = 0;
-
         // load a block of DEM for the current geocoded grid
         isce3::geometry::DEMInterpolator demInterp = isce3::geocode::loadDEM(
                 demRaster, geogrid, lineStart, geoBlockLength, geogrid.width(),
@@ -450,17 +442,17 @@ void Geocode<T>::geocodeInterp(
         std::valarray<double> radarX(blockSize);
         std::valarray<double> radarY(blockSize);
 
-        int localAzimuthFirstLine = radar_grid.length() - 1;
-        int localAzimuthLastLine = 0;
-        int localRangeFirstPixel = radar_grid.width() - 1;
-        int localRangeLastPixel = 0;
+        int azimuthFirstLine = radar_grid.length() - 1;
+        int azimuthLastLine = 0;
+        int rangeFirstPixel = radar_grid.width() - 1;
+        int rangeLastPixel = 0;
 
         // Loop over lines, samples of the output grid
 #pragma omp parallel for reduction(                                            \
         min                                                                    \
-        : localAzimuthFirstLine, localRangeFirstPixel)                         \
+        : azimuthFirstLine, rangeFirstPixel)                         \
         reduction(max                                                          \
-                  : localAzimuthLastLine, localRangeLastPixel)
+                  : azimuthLastLine, rangeLastPixel)
 
         for (size_t kk = 0; kk < geoBlockLength * geogrid.width(); ++kk) {
 
@@ -525,14 +517,14 @@ void Geocode<T>::geocodeInterp(
                     rdrX >= radar_grid.width())
                 continue;
 
-            localAzimuthFirstLine = std::min(
-                    localAzimuthFirstLine, static_cast<int>(std::floor(rdrY)));
-            localAzimuthLastLine = std::max(localAzimuthLastLine,
+            azimuthFirstLine = std::min(
+                    azimuthFirstLine, static_cast<int>(std::floor(rdrY)));
+            azimuthLastLine = std::max(azimuthLastLine,
                     static_cast<int>(std::ceil(rdrY) - 1));
-            localRangeFirstPixel = std::min(
-                    localRangeFirstPixel, static_cast<int>(std::floor(rdrX)));
-            localRangeLastPixel = std::max(
-                    localRangeLastPixel, static_cast<int>(std::ceil(rdrX) - 1));
+            rangeFirstPixel = std::min(
+                    rangeFirstPixel, static_cast<int>(std::floor(rdrX)));
+            rangeLastPixel = std::max(
+                    rangeLastPixel, static_cast<int>(std::ceil(rdrX) - 1));
 
             // store the adjusted X and Y indices
             radarX[blockLine * geogrid.width() + pixel] = rdrX;
@@ -558,11 +550,15 @@ void Geocode<T>::geocodeInterp(
                     geogrid.width(), geoBlockLength, 1);
         }
 
-        // Get min and max swath extents from among all threads
-        azimuthFirstLine = std::min(azimuthFirstLine, localAzimuthFirstLine);
-        azimuthLastLine = std::max(azimuthLastLine, localAzimuthLastLine);
-        rangeFirstPixel = std::min(rangeFirstPixel, localRangeFirstPixel);
-        rangeLastPixel = std::max(rangeLastPixel, localRangeLastPixel);
+        // Add extra margin for interpolation
+        int interp_margin = 5;
+        azimuthFirstLine = std::max(azimuthFirstLine - interp_margin, 0);
+        rangeFirstPixel = std::max(rangeFirstPixel - interp_margin, 0);
+
+        azimuthLastLine = std::min(azimuthLastLine + interp_margin, 
+                                   static_cast<int>(radar_grid.length() - 1));
+        rangeLastPixel = std::min(rangeLastPixel + interp_margin, 
+                                  static_cast<int>(radar_grid.width() - 1));
 
         if (azimuthFirstLine > azimuthLastLine ||
                 rangeFirstPixel > rangeLastPixel)
@@ -735,7 +731,8 @@ inline void Geocode<T>::_interpolate(
 
     size_t length = geoDataBlock.length();
     size_t width = geoDataBlock.width();
-    int extraMargin = isce3::core::SINC_HALF;
+    // Add extra margin for interpolation
+    int interp_margin = 5;
 
     double offsetY =
             azimuthFirstLine / radar_grid.prf() + radar_grid.sensingStart();
@@ -753,9 +750,9 @@ inline void Geocode<T>::_interpolate(
         double rdrY = radarY[i * width + j] - azimuthFirstLine;
         double rdrX = radarX[i * width + j] - rangeFirstPixel;
 
-        if (rdrX < extraMargin || rdrY < extraMargin ||
-                rdrX >= (radarBlockWidth - extraMargin) ||
-                rdrY >= (radarBlockLength - extraMargin)) {
+        if (rdrX < interp_margin || rdrY < interp_margin ||
+                rdrX >= (radarBlockWidth - interp_margin) ||
+                rdrY >= (radarBlockLength - interp_margin)) {
 
             // set NaN values according to T_out, i.e. real (NaN) or complex
             // (NaN, NaN)

cxx/isce3/geocode/geocodeSlc.cpp

@@ -59,14 +59,6 @@ void isce3::geocode::geocodeSlc(
         }
         size_t blockSize = geoBlockLength * geoGrid.width();
 
-        // First and last line of the data block in radar coordinates
-        int azimuthFirstLine = radarGrid.length() - 1;
-        int azimuthLastLine = 0;
-
-        // First and last pixel of the data block in radar coordinates
-        int rangeFirstPixel = radarGrid.width() - 1;
-        int rangeLastPixel = 0;
-
         // get a DEM interpolator for a block of DEM for the current geocoded
         // grid
         isce3::geometry::DEMInterpolator demInterp = isce3::geocode::loadDEM(
@@ -78,18 +70,20 @@ void isce3::geocode::geocodeSlc(
         std::valarray<double> radarX(blockSize);
         std::valarray<double> radarY(blockSize);
 
-        int localAzimuthFirstLine = radarGrid.length() - 1;
-        int localAzimuthLastLine = 0;
-        int localRangeFirstPixel = radarGrid.width() - 1;
-        int localRangeLastPixel = 0;
+        // First and last line of the data block in radar coordinates
+        int azimuthFirstLine = radarGrid.length() - 1;
+        int azimuthLastLine = 0;
+        int rangeFirstPixel = radarGrid.width() - 1;
+        int rangeLastPixel = 0;
+        // First and last pixel of the data block in radar coordinates
 
         size_t geoGridWidth = geoGrid.width();
 // Loop over lines, samples of the output grid
-#pragma omp parallel for reduction(min                                         \
-                                   : localAzimuthFirstLine,                    \
-                                     localRangeFirstPixel)                     \
-        reduction(max                                                          \
-                  : localAzimuthLastLine, localRangeLastPixel)
+#pragma omp parallel for reduction(min                                    \
+                                   : azimuthFirstLine,                    \
+                                     rangeFirstPixel)                     \
+        reduction(max                                                     \
+                  : azimuthLastLine, rangeLastPixel)
         for (size_t kk = 0; kk < geoBlockLength * geoGridWidth; ++kk) {
 
             size_t blockLine = kk / geoGridWidth;
@@ -98,8 +92,8 @@ void isce3::geocode::geocodeSlc(
             // Global line index
             const size_t line = lineStart + blockLine;
 
-            // y coordinate in the out put grid 
-	    // Assuming geoGrid.startY() and geoGrid.startX() represent the top-left 
+            // y coordinate in the out put grid
+	    // Assuming geoGrid.startY() and geoGrid.startX() represent the top-left
 	    // corner of the first pixel, then 0.5 pixel shift is needed to get
 	    // to the center of each pixel
             double y = geoGrid.startY() + geoGrid.spacingY() * (line + 0.5);
@@ -143,27 +137,33 @@ void isce3::geocode::geocodeSlc(
                 not nativeDoppler.contains(aztime, srange))
                 continue;
 
-            localAzimuthFirstLine = std::min(
-                    localAzimuthFirstLine, static_cast<int>(std::floor(rdrY)));
-            localAzimuthLastLine =
-                    std::max(localAzimuthLastLine,
+            azimuthFirstLine = std::min(
+                    azimuthFirstLine, static_cast<int>(std::floor(rdrY)));
+            azimuthLastLine =
+                    std::max(azimuthLastLine,
                              static_cast<int>(std::ceil(rdrY) - 1));
-            localRangeFirstPixel = std::min(localRangeFirstPixel,
+            rangeFirstPixel = std::min(rangeFirstPixel,
                                             static_cast<int>(std::floor(rdrX)));
-            localRangeLastPixel = std::max(
-                    localRangeLastPixel, static_cast<int>(std::ceil(rdrX) - 1));
+            rangeLastPixel = std::max(
+                    rangeLastPixel, static_cast<int>(std::ceil(rdrX) - 1));
 
             // store the adjusted X and Y indices
             radarX[blockLine * geoGrid.width() + pixel] = rdrX;
             radarY[blockLine * geoGrid.width() + pixel] = rdrY;
 
         } // end loops over lines and pixel of output grid
 
+        // Extra margin for interpolation to avoid gaps between blocks in output
+        int interp_margin = 5;
+
         // Get min and max swath extents from among all threads
-        azimuthFirstLine = std::min(azimuthFirstLine, localAzimuthFirstLine);
-        azimuthLastLine = std::max(azimuthLastLine, localAzimuthLastLine);
-        rangeFirstPixel = std::min(rangeFirstPixel, localRangeFirstPixel);
-        rangeLastPixel = std::max(rangeLastPixel, localRangeLastPixel);
+        azimuthFirstLine = std::max(azimuthFirstLine - interp_margin, 0);
+        rangeFirstPixel = std::max(rangeFirstPixel - interp_margin, 0);
+
+        azimuthLastLine = std::min(azimuthLastLine + interp_margin,
+                                   static_cast<int>(radarGrid.length() - 1));
+        rangeLastPixel = std::min(rangeLastPixel + interp_margin,
+                                  static_cast<int>(radarGrid.width() - 1));
 
         if (azimuthFirstLine > azimuthLastLine ||
             rangeFirstPixel > rangeLastPixel)