CUDA geocode for InSAR workflow (#820)

Liang Yu · GitHub Enterprise · commit 758a1a9f5ee9 · 2021-07-01T18:03:45.000-07:00
* CUDA nearest neighbor interp2d commit 
* added unsigned char support to CUDA interp2d
* initial commit of CUDA geocode workflow
CUDA and CPU geocode split into 2 functions due differences in block processing
* temporary fix for CUDA GUNW's inability to directly write to HDF5 with chunks
diff --git a/cxx/isce3/cuda/Sources.cmake b/cxx/isce3/cuda/Sources.cmake
@@ -5,6 +5,7 @@ core/gpuBicubicInterpolator.cu
 core/gpuBilinearInterpolator.cu
 core/gpuLUT1d.cu
 core/gpuLUT2d.cu
+core/gpuNearestNeighborInterpolator.cu
 core/gpuPoly2d.cu
 core/gpuProjections.cu
 core/gpuSinc2dInterpolator.cu
diff --git a/cxx/isce3/cuda/core/InterpolatorHandle.cu b/cxx/isce3/cuda/core/InterpolatorHandle.cu
@@ -20,20 +20,23 @@ __global__ void init_interp(DeviceInterp<T>** interp,
 {
     if (threadIdx.x == 0 && blockIdx.x == 0) {
         // Choose interpolator
-        switch (interp_method) {
-        case isce3::core::BILINEAR_METHOD:
-            (*interp) = new isce3::cuda::core::gpuBilinearInterpolator<T>();
-            break;
-        case isce3::core::BICUBIC_METHOD:
-            (*interp) = new isce3::cuda::core::gpuBicubicInterpolator<T>();
-            break;
-        case isce3::core::BIQUINTIC_METHOD: {
-            size_t order = 6;
-            (*interp) =
-                    new isce3::cuda::core::gpuSpline2dInterpolator<T>(order);
-            break;
-        }
-        default: *unsupported_interp = true; break;
+        switch(interp_method) {
+            case isce3::core::BILINEAR_METHOD:
+                (*interp) = new isce3::cuda::core::gpuBilinearInterpolator<T>();
+                break;
+            case isce3::core::BICUBIC_METHOD:
+                (*interp) = new isce3::cuda::core::gpuBicubicInterpolator<T>();
+                break;
+            case isce3::core::BIQUINTIC_METHOD:
+                size_t order = 6;
+                (*interp) = new isce3::cuda::core::gpuSpline2dInterpolator<T>(order);
+                break;
+            case isce3::core::NEAREST_METHOD:
+                (*interp) = new isce3::cuda::core::gpuNearestNeighborInterpolator<T>();
+                break;
+            default:
+                *unsupported_interp = true;
+                break;
         }
     }
 }
@@ -82,4 +85,5 @@ template class InterpolatorHandle<float>;
 template class InterpolatorHandle<thrust::complex<float>>;
 template class InterpolatorHandle<double>;
 template class InterpolatorHandle<thrust::complex<double>>;
+template class InterpolatorHandle<unsigned char>;
 } // namespace isce3::cuda::core
diff --git a/cxx/isce3/cuda/core/gpuBicubicInterpolator.cu b/cxx/isce3/cuda/core/gpuBicubicInterpolator.cu
@@ -178,6 +178,7 @@ template class gpuBicubicInterpolator<double>;
 template class gpuBicubicInterpolator<thrust::complex<double>>;
 template class gpuBicubicInterpolator<float>;
 template class gpuBicubicInterpolator<thrust::complex<float>>;
+template class gpuBicubicInterpolator<unsigned char>;
 
 template __global__ void gpuInterpolator_g<double>(
         gpuBicubicInterpolator<double> interp, double* x, double* y,
diff --git a/cxx/isce3/cuda/core/gpuBilinearInterpolator.cu b/cxx/isce3/cuda/core/gpuBilinearInterpolator.cu
@@ -114,6 +114,7 @@ template class gpuBilinearInterpolator<double>;
 template class gpuBilinearInterpolator<thrust::complex<double>>;
 template class gpuBilinearInterpolator<float>;
 template class gpuBilinearInterpolator<thrust::complex<float>>;
+template class gpuBilinearInterpolator<unsigned char>;
 
 template __global__ void gpuInterpolator_g<double>(
         gpuBilinearInterpolator<double> interp, double* x, double* y,
diff --git a/cxx/isce3/cuda/core/gpuInterpolator.h b/cxx/isce3/cuda/core/gpuInterpolator.h
@@ -75,4 +75,13 @@ class gpuSinc2dInterpolator : public isce3::cuda::core::gpuInterpolator<U> {
         CUDA_DEV U interpolate(double, double, const U*, size_t, size_t);
         CUDA_HOST void interpolate_h(const Matrix<double>&, Matrix<U>&, double, double, U*);
 };
+
+/** gpuNearestNeighborInterpolator class derived from abstract gpuInterpolator class */
+template <class T>
+class gpuNearestNeighborInterpolator : public isce3::cuda::core::gpuInterpolator<T> {
+    public:
+        CUDA_HOSTDEV gpuNearestNeighborInterpolator(){};
+        CUDA_DEV T interpolate(double x, double y, const T* z, size_t nx, size_t ny = 0);
+};
+
 }}}
diff --git a/cxx/isce3/cuda/core/gpuNearestNeighborInterpolator.cu b/cxx/isce3/cuda/core/gpuNearestNeighborInterpolator.cu
@@ -0,0 +1,21 @@
+#include "gpuInterpolator.h"
+
+#include <thrust/complex.h>
+
+using isce3::cuda::core::gpuNearestNeighborInterpolator;
+
+template<class T>
+__device__ T gpuNearestNeighborInterpolator<T>::interpolate(
+        double x, double y, const T* z, size_t nx, size_t /*ny*/)
+{
+    size_t x_round = static_cast<size_t>(std::round(x));
+    size_t y_round = static_cast<size_t>(std::round(y));
+
+    return z[y_round * nx + x_round];
+}
+
+template class gpuNearestNeighborInterpolator<double>;
+template class gpuNearestNeighborInterpolator<thrust::complex<double>>;
+template class gpuNearestNeighborInterpolator<float>;
+template class gpuNearestNeighborInterpolator<thrust::complex<float>>;
+template class gpuNearestNeighborInterpolator<unsigned char>;
diff --git a/cxx/isce3/cuda/core/gpuSpline2dInterpolator.cu b/cxx/isce3/cuda/core/gpuSpline2dInterpolator.cu
@@ -153,6 +153,7 @@ template class gpuSpline2dInterpolator<double>;
 template class gpuSpline2dInterpolator<thrust::complex<double>>;
 template class gpuSpline2dInterpolator<float>;
 template class gpuSpline2dInterpolator<thrust::complex<float>>;
+template class gpuSpline2dInterpolator<unsigned char>;
 
 template __global__ void gpuInterpolator_g<double>(
         gpuSpline2dInterpolator<double> interp, double* x, double* y,
diff --git a/cxx/isce3/cuda/geocode/Geocode.cu b/cxx/isce3/cuda/geocode/Geocode.cu
@@ -172,26 +172,35 @@ __global__ void interpolate(T* geo_data_block, const double* __restrict__ rdr_x,
     geo_data_block[tid] = interp_val;
 }
 
-__host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid,
-        const isce3::container::RadarGeometry& rdr_geom,
-        const Raster& dem_raster, const double dem_margin,
-        const size_t lines_per_block,
-        const isce3::core::dataInterpMethod data_interp_method,
-        const isce3::core::dataInterpMethod dem_interp_method,
-        const double threshold, const int maxiter, const double dr,
-        const float invalid_value)
-    : _geogrid(geogrid), _rdr_geom(rdr_geom),
-      _ellipsoid(isce3::core::makeProjection(_geogrid.epsg())->ellipsoid()),
-      _lines_per_block(lines_per_block), _geo_block_length(_lines_per_block),
-      _n_blocks((geogrid.length() + _lines_per_block - 1) / _lines_per_block),
-      _az_first_line(_rdr_geom.radarGrid().length() - 1), _az_last_line(0),
-      _range_first_pixel(_rdr_geom.radarGrid().width() - 1),
-      _range_last_pixel(0), _dem_raster(dem_raster), _dem_margin(dem_margin),
-      _interp_float_handle(data_interp_method),
-      _interp_cfloat_handle(data_interp_method),
-      _interp_double_handle(data_interp_method),
-      _interp_cdouble_handle(data_interp_method), _proj_handle(geogrid.epsg()),
-      _dem_interp_method(dem_interp_method)
+__host__
+Geocode::Geocode(const isce3::product::GeoGridParameters & geogrid,
+                const isce3::container::RadarGeometry & rdr_geom,
+                const Raster & dem_raster,
+                const double dem_margin,
+                const size_t lines_per_block,
+                const isce3::core::dataInterpMethod data_interp_method,
+                const isce3::core::dataInterpMethod dem_interp_method,
+                const double threshold, const int maxiter,
+                const double dr, const float invalid_value) :
+    _geogrid(geogrid),
+    _rdr_geom(rdr_geom),
+    _ellipsoid(isce3::core::makeProjection(_geogrid.epsg())->ellipsoid()),
+    _lines_per_block(lines_per_block),
+    _geo_block_length(_lines_per_block),
+    _n_blocks((geogrid.length() + _lines_per_block -1) / _lines_per_block),
+    _az_first_line(_rdr_geom.radarGrid().length() - 1),
+    _az_last_line(0),
+    _range_first_pixel(_rdr_geom.radarGrid().width() - 1),
+    _range_last_pixel(0),
+    _dem_raster(dem_raster),
+    _dem_margin(dem_margin),
+    _interp_float_handle(data_interp_method),
+    _interp_cfloat_handle(data_interp_method),
+    _interp_double_handle(data_interp_method),
+    _interp_cdouble_handle(data_interp_method),
+    _interp_unsigned_char_handle(data_interp_method),
+    _proj_handle(geogrid.epsg()),
+    _dem_interp_method(dem_interp_method)
 {
     // init light weight radar grid
     _radar_grid.sensing_start = _rdr_geom.radarGrid().sensingStart();
@@ -219,9 +228,11 @@ __host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid,
     if (std::isnan(invalid_value)) {
         _invalid_float = std::numeric_limits<float>::quiet_NaN();
         _invalid_double = std::numeric_limits<double>::quiet_NaN();
+        _invalid_unsigned_char = 255;
     } else {
         _invalid_float = invalid_value;
-        _invalid_double = (double) invalid_value;
+        _invalid_double = static_cast<double>(invalid_value);
+        _invalid_unsigned_char = static_cast<unsigned char>(invalid_value);
     }
 }
 
@@ -334,8 +345,11 @@ void Geocode::geocodeRasterBlock(Raster& output_raster, Raster& input_raster)
         invalid_value = _invalid_double;
     } else if constexpr (std::is_same_v<T, thrust::complex<double>>) {
         interp = _interp_cdouble_handle.getInterp();
-        invalid_value =
-                thrust::complex<double>(_invalid_double, _invalid_double);
+        invalid_value = thrust::complex<double>(_invalid_double,
+                                                _invalid_double);
+    } else if constexpr (std::is_same_v<T, unsigned char>) {
+        interp = _interp_unsigned_char_handle.getInterp();
+        invalid_value = _invalid_unsigned_char;
     }
 
     // 0 width indicates current block is out of bounds
@@ -396,4 +410,5 @@ EXPLICIT_INSTATIATION(float);
 EXPLICIT_INSTATIATION(thrust::complex<float>);
 EXPLICIT_INSTATIATION(double);
 EXPLICIT_INSTATIATION(thrust::complex<double>);
+EXPLICIT_INSTATIATION(unsigned char);
 } // namespace isce3::cuda::geocode
diff --git a/cxx/isce3/cuda/geocode/Geocode.h b/cxx/isce3/cuda/geocode/Geocode.h
@@ -163,10 +163,13 @@ class Geocode {
     isce3::cuda::core::InterpolatorHandle<double> _interp_double_handle;
     isce3::cuda::core::InterpolatorHandle<thrust::complex<double>>
             _interp_cdouble_handle;
+    isce3::cuda::core::InterpolatorHandle<unsigned char>
+            _interp_unsigned_char_handle;
 
     // value applied to invalid geogrid pixels
     float _invalid_float;
     double _invalid_double;
+    unsigned char _invalid_unsigned_char;
 
     isce3::core::dataInterpMethod _dem_interp_method;
 };
diff --git a/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp b/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp
@@ -72,40 +72,40 @@ void addbinding(pybind11::class_<Geocode>& pyGeocode)
                 Index of block of raster where radar grid coordinates are calculated
                 and set.
             )")
-            .def(
-                    "geocode_raster_block",
-                    [](Geocode& self, isce3::io::Raster& output_raster,
-                            isce3::io::Raster& input_raster) {
-                        const int dtype = input_raster.dtype();
-                        switch (dtype) {
-                        case GDT_Float32: {
-                            self.geocodeRasterBlock<float>(
-                                    output_raster, input_raster);
-                            break;
-                        }
-                        case GDT_CFloat32: {
-                            self.geocodeRasterBlock<thrust::complex<float>>(
-                                    output_raster, input_raster);
-                            break;
-                        }
-                        case GDT_Float64: {
-                            self.geocodeRasterBlock<double>(
-                                    output_raster, input_raster);
-                            break;
-                        }
-                        case GDT_CFloat64: {
-                            self.geocodeRasterBlock<thrust::complex<double>>(
-                                    output_raster, input_raster);
-                            break;
-                        }
-                        default: {
-                            throw isce3::except::RuntimeError(
-                                    ISCE_SRCINFO(), "unsupported datatype");
-                        }
-                        }
-                    },
-                    py::arg("output_raster"), py::arg("input_raster"),
-                    R"(
+    .def("geocode_raster_block", [](Geocode & self,
+                                    isce3::io::Raster & output_raster,
+                                    isce3::io::Raster & input_raster) {
+                const int dtype =  input_raster.dtype();
+                switch (dtype) {
+                    case GDT_Float32:   {
+                        self.geocodeRasterBlock<float>(
+                                output_raster, input_raster);
+                        break; }
+                    case GDT_CFloat32:  {
+                        self.geocodeRasterBlock<thrust::complex<float>>(
+                                output_raster, input_raster);
+                        break;}
+                    case GDT_Float64:   {
+                        self.geocodeRasterBlock<double>(
+                                output_raster, input_raster);
+                        break; }
+                    case GDT_CFloat64:  {
+                        self.geocodeRasterBlock<thrust::complex<double>>(
+                                output_raster, input_raster);
+                        break;}
+                    case GDT_Byte:  {
+                        self.geocodeRasterBlock<unsigned char>(
+                                output_raster, input_raster);
+                        break;}
+                    default: {
+                        throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                                "unsupported datatype");
+                             }
+                }
+            },
+            py::arg("output_raster"),
+            py::arg("input_raster"),
+            R"(
             Geocode raster according to block specified in make_radar_grid_coordinates.
 
             Parameters
diff --git a/python/packages/pybind_nisar/workflows/geocode_insar.py b/python/packages/pybind_nisar/workflows/geocode_insar.py
diff --git a/python/packages/pybind_nisar/workflows/h5_prep.py b/python/packages/pybind_nisar/workflows/h5_prep.py
diff --git a/python/packages/pybind_nisar/workflows/insar_runconfig.py b/python/packages/pybind_nisar/workflows/insar_runconfig.py
diff --git a/tests/python/packages/pybind_nisar/workflows/geocode_insar.py b/tests/python/packages/pybind_nisar/workflows/geocode_insar.py
