GSLC carriers (#912)

* add az and range phase carrier params to geocode SLC

* fix pybind

* added missing docstring info

* removed extra semicolon

* Remove az and rg carriers before interpolation

* reuse carrier_phase computation

* GSLC carrier removal like resample SLC

* phase to chipPhase

* add offset in range and slight shuffle/reorg

* correct carrier restoration indices

* fix misordered eval inputs

* parallelize carrier removal

* simplify indexing

* GSLC frequency removal/reintroduction made modular

* moved doppler calculation back to interpolate
added files for carrier remove/restore

* updated documentation + comments and some cleanup

* geocodeSlc templated

* move ramping/interpolate functions inside geocodeSlc.cpp

* correct frequency and phrase usage

* self describing template name

* replace std::valarray with isce3::core::matrix

* added more comments for clarity
renamed variables for clarity
cleaned up commented out code

* carriers eval with az/range instead of indices

* fix sign bug
rename carrierRampAndFlatten to carrierPhaseRampAndFlatten for consistency
documentation updates and fixes

* fix documentation and typo

* refactor double loops to single for speed up

* renames for clarity and docstring fix

* add pybind docstring

* apply doppler frac to chip instead of interped value

* fill invalid pixels to NaN

Co-authored-by: vbrancat <[email protected]>

Author: 2 people

cxx/isce3/Headers.cmake

@@ -83,7 +83,6 @@ focus/Presum.icc
 focus/RangeComp.h
 geocode/baseband.h
 geocode/geocodeSlc.h
-geocode/interpolate.h
 geocode/loadDem.h
 geometry/DEMInterpolator.h
 geometry/forward.h

cxx/isce3/Sources.cmake

@@ -43,7 +43,6 @@ focus/Presum.cpp
 focus/RangeComp.cpp
 geocode/baseband.cpp
 geocode/geocodeSlc.cpp
-geocode/interpolate.cpp
 geocode/loadDem.cpp
 geometry/DEMInterpolator.cpp
 geometry/Geo2rdr.cpp

cxx/isce3/geocode/geocodeSlc.cpp

@@ -1,6 +1,7 @@
 #pragma once
 #include <cstddef>
 #include <isce3/core/forward.h>
+#include <isce3/core/Poly2d.h>
 #include <isce3/io/forward.h>
 #include <isce3/product/forward.h>
 
@@ -14,15 +15,19 @@ namespace isce3 { namespace geocode {
  * \param[in]  radarGrid     radar grid parameters
  * \param[in]  geoGrid       geo grid parameters
  * \param[in]  orbit             orbit
- * \param[in]  nativeDoppler     2D LUT Doppler of the SLC image 
- * \param[in]  imageGridDoppler  2D LUT Doppler of the image grid 
- * \param[in]  ellipsoid         ellipsoid object 
- * \param[in]  thresholdGeo2rdr  threshold for geo2rdr computations 
+ * \param[in]  nativeDoppler     2D LUT Doppler of the SLC image
+ * \param[in]  imageGridDoppler  2D LUT Doppler of the image grid
+ * \param[in]  ellipsoid         ellipsoid object
+ * \param[in]  thresholdGeo2rdr  threshold for geo2rdr computations
  * \param[in]  numiterGeo2rdr    maximum number of iterations for Geo2rdr convergence
  * \param[in]  linesPerBlock     number of lines in each block
  * \param[in]  demBlockMargin    margin of a DEM block in degrees
  * \param[in]  flatten           flag to flatten the geocoded SLC
+ * \param[in]  azCarrier        azimuth carrier phase of the SLC data, in radians, as a function of azimuth and range
+ * \param[in]  rgCarrier        range carrier phase of the SLC data, in radians, as a function of azimuth and range
+ * \tparam[in]  AzRgFunc  2-D real-valued function of azimuth and range
  */
+template<typename AzRgFunc = isce3::core::Poly2d>
 void geocodeSlc(isce3::io::Raster& outputRaster, isce3::io::Raster& inputRaster,
                 isce3::io::Raster& demRaster,
                 const isce3::product::RadarGridParameters& radarGrid,
@@ -33,6 +38,8 @@ void geocodeSlc(isce3::io::Raster& outputRaster, isce3::io::Raster& inputRaster,
                 const isce3::core::Ellipsoid& ellipsoid,
                 const double& thresholdGeo2rdr, const int& numiterGeo2rdr,
                 const size_t& linesPerBlock, const double& demBlockMargin,
-                const bool flatten = true);
+                const bool flatten = true,
+                const AzRgFunc& azCarrier = AzRgFunc(),
+                const AzRgFunc& rgCarrier = AzRgFunc());
 
 }} // namespace isce3::geocode

cxx/isce3/geocode/geocodeSlc.h

@@ -3,6 +3,7 @@
 #include <isce3/core/Ellipsoid.h>
 #include <isce3/core/LUT2d.h>
 #include <isce3/core/Orbit.h>
+#include <isce3/core/Poly2d.h>
 #include <isce3/io/Raster.h>
 #include <isce3/product/RadarGridParameters.h>
 #include <isce3/product/GeoGridParameters.h>
@@ -11,9 +12,10 @@
 
 namespace py = pybind11;
 
+template<typename AzRgFunc = isce3::core::Poly2d>
 void addbinding_geocodeslc(py::module & m)
 {
-    m.def("geocode_slc", &isce3::geocode::geocodeSlc,
+    m.def("geocode_slc", &isce3::geocode::geocodeSlc<AzRgFunc>,
         py::arg("output_raster"),
         py::arg("input_raster"),
         py::arg("dem_raster"),
@@ -27,5 +29,46 @@ void addbinding_geocodeslc(py::module & m)
         py::arg("numiter_geo2rdr") = 25,
         py::arg("lines_per_block") = 1000,
         py::arg("dem_block_margin") = 0.1,
-        py::arg("flatten") = true);
+        py::arg("flatten") = true,
+        py::arg("azimuth_carrier") = AzRgFunc(),
+        py::arg("range_carrier") = AzRgFunc(),
+        R"(
+        Geocode a SLC
+
+        Parameters
+        ----------
+        outputRaster: Raster
+            Output raster containing geocoded SLC
+        inputRaster: Raster
+            Input raster of the SLC in radar coordinates
+        demRaster: Raster
+            Raster of the DEM
+        radargrid: RadarGridParameters
+            Radar grid parameters of input SLC raster
+        geogrid: GeoGridParameters
+            Geo grid parameters of output raster
+        native_doppler: LUT2d
+            2D LUT doppler of the SLC image
+        image_grid_doppler: LUT2d
+            2d LUT doppler of the image grid
+        ellipsoid: Ellipsoid
+            Ellipsoid object
+        threshold_geo2rdr: float
+            Threshold for geo2rdr computations
+        numiter_geo2rdr: int
+            Maximum number of iterations for geo2rdr convergence
+        lines_per_block: int
+            Number of lines per block
+        dem_block_margin: float
+            Margin of a DEM block in degrees
+        flatten: bool
+            Flag to flatten the geocoded SLC
+        azimuth_carrier: [LUT2d, Poly2d]
+            Azimuth carrier phase of the SLC data, in radians, as a function of azimuth and range
+        range_carrier: [LUT2d, Poly2d]
+            Range carrier phase of the SLC data, in radians, as a function of azimuth and range
+        )");
 }
+
+template void addbinding_geocodeslc<isce3::core::LUT2d<double>>(py::module & m);
+template void addbinding_geocodeslc<isce3::core::Poly2d>(py::module & m);

cxx/isce3/geocode/interpolate.cpp

@@ -2,4 +2,5 @@
 
 #include <pybind11/pybind11.h>
 
+template<typename T>
 void addbinding_geocodeslc(pybind11::module&);

cxx/isce3/geocode/interpolate.h

@@ -10,7 +10,8 @@ void addsubmodule_geocode(py::module & m)
 {
     py::module geocode = m.def_submodule("geocode");
 
-    addbinding_geocodeslc(geocode);
+    addbinding_geocodeslc<isce3::core::LUT2d<double>>(geocode);
+    addbinding_geocodeslc<isce3::core::Poly2d>(geocode);
 
     // forward declare bound classes
     py::class_<isce3::geocode::Geocode<float>>

python/extensions/pybind_isce3/geocode/GeocodeSlc.cpp

@@ -11,6 +11,7 @@
 #include <isce3/core/LUT2d.h>
 #include <isce3/core/Metadata.h>
 #include <isce3/core/Orbit.h>
+#include <isce3/core/Poly2d.h>
 #include <isce3/geocode/geocodeSlc.h>
 #include <isce3/geocode/GeocodeCov.h>
 #include <isce3/geometry/Topo.h>
@@ -185,7 +186,7 @@ TEST(GeocodeTest, TestGeocodeCov) {
                     offset_az_raster, offset_rg_raster, geocode_memory_mode_1,
                     min_block_size, max_block_size);
         }
-    } 
+    }
 
     // Test generation of full-covariance elements and block processing
 
@@ -265,7 +266,7 @@ TEST(GeocodeTest, CheckGeocodeCovFullCovResults) {
     double err_y = 0.0;
     double err_x_conj_y = 0.0;
     double max_err_x = 0.0;
-    double max_err_y = 0.0; 
+    double max_err_y = 0.0;
     double max_err_x_conj_y = 0.0;
 
     std::valarray<double> geocoded_diag_array_x(length * width);
@@ -352,17 +353,17 @@ TEST(GeocodeTest, CheckGeocodeCovFullCovResults) {
             nvalid_x_conj_y++;
             // }
 
-            max_err_x_conj_y = std::max(max_err_x_conj_y, 
+            max_err_x_conj_y = std::max(max_err_x_conj_y,
                                         std::abs(err_x_conj_y));
-    
+
         }
     }
 
     stats.mean = total / static_cast<double>(stats.n_valid);
     double mean_real_valued = total_real_valued / stats.n_valid;
 
     double sample_stddev_factor = (
-        static_cast<double>(stats.n_valid)) / 
+        static_cast<double>(stats.n_valid)) /
         (static_cast<double>(stats.n_valid - 1));
 
     stats.sample_stddev = std::sqrt(sample_stddev_factor *
@@ -488,7 +489,7 @@ TEST(GeocodeTest, CheckGeocodeCovResults) {
         stats_y.mean = total_y / stats_y.n_valid;
 
         double sample_stddev_factor = (
-            static_cast<double>(stats_x.n_valid)) / 
+            static_cast<double>(stats_x.n_valid)) /
             (static_cast<double>(stats_x.n_valid - 1));
 
         stats_x.sample_stddev = std::sqrt(sample_stddev_factor *
@@ -510,7 +511,7 @@ TEST(GeocodeTest, CheckGeocodeCovResults) {
         std::cout << "  maxErrY: " << maxErrY << std::endl;
         std::cout << "  dx: " << dx << std::endl;
         std::cout << "  dy: " << dy << std::endl;
-    
+
         ASSERT_GE(stats_x.n_valid, 800);
         ASSERT_GE(stats_y.n_valid, 800);
 
@@ -603,20 +604,20 @@ TEST(GeocodeTest, TestGeocodeSlc)
 
     bool flatten = false;
 
-    isce3::geocode::geocodeSlc(geocodedSlc, inputSlc, demRaster, radarGrid,
-                              geoGrid, orbit, nativeDoppler, imageGridDoppler,
-                              ellipsoid, thresholdGeo2rdr, numiterGeo2rdr,
-                              linesPerBlock, demBlockMargin, flatten);
+    isce3::geocode::geocodeSlc<isce3::core::Poly2d> (
+            geocodedSlc, inputSlc, demRaster, radarGrid, geoGrid, orbit,
+            nativeDoppler, imageGridDoppler, ellipsoid, thresholdGeo2rdr,
+            numiterGeo2rdr, linesPerBlock, demBlockMargin, flatten);
 
     isce3::io::Raster inputSlcY("yslc_rdr.bin", GA_ReadOnly);
 
     isce3::io::Raster geocodedSlcY("yslc_geo.bin", geoGridWidth, geoGridLength, 1,
                                   GDT_CFloat32, "ENVI");
 
-    isce3::geocode::geocodeSlc(geocodedSlcY, inputSlcY, demRaster, radarGrid,
-                              geoGrid, orbit, nativeDoppler, imageGridDoppler,
-                              ellipsoid, thresholdGeo2rdr, numiterGeo2rdr,
-                              linesPerBlock, demBlockMargin, flatten);
+    isce3::geocode::geocodeSlc<isce3::core::LUT2d<double>> (
+            geocodedSlcY, inputSlcY, demRaster, radarGrid, geoGrid, orbit,
+            nativeDoppler, imageGridDoppler, ellipsoid, thresholdGeo2rdr,
+            numiterGeo2rdr, linesPerBlock, demBlockMargin, flatten);
 
     double* _geoTrans = new double[6];
     _geoTrans[0] = geoGridStartX;
@@ -726,7 +727,7 @@ void createZeroDem()
 
 
 template<class T>
-void checkStatsReal(isce3::math::Stats<T> computed_stats, 
+void checkStatsReal(isce3::math::Stats<T> computed_stats,
                     isce3::io::Raster &raster) {
 
         int band = 0, approx_ok = 0;
@@ -798,7 +799,7 @@ void checkStatsComplex(isce3::math::Stats<T> computed_stats,
     ASSERT_LT(std::abs(isce3_stats.max - computed_stats.max), 1.0e-15);
 
     if (!isnan(computed_stats.sample_stddev)) {
-        ASSERT_LT(std::abs(isce3_stats.sample_stddev - 
+        ASSERT_LT(std::abs(isce3_stats.sample_stddev -
                            computed_stats.sample_stddev), 1.0e-8);
     }