Merge pull request #11 from isce-3/develop

rebase develop

Author: jungkyoJung

cxx/isce3/core/Attitude.h

@@ -59,6 +59,11 @@ class Attitude {
     /** Time of last measurement relative to reference epoch (s) */
     double endTime() const { return _time[size() - 1]; }
 
+    /** Check if time falls in the valid interpolation domain. */
+    bool contains(double time) const {
+        return (startTime() <= time) && (time <= endTime());
+    }
+
     /** UTC time of first measurement */
     DateTime startDateTime() const
     {

cxx/isce3/core/Orbit.h

@@ -76,7 +76,12 @@ class Orbit {
     /** Reference epoch (UTC) */
     const DateTime & referenceEpoch() const { return _reference_epoch; }
 
-    /** Set reference epoch (UTC) */
+    /** Set reference epoch (UTC)
+     *
+     * Also updates relative time tags so that
+     *      referenceEpoch() + TimeDelta(time()[i])
+     * results in the same absolute time tags before and after this call.
+     */
     void referenceEpoch(const DateTime &);
 
     /** Interpolation method */
@@ -103,6 +108,11 @@ class Orbit {
     /** UTC time of last state vector */
     DateTime endDateTime() const { return _reference_epoch + TimeDelta(endTime()); }
 
+    /** Check if time falls in the valid interpolation domain. */
+    bool contains(double time) const {
+        return (startTime() <= time) && (time <= endTime());
+    }
+
     /** Time interval between state vectors (s) */
     double spacing() const { return _time.spacing(); }
 

cxx/isce3/cuda/core/InterpolatorHandle.cu

@@ -16,28 +16,23 @@ namespace isce3::cuda::core {
 
 template<class T>
 __global__ void init_interp(DeviceInterp<T>** interp,
-        isce3::core::dataInterpMethod interp_method, bool* unsupported_interp)
+        isce3::core::dataInterpMethod interp_method)
 {
     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;
-            case isce3::core::NEAREST_METHOD:
-                (*interp) = new isce3::cuda::core::gpuNearestNeighborInterpolator<T>();
-                break;
-            default:
-                *unsupported_interp = true;
-                break;
+        if (interp_method == isce3::core::BILINEAR_METHOD)
+            (*interp) = new isce3::cuda::core::gpuBilinearInterpolator<T>();
+
+        if (interp_method == isce3::core::BICUBIC_METHOD)
+            (*interp) = new isce3::cuda::core::gpuBicubicInterpolator<T>();
+
+        if (interp_method == isce3::core::BIQUINTIC_METHOD) {
+            size_t order = 6;
+            (*interp) = new isce3::cuda::core::gpuSpline2dInterpolator<T>(order);
         }
+
+        if (interp_method == isce3::core::NEAREST_METHOD)
+            (*interp) = new isce3::cuda::core::gpuNearestNeighborInterpolator<T>();
     }
 }
 
@@ -55,21 +50,23 @@ InterpolatorHandle<T>::InterpolatorHandle(
 {
     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) {
+    if (interp_method != isce3::core::BILINEAR_METHOD
+            && interp_method != isce3::core::BICUBIC_METHOD
+            && interp_method != isce3::core::BIQUINTIC_METHOD
+            && interp_method != isce3::core::NEAREST_METHOD)
+    {
         pyre::journal::error_t error(
-                "isce.cuda.core.InterpolatorHandle.InterpolatorHandle");
+                "isce3.cuda.core.InterpolatorHandle.InterpolatorHandle");
         error << "Unsupported interpolator method provided."
               << pyre::journal::endl;
         throw isce3::except::InvalidArgument(
                 ISCE_SRCINFO(), "Unsupported interpolator method provided.");
     }
+
+    thrust::device_vector<bool> d_unsupported_interp(1, false);
+    init_interp<<<1, 1>>>(_interp, interp_method);
+    checkCudaErrors(cudaPeekAtLastError());
+    checkCudaErrors(cudaDeviceSynchronize());
 }
 
 template<class T>
@@ -86,4 +83,5 @@ template class InterpolatorHandle<thrust::complex<float>>;
 template class InterpolatorHandle<double>;
 template class InterpolatorHandle<thrust::complex<double>>;
 template class InterpolatorHandle<unsigned char>;
+template class InterpolatorHandle<unsigned int>;
 } // namespace isce3::cuda::core

cxx/isce3/cuda/core/gpuBicubicInterpolator.cu

@@ -179,6 +179,7 @@ template class gpuBicubicInterpolator<thrust::complex<double>>;
 template class gpuBicubicInterpolator<float>;
 template class gpuBicubicInterpolator<thrust::complex<float>>;
 template class gpuBicubicInterpolator<unsigned char>;
+template class gpuBicubicInterpolator<unsigned int>;
 
 template __global__ void gpuInterpolator_g<double>(
         gpuBicubicInterpolator<double> interp, double* x, double* y,

cxx/isce3/cuda/core/gpuBilinearInterpolator.cu

@@ -115,6 +115,7 @@ template class gpuBilinearInterpolator<thrust::complex<double>>;
 template class gpuBilinearInterpolator<float>;
 template class gpuBilinearInterpolator<thrust::complex<float>>;
 template class gpuBilinearInterpolator<unsigned char>;
+template class gpuBilinearInterpolator<unsigned int>;
 
 template __global__ void gpuInterpolator_g<double>(
         gpuBilinearInterpolator<double> interp, double* x, double* y,

cxx/isce3/cuda/core/gpuNearestNeighborInterpolator.cu

@@ -19,3 +19,4 @@ template class gpuNearestNeighborInterpolator<thrust::complex<double>>;
 template class gpuNearestNeighborInterpolator<float>;
 template class gpuNearestNeighborInterpolator<thrust::complex<float>>;
 template class gpuNearestNeighborInterpolator<unsigned char>;
+template class gpuNearestNeighborInterpolator<unsigned int>;

cxx/isce3/cuda/core/gpuSpline2dInterpolator.cu

@@ -154,6 +154,7 @@ template class gpuSpline2dInterpolator<thrust::complex<double>>;
 template class gpuSpline2dInterpolator<float>;
 template class gpuSpline2dInterpolator<thrust::complex<float>>;
 template class gpuSpline2dInterpolator<unsigned char>;
+template class gpuSpline2dInterpolator<unsigned int>;
 
 template __global__ void gpuInterpolator_g<double>(
         gpuSpline2dInterpolator<double> interp, double* x, double* y,

cxx/isce3/cuda/geocode/Geocode.cu

@@ -201,8 +201,10 @@ Geocode::Geocode(const isce3::product::GeoGridParameters & geogrid,
     _interp_double_handle(data_interp_method),
     _interp_cdouble_handle(data_interp_method),
     _interp_unsigned_char_handle(data_interp_method),
+    _interp_unsigned_int_handle(data_interp_method),
     _proj_handle(geogrid.epsg()),
-    _dem_interp_method(dem_interp_method)
+    _dem_interp_method(dem_interp_method),
+    _data_interp_method(data_interp_method)
 {
     // init light weight radar grid
     _radar_grid.sensing_start = _rdr_geom.radarGrid().sensingStart();
@@ -231,10 +233,12 @@ Geocode::Geocode(const isce3::product::GeoGridParameters & geogrid,
         _invalid_float = std::numeric_limits<float>::quiet_NaN();
         _invalid_double = std::numeric_limits<double>::quiet_NaN();
         _invalid_unsigned_char = 255;
+        _invalid_unsigned_int = 4294967295;
     } else {
         _invalid_float = invalid_value;
         _invalid_double = static_cast<double>(invalid_value);
         _invalid_unsigned_char = static_cast<unsigned char>(invalid_value);
+        _invalid_unsigned_int = static_cast<unsigned int>(invalid_value);
     }
 }
 
@@ -337,9 +341,21 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number)
     }
 }
 
+void Geocode::rasterDtypeInterpCheck(const int dtype) const
+{
+    if ((dtype == GDT_Byte || dtype == GDT_UInt32) &&
+        _data_interp_method != isce3::core::NEAREST_METHOD)
+    {
+        std::string err_str {"int type of raster can only use nearest neighbor interp"};
+        throw isce3::except::InvalidArgument(ISCE_SRCINFO(), err_str);
+    }
+}
+
 template<class T>
 void Geocode::geocodeRasterBlock(Raster& output_raster, Raster& input_raster)
 {
+    rasterDtypeInterpCheck(input_raster.dtype());
+
     // determine number of elements in output vector
     const auto n_elem_out = _geo_block_length * _geogrid.width();
 
@@ -362,6 +378,9 @@ void Geocode::geocodeRasterBlock(Raster& output_raster, Raster& input_raster)
     } else if constexpr (std::is_same_v<T, unsigned char>) {
         interp = _interp_unsigned_char_handle.getInterp();
         invalid_value = _invalid_unsigned_char;
+    } else if constexpr (std::is_same_v<T, unsigned int>) {
+        interp = _interp_unsigned_int_handle.getInterp();
+        invalid_value = _invalid_unsigned_int;
     }
 
     // 0 width indicates current block is out of bounds
@@ -427,6 +446,10 @@ void Geocode::geocodeRasters(
 
     const auto n_raster_pairs = output_rasters.size();
 
+    // check if raster types consistent with data interp method
+    for (size_t i_raster = 0; i_raster < n_raster_pairs; ++i_raster)
+        rasterDtypeInterpCheck(input_rasters[i_raster].get().dtype());
+
     // iterate over blocks
     for (size_t i_block = 0; i_block < _n_blocks; ++i_block) {
 
@@ -457,6 +480,10 @@ void Geocode::geocodeRasters(
                     geocodeRasterBlock<unsigned char>(
                             output_rasters[i_raster], input_rasters[i_raster]);
                     break;}
+                case GDT_UInt32:  {
+                    geocodeRasterBlock<unsigned int>(
+                            output_rasters[i_raster], input_rasters[i_raster]);
+                    break;}
                 default: {
                     throw isce3::except::RuntimeError(ISCE_SRCINFO(),
                             "unsupported datatype");
@@ -474,4 +501,5 @@ EXPLICIT_INSTATIATION(thrust::complex<float>);
 EXPLICIT_INSTATIATION(double);
 EXPLICIT_INSTATIATION(thrust::complex<double>);
 EXPLICIT_INSTATIATION(unsigned char);
+EXPLICIT_INSTATIATION(unsigned int);
 } // namespace isce3::cuda::geocode

cxx/isce3/cuda/geocode/Geocode.h

@@ -178,12 +178,18 @@ class Geocode {
             _interp_cdouble_handle;
     isce3::cuda::core::InterpolatorHandle<unsigned char>
             _interp_unsigned_char_handle;
+    isce3::cuda::core::InterpolatorHandle<unsigned int>
+            _interp_unsigned_int_handle;
 
     // value applied to invalid geogrid pixels
     float _invalid_float;
     double _invalid_double;
     unsigned char _invalid_unsigned_char;
+    unsigned int _invalid_unsigned_int;
 
     isce3::core::dataInterpMethod _dem_interp_method;
+    isce3::core::dataInterpMethod _data_interp_method;
+
+    void rasterDtypeInterpCheck(const int dtype) const;
 };
 } // namespace isce3::cuda::geocode

cxx/isce3/cuda/geometry/Geo2rdr.cpp

@@ -102,12 +102,9 @@ geo2rdr(isce3::io::Raster & topoRaster,
     // Interpolate orbit to middle of the scene as a test
     _checkOrbitInterpolation(tmid);
 
-    // Compute number of lines per block
-    computeLinesPerBlock();
-
     // Compute number of DEM blocks needed to process image
-    size_t nBlocks = demLength / _linesPerBlock;
-    if ((demLength % _linesPerBlock) != 0)
+    size_t nBlocks = demLength / linesPerBlock();
+    if ((demLength % linesPerBlock()) != 0)
         nBlocks += 1;
 
     // Cache near, mid, far ranges for diagnostics on Doppler
@@ -121,11 +118,11 @@ geo2rdr(isce3::io::Raster & topoRaster,
 
         // Get block extents
         size_t lineStart, blockLength;
-        lineStart = block * _linesPerBlock;
+        lineStart = block * linesPerBlock();
         if (block == (nBlocks - 1)) {
             blockLength = demLength - lineStart;
         } else {
-            blockLength = _linesPerBlock;
+            blockLength = linesPerBlock();
         }
         size_t blockSize = blockLength * demWidth;
 
@@ -194,33 +191,3 @@ _checkOrbitInterpolation(double aztime) {
     isce3::core::cartesian_t satxyz, satvel;
     this->orbit().interpolate(&satxyz, &satvel, aztime);
 }
-
-// Compute number of lines per block dynamically from GPU memory
-void isce3::cuda::geometry::Geo2rdr::
-computeLinesPerBlock() {
-
-    // Compute GPU memory
-    const size_t nGPUBytes = getDeviceMem();
-
-    // 2 GB buffer for safeguard for large enough devices (> 6 GB)
-    size_t gpuBuffer;
-    if (nGPUBytes > 6e9) {
-        gpuBuffer = 2e9;
-    } else {
-        // Else use 500 MB buffer
-        gpuBuffer = 500e6;
-    }
-
-    // Compute pixels per Block (3 double input layers and 2 float output layers)
-    size_t pixelsPerBlock = (nGPUBytes - gpuBuffer) /
-                            (3 * sizeof(double) + 2 * sizeof(float));
-    // Round down to nearest 10 million pixels
-    pixelsPerBlock = (pixelsPerBlock / 10000000) * 10000000;
-
-    // Compute number of lines per block
-    _linesPerBlock = pixelsPerBlock / this->radarGridParameters().width();
-    // Round down to nearest 500 lines
-    _linesPerBlock = (_linesPerBlock / 500) * 500;
-}
-
-// end of file