resample SLC reference radar grid constructor, carrier indexing, and CUDA doppler fix (#969)

* untested: modify reference grid constructor
* carrier eval with az and range instead of index
* removed unnecessary modulos
* added docstrings

Author: Liang Yu

cxx/isce3/cuda/image/ResampSlc.cpp

@@ -1,10 +1,3 @@
-//-*- C++ -*-
-//-*- coding: utf-8 -*-
-//
-// Author: Joshua Cohen, Bryan V. Riel, Liang Yu
-// Copyright 2017-2018
-//
-
 #include "ResampSlc.h"
 
 #include <algorithm>
@@ -14,7 +7,7 @@
 
 // isce3::core
 #include <isce3/core/Constants.h>
-#include <isce3/core/LUT1d.h>
+#include <isce3/core/LUT2d.h>
 
 #include <isce3/image/Tile.h>
 
@@ -30,9 +23,9 @@ resamp(const std::string & inputFilename,          // filename of input SLC
        const std::string & outputFilename,         // filename of output resampled SLC
        const std::string & rgOffsetFilename,       // filename of range offsets
        const std::string & azOffsetFilename,       // filename of azimuth offsets
-       int inputBand, bool flatten, bool isComplex, int rowBuffer,
-       int chipSize) {
-
+       int inputBand, bool flatten, int rowBuffer,
+       int chipSize)
+{
     // Make input rasters
     Raster inputSlc(inputFilename, GA_ReadOnly);
     Raster rgOffsetRaster(rgOffsetFilename, GA_ReadOnly);
@@ -45,21 +38,16 @@ resamp(const std::string & inputFilename,          // filename of input SLC
 
     // Call generic resamp
     resamp(inputSlc, outputSlc, rgOffsetRaster, azOffsetRaster, inputBand, flatten,
-           isComplex, rowBuffer, chipSize);
+           rowBuffer, chipSize);
 }
 
 // Generic resamp entry point from externally created rasters
 void isce3::cuda::image::ResampSlc::
 resamp(isce3::io::Raster & inputSlc, isce3::io::Raster & outputSlc,
        isce3::io::Raster & rgOffsetRaster, isce3::io::Raster & azOffsetRaster,
-       int inputBand, bool flatten, bool isComplex, int rowBuffer,
-       int chipSize) {
-
-    // Check if data are not complex
-    if (!isComplex) {
-        std::cout << "Real data interpolation not implemented yet.\n";
-        return;
-    }
+       int inputBand, bool flatten, int rowBuffer,
+       int chipSize)
+{
     // Set the band number for input SLC
     _inputBand = inputBand;
     // Cache width of SLC image
@@ -70,17 +58,18 @@ resamp(isce3::io::Raster & inputSlc, isce3::io::Raster & outputSlc,
     const int outWidth = rgOffsetRaster.width();
 
     // Check if reference data is available
-    if (!this->haveRefData()) {
-        flatten = false;
+    if (flatten && !this->haveRefData()) {
+        std::string error_msg{"Unable to flatten; reference data not provided."};
+        throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg);
     }
 
     // initialize interpolator
     isce3::cuda::core::gpuSinc2dInterpolator<thrust::complex<float>> interp(chipSize-1, isce3::core::SINC_SUB);
 
     // Determine number of tiles needed to process image
     const int nTiles = _computeNumberOfTiles(outLength, _linesPerTile);
-    std::cout << 
-        "GPU resampling using " << nTiles << " tiles of " << _linesPerTile 
+    std::cout <<
+        "GPU resampling using " << nTiles << " tiles of " << _linesPerTile
         << " lines per tile\n";
     // Start timer
     auto timerStart = std::chrono::steady_clock::now();
@@ -106,15 +95,16 @@ resamp(isce3::io::Raster & inputSlc, isce3::io::Raster & outputSlc,
 
         // Get corresponding image indices
         std::cout << "Reading in image data for tile " << tileCount << std::endl;
-        _initializeTile(tile, inputSlc, azOffTile, outLength, rowBuffer, chipSize/2); 
+        _initializeTile(tile, inputSlc, azOffTile, outLength, rowBuffer, chipSize/2);
 
         // Perform interpolation
         std::cout << "Interpolating tile " << tileCount << std::endl;
-        gpuTransformTile(tile, outputSlc, rgOffTile, azOffTile, _rgCarrier, _azCarrier, 
-                isce3::core::avgLUT2dToLUT1d<double>(_dopplerLUT), interp, inWidth, inLength, this->startingRange(),
-                this->rangePixelSpacing(), this->prf(), this->wavelength(),
-                this->refStartingRange(), this->refRangePixelSpacing(),
-                this->refWavelength(), flatten, chipSize);
+        gpuTransformTile(tile, outputSlc, rgOffTile, azOffTile, _rgCarrier, _azCarrier,
+                _dopplerLUT, interp, inWidth, inLength, this->startingRange(),
+                this->rangePixelSpacing(), this->sensingStart(), this->prf(),
+                this->wavelength(), this->refStartingRange(),
+                this->refRangePixelSpacing(), this->refWavelength(), flatten,
+                chipSize, _invalid_value);
 
     }
 
@@ -124,5 +114,3 @@ resamp(isce3::io::Raster & inputSlc, isce3::io::Raster & outputSlc,
         timerEnd - timerStart).count();
     std::cout << "Elapsed processing time: " << elapsed << " sec" << "\n";
 }
-
-// end of file

cxx/isce3/cuda/image/ResampSlc.h

@@ -11,54 +11,84 @@ class isce3::cuda::image::ResampSlc : public isce3::image::ResampSlc {
     public:
         // Meta-methods
         // Constructor from an isce3::product::RadarGridProduct
-        inline ResampSlc(const isce3::product::RadarGridProduct &product, char frequency = 'A') :
-            isce3::image::ResampSlc(product, frequency) {}
+        inline ResampSlc(const isce3::product::RadarGridProduct &product, char frequency = 'A',
+                         std::complex<float> invalid_value = std::complex<float>(0.0, 0.0)) :
+            isce3::image::ResampSlc(product, frequency, invalid_value) {}
 
         // Constructor from an isce3::product::RadarGridProduct and reference product (flattening)
         inline ResampSlc(const isce3::product::RadarGridProduct & product,
                          const isce3::product::RadarGridProduct & refProduct,
-                         char frequency = 'A') :
-            isce3::image::ResampSlc(product, refProduct, frequency) {}
+                         char frequency = 'A',
+                         std::complex<float> invalid_value = std::complex<float>(0.0, 0.0)) :
+            isce3::image::ResampSlc(product, refProduct, frequency, invalid_value) {}
 
         /** Constructor from an isce3::product::RadarGridParameters (no flattening) */
         inline ResampSlc(const isce3::product::RadarGridParameters & rdr_grid,
-                         const isce3::core::LUT2d<double> & doppler) :
-            isce3::image::ResampSlc(rdr_grid, doppler) {}
+                         const isce3::core::LUT2d<double> & doppler,
+                         std::complex<float> invalid_value = std::complex<float>(0.0, 0.0)) :
+            isce3::image::ResampSlc(rdr_grid, doppler, invalid_value) {}
 
         /** Constructor from an isce3::product::RadarGridParameters and reference radar grid (flattening) */
         inline ResampSlc(const isce3::product::RadarGridParameters & rdr_grid,
                          const isce3::product::RadarGridParameters & ref_rdr_grid,
                          const isce3::core::LUT2d<double> & doppler,
-                         double wvl, double ref_wvl) :
-            isce3::image::ResampSlc(rdr_grid, ref_rdr_grid, doppler, wvl, ref_wvl) {}
+                         std::complex<float> invalid_value = std::complex<float>(0.0, 0.0)) :
+            isce3::image::ResampSlc(rdr_grid, ref_rdr_grid, doppler, invalid_value) {}
 
         // Constructor from individual components (no flattening)
         inline ResampSlc(const isce3::core::LUT2d<double> & doppler,
                          double startingRange, double rangePixelSpacing,
-                         double sensingStart, double prf, double wvl) :
+                         double sensingStart, double prf, double wvl,
+                         std::complex<float> invalid_value = std::complex<float>(0.0, 0.0)) :
             isce3::image::ResampSlc(doppler, startingRange, rangePixelSpacing, sensingStart,
-                                   prf, wvl) {}
+                                   prf, wvl, invalid_value) {}
 
         // Constructor from individual components (flattening)
         inline ResampSlc(const isce3::core::LUT2d<double> & doppler,
                          double startingRange, double rangePixelSpacing,
                          double sensingStart, double prf, double wvl,
                          double refStartingRange, double refRangePixelSpacing,
-                         double refWvl) :
+                         double refWvl,
+                         std::complex<float> invalid_value = std::complex<float>(0.0, 0.0)) :
             isce3::image::ResampSlc(doppler, startingRange, rangePixelSpacing, sensingStart,
-                                   prf, wvl, refStartingRange, refRangePixelSpacing, refWvl) {}
+                                   prf, wvl, refStartingRange, refRangePixelSpacing, refWvl,
+                                   invalid_value) {}
 
-        // All resamp need? to be redefined to ensure derived functions used
-        // Generic resamp entry point from externally created rasters
+        /* Generic resamp entry point from externally created rasters
+         *
+         * \param[in] inputSlc          raster of SLC to be resampled
+         * \param[in] outputSlc         raster of resampled SLC
+         * \param[in] rgOffsetRaster    raster of range shift to be applied
+         * \param[in] azOffsetRaster    raster of azimuth shift to be applied
+         * \param[in] inputBand         band of input raster to resample
+         * \param[in] flatten           flag to flatten resampled SLC
+         * \param[in] rowBuffer         number of rows excluded from top/bottom of azimuth
+         *                              raster while searching for min/max indices of
+         *                              resampled SLC
+         * \param[in] chipSize          size of chip used in sinc interpolation
+         */
         void resamp(isce3::io::Raster & inputSlc, isce3::io::Raster & outputSlc,
                     isce3::io::Raster & rgOffsetRaster, isce3::io::Raster & azOffsetRaster,
-                    int inputBand=1, bool flatten=false, bool isComplex=true, int rowBuffer=40,
+                    int inputBand=1, bool flatten=false, int rowBuffer=40,
                     int chipSize=isce3::core::SINC_ONE);
 
-        // Generic resamp entry point: use filenames to create rasters
+        /* Generic resamp entry point: use filenames to create rasters
+         * internally in function.
+         *
+         * \param[in] inputFilename     path of file containing SLC to be resampled
+         * \param[in] outputFilename    path of file containing resampled SLC
+         * \param[in] rgOffsetFilename  path of file containing range shift to be applied
+         * \param[in] azOffsetFilename  path of file containing azimuth shift to be applied
+         * \param[in] inputBand         band of input raster to resample
+         * \param[in] flatten           flag to flatten resampled SLC
+         * \param[in] rowBuffer         number of rows excluded from top/bottom of azimuth
+         *                              raster while searching for min/max indices of
+         *                              resampled SLC
+         * \param[in] chipSize          size of chip used in sinc interpolation
+         */
         void resamp(const std::string & inputFilename, const std::string & outputFilename,
                     const std::string & rgOffsetFilename, const std::string & azOffsetFilename,
-                    int inputBand=1, bool flatten=false, bool isComplex=true, int rowBuffer=40,
+                    int inputBand=1, bool flatten=false, int rowBuffer=40,
                     int chipSize=isce3::core::SINC_ONE);
 
 };