Fix just geocode insar block processing (#955)

* added CUDA geocoding block processing done on C++

* pybind for block processing C++ function

* CPU geocode_insar discontinuity errors fixed
CUDA geocode_insar less nested
created shared input/output raster creator

Author: Liang Yu

cxx/isce3/cuda/geocode/Geocode.cu

@@ -5,6 +5,7 @@
 #include <thrust/functional.h>
 #include <thrust/host_vector.h>
 
+#include <gdal_priv.h>
 #include <pyre/journal.h>
 
 #include <isce3/core/Ellipsoid.h>
@@ -414,6 +415,57 @@ void Geocode::geocodeRasterBlock(Raster& output_raster, Raster& input_raster)
             _geogrid.width(), _geo_block_length, 1);
 }
 
+void Geocode::geocodeRasters(
+        std::vector<std::reference_wrapper<isce3::io::Raster>> output_rasters,
+        std::vector<std::reference_wrapper<isce3::io::Raster>> input_rasters)
+{
+    // check if vectors are of same length
+    if (output_rasters.size() != input_rasters.size()) {
+        throw isce3::except::LengthError(
+                ISCE_SRCINFO(), "number of input and output rasters not equal");
+    }
+
+    const auto n_raster_pairs = output_rasters.size();
+
+    // iterate over blocks
+    for (size_t i_block = 0; i_block < _n_blocks; ++i_block) {
+
+        // set radar coords for each geocode obj for curret block
+        setBlockRdrCoordGrid(i_block);
+
+        for (size_t i_raster = 0; i_raster < n_raster_pairs; ++i_raster)
+        {
+            const int dtype = input_rasters[i_raster].get().dtype();
+            switch (dtype) {
+                case GDT_Float32:   {
+                    geocodeRasterBlock<float>(
+                            output_rasters[i_raster], input_rasters[i_raster]);
+                    break; }
+                case GDT_CFloat32:  {
+                    geocodeRasterBlock<thrust::complex<float>>(
+                            output_rasters[i_raster], input_rasters[i_raster]);
+                    break;}
+                case GDT_Float64:   {
+                    geocodeRasterBlock<double>(
+                            output_rasters[i_raster], input_rasters[i_raster]);
+                    break; }
+                case GDT_CFloat64:  {
+                    geocodeRasterBlock<thrust::complex<double>>(
+                            output_rasters[i_raster], input_rasters[i_raster]);
+                    break;}
+                case GDT_Byte:  {
+                    geocodeRasterBlock<unsigned char>(
+                            output_rasters[i_raster], input_rasters[i_raster]);
+                    break;}
+                default: {
+                    throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                            "unsupported datatype");
+                         }
+            }
+        }
+    }
+}
+
 #define EXPLICIT_INSTATIATION(T)                                               \
     template void Geocode::geocodeRasterBlock<T>(                              \
             Raster & output_raster, Raster & input_raster);

cxx/isce3/cuda/geocode/Geocode.h

@@ -1,5 +1,8 @@
 #pragma once
 
+#include <functional>
+#include <vector>
+
 #include <isce3/core/forward.h>
 #include <isce3/geometry/forward.h>
 
@@ -104,6 +107,16 @@ class Geocode {
     void geocodeRasterBlock(
             isce3::io::Raster& output_raster, isce3::io::Raster& input_raster);
 
+    /** Geocode rasters with a shared geogrid. Block processing handled
+     * internally in function.
+     *
+     * \param[in] output_rasters    Geocoded rasters
+     * \param[in] input_rasters     Rasters to be geocoded
+     */
+    void geocodeRasters(
+            std::vector<std::reference_wrapper<isce3::io::Raster>> output_rasters,
+            std::vector<std::reference_wrapper<isce3::io::Raster>> input_rasters);
+
     size_t numBlocks() const { return _n_blocks; }
     size_t linesPerBlock() const { return _lines_per_block; }
 

python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp

@@ -1,7 +1,7 @@
 #include "cuGeocode.h"
-#include <isce3/core/forward.h>
 
 #include <gdal_priv.h>
+#include <pybind11/stl.h>
 
 #include <isce3/container/RadarGeometry.h>
 #include <isce3/core/Constants.h>
@@ -119,7 +119,20 @@ void addbinding(pybind11::class_<Geocode>& pyGeocode)
             input_raster: io::Raster
                 Raster to be geocoded
         )")
-            .def_property_readonly("n_blocks", &Geocode::numBlocks)
-            .def_property_readonly("lines_per_block", &Geocode::linesPerBlock);
+    .def("geocode_rasters", &Geocode::geocodeRasters,
+            py::arg("output_rasters"),
+            py::arg("input_rasters"),
+            R"(
+            Geocode rasters with a shared geogrid with block processing handled internally.
+
+            Parameters
+            ----------
+            output_rasters: list(io::Raster)
+                List of geocoded rasters.
+            input_rasters: list(io::Raster)
+                List of rasters to be geocoded.
+        )")
+    .def_property_readonly("n_blocks", &Geocode::numBlocks)
+    .def_property_readonly("lines_per_block", &Geocode::linesPerBlock);
     ;
 }