CUDA interpolated geocode (#760)

* initial commit of CUDA geocode; interp mode only

* test code for pybind bindings

* moved to masked minmax to own files

* added explicit instantiations for thrust::complex

* added masked minmax to index calcuation and interpolation

* made line_start class member and tweaked mask value assignment

* making line_per_block adjustable

* added bindings for geocoding, n_blocks+lines_per_block getters, and documentation

* index and mask device vectors resized according to block size

* fix max_compare and add comments

* add missing cudaMalloc/Free for **interp and some reshuffling

* added 1st line/pixel offsets, extra margin, and row/col index swap to interp

* added projection as a class member to replace insufficient inverse call in index kernal

* added class for common params and block processing testing

* templated geocode and adopted templated pybind and condensed test run

* Properly populated pybind docstrings

* made read only data in kernel const

* added wavelength getter for consistency

* fixed cuda memcpy dst pointers, bicubic explicit instantiation, and clang-format

* one Rdr2Geo/Geo2RdrParams to rule them all

* added handle for CUDA interpolator and projection base and incorporated into CUDA geocode

* add user defined invalid pixel value; default 0.0

* expose geo2rdr parameters in constructor

* InterpolatorHandle constructor constructs full object

* introduced dem interp method

* added unit test for masked minmax

* added cuda::geocode to isce3 package and updated unit test accordingly

* Reset data interpolation method back to bilinear

Co-authored-by: Liang Yu <[email protected]>
Co-authored-by: Geoffrey M Gunter <[email protected]>

Author: 3 people

cxx/isce3/container/RadarGeometry.h

@@ -54,6 +54,9 @@ class RadarGeometry {
     /** Get radar look side */
     LookSide lookSide() const { return radarGrid().lookSide(); }
 
+    /** Get radar wavelength */
+    double wavelength() const { return radarGrid().wavelength(); }
+
 private:
     RadarGridParameters _radar_grid;
     Orbit _orbit;

cxx/isce3/cuda/Sources.cmake

@@ -9,11 +9,15 @@ core/gpuPoly2d.cu
 core/gpuProjections.cu
 core/gpuSinc2dInterpolator.cu
 core/gpuSpline2dInterpolator.cu
+core/InterpolatorHandle.cu
 core/Orbit.cu
+core/ProjectionBaseHandle.cu
 core/Stream.cu
 except/Error.cpp
 focus/Backproject.cu
 fft/detail/CufftWrapper.cu
+geocode/Geocode.cu
+geocode/MaskedMinMax.cu
 geometry/Geo2rdr.cpp
 geometry/gpuDEMInterpolator.cu
 geometry/gpuGeo2rdr.cu

cxx/isce3/cuda/container/RadarGeometry.h

@@ -73,6 +73,9 @@ class RadarGeometry {
     /** Get radar look side */
     LookSide lookSide() const { return radarGrid().lookSide(); }
 
+    /** Get radar wavelength */
+    double wavelength() const { return radarGrid().wavelength(); }
+
 private:
     RadarGridParameters _radar_grid;
     DeviceOrbit _orbit;

cxx/isce3/cuda/core/InterpolatorHandle.cu

@@ -0,0 +1,88 @@
+#include "InterpolatorHandle.h"
+
+#include <isce3/cuda/except/Error.h>
+#include <isce3/except/Error.h>
+
+#include <pyre/journal.h>
+
+#include <cuda_runtime.h>
+#include <thrust/complex.h>
+#include <thrust/device_vector.h>
+
+template<class T>
+using DeviceInterp = isce3::cuda::core::gpuInterpolator<T>;
+
+namespace isce3::cuda::core {
+
+template<class T>
+__global__ void init_interp(
+        DeviceInterp<T>** interp, isce3::core::dataInterpMethod interp_method,
+        bool * unsupported_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;
+        }
+    }
+}
+
+template<class T>
+__global__ void finalize_interp(DeviceInterp<T>** interp)
+{
+    if (threadIdx.x == 0 && blockIdx.x == 0) {
+        delete *interp;
+    }
+}
+
+template<class T>
+InterpolatorHandle<T>::InterpolatorHandle(
+        isce3::core::dataInterpMethod interp_method)
+{
+    checkCudaErrors(cudaMalloc(&_interp, sizeof(DeviceInterp<T>**)));
+
+    thrust::device_vector<bool> d_unsupported_interp(1, false);
+    init_interp<<<1, 1>>>(_interp, interp_method, d_unsupported_interp.data().get());
+    checkCudaErrors(cudaPeekAtLastError());
+    checkCudaErrors(cudaDeviceSynchronize());
+
+    bool unsupported_interp = d_unsupported_interp[0];
+    if (unsupported_interp)
+    {
+        pyre::journal::error_t error(
+                "isce.cuda.core.InterpolatorHandle.InterpolatorHandle");
+        error << "Unsupported interpolator method provided."
+              << pyre::journal::endl;
+        throw isce3::except::InvalidArgument(ISCE_SRCINFO(),
+                                         "Unsupported interpolator method provided.");
+    }
+}
+
+template<class T>
+InterpolatorHandle<T>::~InterpolatorHandle()
+{
+    finalize_interp<<<1, 1>>>(_interp);
+    checkCudaErrors(cudaPeekAtLastError());
+    checkCudaErrors(cudaDeviceSynchronize());
+    checkCudaErrors(cudaFree(_interp));
+}
+
+template class InterpolatorHandle<float>;
+template class InterpolatorHandle<thrust::complex<float>>;
+template class InterpolatorHandle<double>;
+template class InterpolatorHandle<thrust::complex<double>>;
+} // namespace isce3::cuda::core

cxx/isce3/cuda/core/InterpolatorHandle.h

@@ -0,0 +1,40 @@
+#pragma once
+
+#include <isce3/core/Constants.h>
+#include <isce3/cuda/core/gpuInterpolator.h>
+
+namespace isce3::cuda::core {
+
+/** Class that handles device gpuInterpolator double pointers on device.
+ *
+ * This handle class ensures that pointers are properly allocated
+ * and deallocated.
+ *
+ */
+template<class T>
+class InterpolatorHandle {
+private:
+    // double pointer to gpuInterpolator on device
+    // 1st pointer is the gpuInterpolator location on device
+    // 2nd pointer is the gpuInterpolator object on device
+    isce3::cuda::core::gpuInterpolator<T>** _interp = nullptr;
+
+public:
+    /** Class constructor. Mallocs 1st pointer and creates gpuInterpolator
+     * object on device.
+     */
+    InterpolatorHandle(isce3::core::dataInterpMethod interp_method);
+
+    /** Destructor that frees and deletes pointers accordingly. */
+    ~InterpolatorHandle();
+
+    /** Disabling copy constructor and assignment operator to prever misuse */
+    InterpolatorHandle(const InterpolatorHandle&) = delete;
+    InterpolatorHandle& operator=(const InterpolatorHandle&) = delete;
+
+    isce3::cuda::core::gpuInterpolator<T>** getInterp() const
+    {
+        return _interp;
+    };
+};
+} // namespace isce3::cuda::core

cxx/isce3/cuda/core/ProjectionBaseHandle.cu

@@ -0,0 +1,63 @@
+#include "ProjectionBaseHandle.h"
+
+#include <isce3/cuda/except/Error.h>
+#include <isce3/except/Error.h>
+
+#include <pyre/journal.h>
+
+#include <cuda_runtime.h>
+#include <thrust/device_vector.h>
+
+namespace isce3::cuda::core {
+
+using isce3::cuda::core::ProjectionBase;
+
+__global__ void init_proj(ProjectionBase** proj, int epsg_code,
+                          bool * proj_invalid)
+{
+
+    if (threadIdx.x == 0 && blockIdx.x == 0) {
+        (*proj) = isce3::cuda::core::createProj(epsg_code);
+        if (!*proj)
+            *proj_invalid = true;
+    }
+}
+
+__global__ void finalize_proj(ProjectionBase** proj)
+{
+
+    if (threadIdx.x == 0 && blockIdx.x == 0) {
+        delete *proj;
+    }
+}
+
+ProjectionBaseHandle::ProjectionBaseHandle(int epsg)
+{
+    checkCudaErrors(cudaMalloc(&_proj, sizeof(ProjectionBase**)));
+
+    thrust::device_vector<bool> d_proj_invalid(1, false);
+    init_proj<<<1, 1>>>(_proj, epsg, d_proj_invalid.data().get());
+    checkCudaErrors(cudaPeekAtLastError());
+    checkCudaErrors(cudaDeviceSynchronize());
+
+    bool proj_invalid = d_proj_invalid[0];
+    if (proj_invalid)
+    {
+        pyre::journal::error_t error(
+                "isce.cuda.core.ProjectionBaseHandle.ProjectionBaseHandle");
+        error << "Unsupported EPSG provided."
+              << pyre::journal::endl;
+        throw isce3::except::InvalidArgument(ISCE_SRCINFO(),
+                                         "Unsupported ESPG provided.");
+    }
+}
+
+ProjectionBaseHandle::~ProjectionBaseHandle()
+{
+    finalize_proj<<<1, 1>>>(_proj);
+    checkCudaErrors(cudaPeekAtLastError());
+    checkCudaErrors(cudaDeviceSynchronize());
+    checkCudaErrors(cudaFree(_proj));
+}
+
+} // namespace isce3::cuda::core

cxx/isce3/cuda/core/ProjectionBaseHandle.h

@@ -0,0 +1,38 @@
+#pragma once
+
+#include <isce3/cuda/core/gpuProjections.h>
+
+namespace isce3::cuda::core {
+
+/** Class that handles device ProjectionBase double pointers on device.
+ *
+ * This handle class ensures that pointers are properly allocated
+ * and deallocated.
+ *
+ */
+class ProjectionBaseHandle {
+private:
+    // double pointer to Projection Base on device
+    // 1st pointer is the ProjectionBase location on device
+    // 2nd pointer is the ProjectionBase object on device
+    isce3::cuda::core::ProjectionBase** _proj = nullptr;
+
+public:
+    /** Class constructor. Mallocs 1st pointer and creates ProjectionBase
+     *  object on device.
+     *
+     * \param[in] epsg  EPSG of ProjectionBase to be created
+     *
+     * */
+    ProjectionBaseHandle(int epsg);
+
+    /** Destructor that frees and deletes pointers accordingly. */
+    ~ProjectionBaseHandle();
+
+    /** Disabling copy constructor and assignment operator to prever misuse */
+    ProjectionBaseHandle(const ProjectionBaseHandle&) = delete;
+    ProjectionBaseHandle& operator=(const ProjectionBaseHandle&) = delete;
+
+    isce3::cuda::core::ProjectionBase** get_proj() const { return _proj; }
+};
+} // namespace isce3::cuda::core