Integrate SNAPHU in InSAR workflow (#968)

Virginia Brancato · Liang Yu · GitHub Enterprise · commit 3f6270964131 · 2022-04-21T12:41:16.000-07:00
* Move correlation_thresh_increments to common ICU, Phass parameters

* Add snaphu parameters in InSAR schema

* Re-organize unwrap.py to include snaphu unwrapper

* Integrate SNAPHU in InSAR workflow

* Correct rg bandwidth, conn components to uint32, fix stats

* Add functionality to clip nan from correlation data

* add support for unsigned int rasters to CUDA geocode

* Fix computation azimuth resolution. Log ICU and Phass attributes

* Correct typo in clean cfg function

* Rename some functionalities for clarity

* Expose and include MCF initializer

* add raster data type and interpolation method check in CUDA geocode

* Remove default value comments for SNAPHU in InSAR schema

* Remove correlation and phase std model from schema and unwrap.py

* Fix identation doc strings

* Change Snaphu schema params name to match Python wrapper

* Use ternary if conditions and fix params in schema

* Resolve param name diff between Python interface and InSAR schema

* limit int type interp to nearest neighbor

* Move unwrap default values to InSAR default

* Fix pep8 notation mistakes

* remove erroneously introduced file

* fix incorrect init logic

* added int type and interp type check
whitespace cleanup

* Include SNAPHU in Rosamond insar test case

* Remove function cleaning NaNs

* Correct bug on band from where to read data

* Move unwrap-related comments from schema to InSAR default runconfig

* Correct get function for cfg info extraction

* Fix bug of connected component type

Co-authored-by: Liang Yu &lt;Liang.Yu@jpl.nasa.gov&gt;
diff --git a/cxx/isce3/cuda/core/InterpolatorHandle.cu b/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
diff --git a/cxx/isce3/cuda/core/gpuBicubicInterpolator.cu b/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,
diff --git a/cxx/isce3/cuda/core/gpuBilinearInterpolator.cu b/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,
diff --git a/cxx/isce3/cuda/core/gpuNearestNeighborInterpolator.cu b/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>;
diff --git a/cxx/isce3/cuda/core/gpuSpline2dInterpolator.cu b/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,
diff --git a/cxx/isce3/cuda/geocode/Geocode.cu b/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
diff --git a/cxx/isce3/cuda/geocode/Geocode.h b/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
diff --git a/cxx/isce3/geocode/GeocodeCov.cpp b/cxx/isce3/geocode/GeocodeCov.cpp
@@ -105,7 +105,7 @@ void getBlocksNumberAndLength(const int array_length, const int array_width,
         const long long max_block_size)
 {
 
-    const int max_block_length = max_block_size / 
+    const int max_block_length = max_block_size /
         (nbands * array_width * type_size);
 
     int _block_length = std::min(array_length, max_block_length);
@@ -290,6 +290,23 @@ void Geocode<T>::geocodeInterp(
     std::unique_ptr<isce3::core::ProjectionBase> proj(
             isce3::core::createProj(_epsgOut));
 
+    // make sure int type rasters only used with nearest neighbor
+    for (int band = 0; band < nbands; ++band)
+    {
+        const auto dtype = inputRaster.dtype(band + 1);
+        if ((dtype == GDT_Byte || dtype == GDT_UInt32)
+                && _data_interp_method != isce3::core::NEAREST_METHOD)
+        {
+            pyre::journal::error_t error(
+                    "isce3.geocode.GeocodeCov.geocodeInterp");
+            error << "int type of raster can only use nearest neighbor interp"
+                << pyre::journal::endl;
+            std::string err_str {
+                "int type of raster can only use nearest neighbor interp"};
+            throw isce3::except::InvalidArgument(ISCE_SRCINFO(), err_str);
+        }
+    }
+
     // create data interpolator
     std::unique_ptr<isce3::core::Interpolator<T_out>> interp {
             isce3::core::createInterpolator<T_out>(_data_interp_method)};
@@ -374,7 +391,7 @@ void Geocode<T>::geocodeInterp(
             if (std::isnan(rtc_geogrid_upsampling))
                 rtc_geogrid_upsampling = 1;
 
-            isce3::core::MemoryModeBlockY rtc_memory_mode = 
+            isce3::core::MemoryModeBlockY rtc_memory_mode =
                 isce3::core::MemoryModeBlockY::AutoBlocksY;
             int radar_grid_nlooks = 1;
             computeRtc(demRaster, *rtc_raster, radar_grid, _orbit, _doppler,
@@ -555,9 +572,9 @@ void Geocode<T>::geocodeInterp(
         azimuthFirstLine = std::max(azimuthFirstLine - interp_margin, 0);
         rangeFirstPixel = std::max(rangeFirstPixel - interp_margin, 0);
 
-        azimuthLastLine = std::min(azimuthLastLine + interp_margin, 
+        azimuthLastLine = std::min(azimuthLastLine + interp_margin,
                                    static_cast<int>(radar_grid.length() - 1));
-        rangeLastPixel = std::min(rangeLastPixel + interp_margin, 
+        rangeLastPixel = std::min(rangeLastPixel + interp_margin,
                                   static_cast<int>(radar_grid.width() - 1));
 
         if (azimuthFirstLine > azimuthLastLine ||
@@ -706,7 +723,7 @@ void Geocode<T>::geocodeInterp(
     auto elapsed_time_milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(
         std::chrono::high_resolution_clock::now() - start_time);
     float elapsed_time = ((float) elapsed_time_milliseconds.count()) / 1e3;
-    info << "elapsed time (GEO-IN) [s]: " << elapsed_time << pyre::journal::endl;    
+    info << "elapsed time (GEO-IN) [s]: " << elapsed_time << pyre::journal::endl;
 }
 
 template<class T>
@@ -1048,7 +1065,7 @@ void _processUpsampledBlock(isce3::core::Matrix<T_out>* mat, size_t block,
     In this case, the input type T is expected to be complex and the output
     type T_out can be real or complex.
     The conversion from complex (e.g. SLC) to real (SAR backscatter)
-    in the context of geocoding is considered as the geocoding of 
+    in the context of geocoding is considered as the geocoding of
     the covariance matrix (diagonal elements) is done by
     squaring the modulus of the complex input. The conversion between
     variables of same type is considered as regular geocoding and
@@ -1115,17 +1132,17 @@ void _getUpsampledBlock(
                     std::cout.flush();
                 }
                 auto ptr = rdrData[band]->data();
-                input_raster.getBlock(ptr + 
-                                      block * radar_block_length * size_x, 
+                input_raster.getBlock(ptr +
+                                      block * radar_block_length * size_x,
                                       xidx, block * radar_block_length + yidx, size_x,
                                       this_radar_block_length, band + 1);
             }
 
             if (radargrid_nblocks > 1) {
-                std::cout << "reading band " << band + 1 
+                std::cout << "reading band " << band + 1
                           << " progress: 100%" << std::endl;
             }
-        } 
+        }
         else if (!flag_upsample_radar_grid && std::is_same<T, T_out>::value) {
             /*
             Enter here if:
@@ -1139,8 +1156,8 @@ void _getUpsampledBlock(
             {
             input_raster.getBlock(rdrData[band]->data(), xidx, yidx, size_x,
                                   size_y, band + 1);
-            }    
-        } 
+            }
+        }
         else if (!flag_upsample_radar_grid) {
             /*
             Enter here if:
@@ -1901,7 +1918,7 @@ void Geocode<T>::geocodeAreaProj(
     auto elapsed_time_milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(
         std::chrono::high_resolution_clock::now() - start_time);
     float elapsed_time = ((float) elapsed_time_milliseconds.count()) / 1e3;
-    info << "elapsed time (GEO-AP) [s]: " << elapsed_time << pyre::journal::endl;     
+    info << "elapsed time (GEO-AP) [s]: " << elapsed_time << pyre::journal::endl;
 }
 
 template<class T>
@@ -2636,7 +2653,7 @@ void Geocode<T>::_runBlock(
             }
 
             if (flag_self_intersecting_area_element) {
-                /* 
+                /*
                 If self-intersecting, divide area element (geogrid pixel) into
                 two triangles and integrate them separately.
                 */
@@ -2659,7 +2676,7 @@ void Geocode<T>::_runBlock(
                         x10_cut, size_y, size_x, w_arr_2, w_total_2, plane_orientation);
                 isce3::geometry::areaProjIntegrateSegment(y10_cut, y00_cut, x10_cut,
                         x00_cut, size_y, size_x, w_arr_2, w_total_2, plane_orientation);
-                
+
                 w_total = 0;
                 /*
                 The new weight array `w_arr` is the sum of the absolute values of both
@@ -2882,7 +2899,7 @@ void Geocode<T>::_runBlock(
         for (int band = 0; band < nbands_off_diag_terms; ++band) {
             for (int i = 0; i < this_block_size_y; ++i) {
                 for (int j = 0; j < this_block_size_x; ++j) {
-                    /* 
+                    /*
                     Since std::numeric_limits<T_out>::quiet_NaN() with
                     complex T_out is (or may be) undefined, we take the "real type"
                     if T_out (i.e. float or double) to create the NaN value and
diff --git a/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp b/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp
@@ -101,6 +101,10 @@ void addbinding(pybind11::class_<Geocode>& pyGeocode)
                         self.geocodeRasterBlock<unsigned char>(
                                 output_raster, input_raster);
                         break;}
+                    case GDT_UInt32:  {
+                        self.geocodeRasterBlock<unsigned int>(
+                                output_raster, input_raster);
+                        break;}
                     default: {
                         throw isce3::except::RuntimeError(ISCE_SRCINFO(),
                                 "unsupported datatype");
diff --git a/python/packages/nisar/workflows/h5_prep.py b/python/packages/nisar/workflows/h5_prep.py
@@ -746,7 +746,7 @@ def prep_ds_insar(pcfg, dst, dst_h5):
                                         units="radians", grids=grids_val,
                                         long_name='unwrapped phase')
                       descr = f"Connected components for {pol} layer"
-                      _create_datasets(dst_h5[pol_path], igram_shape, np.uint8,
+                      _create_datasets(dst_h5[pol_path], igram_shape, np.uint32,
                                       'connectedComponents', descr=descr, units=" ",
                                       grids=grids_val,
                                       long_name='connected components')
diff --git a/python/packages/nisar/workflows/unwrap.py b/python/packages/nisar/workflows/unwrap.py
diff --git a/share/nisar/defaults/insar.yaml b/share/nisar/defaults/insar.yaml
diff --git a/share/nisar/schemas/insar.yaml b/share/nisar/schemas/insar.yaml
diff --git a/tests/data/bandpass_test.yaml b/tests/data/bandpass_test.yaml
diff --git a/tools/imagesets/runconfigs/end2end_pair_ALPSRP110160680_ALPSRP116870680_Rosamond_insar.yaml b/tools/imagesets/runconfigs/end2end_pair_ALPSRP110160680_ALPSRP116870680_Rosamond_insar.yaml
