Add polarimetric symmetrization to the GCOV workflow (#752)

* add polsar namespace and symmetrizeCrossPolChannels() function

* add python bindings

* updating compiling and linking files

* add c++ unitest

* fix polsar import

* make complex multiply more general

* update GCOV SAS workflow

* update GCOV YAML defaults

* update GCOV YAML schema

* add python unitest

* substitute unitest extension .rdr with .bin

* add python test to CMakeLists.txt

* update calculation of nblocks and add docstrings

* update calculation of nblocks and add docstrings (2)

* update directory for unitest files

* merge develop into update_gcov_projection (2)

* Update cxx/isce3/polsar/symmetrize.h

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update tests/cxx/isce3/polsar/symmetrize.cpp

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update cxx/isce3/polsar/symmetrize.cpp

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update python/extensions/pybind_isce3/polsar/symmetrize.cpp

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update python/packages/pybind_nisar/workflows/gcov.py

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update tests/python/extensions/pybind/polsar/symmetrize.py

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update tests/python/extensions/pybind/polsar/symmetrize.py

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update tests/cxx/isce3/polsar/symmetrize.cpp

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* Update tests/cxx/isce3/polsar/symmetrize.cpp

Co-Authored-By: Geoffrey M Gunter <[email protected]>

* make non-public local functions static

* remove unnecessary white spaces

* improve type conversions

* add band arguments for cross-pol channels; improve docstrings

* improve comments in gcov.py

* improve type aliases in symmetrize unitest

* move memory block mode to cxx/isce3/core/Constants.h

* improve comments to memory block mode in the Y-direction

* Improve enum memoryModeBlockY docstrings

* pass std::function by const reference

* substitute enum with enum class

* update pybind11 enum class memory_mode_block_y docstrings

* Update cxx/isce3/geometry/RTC.cpp

Co-authored-by: Geoffrey M Gunter <[email protected]>

* improve docstrings, add output_raster_band

* add output_raster_band (2), improve dtype validation

* improve comments in runconfig YAML default and schema

* Update cxx/isce3/polsar/symmetrize.cpp

Co-authored-by: Geoffrey M Gunter <[email protected]>

* rename memoryModeBlockY to MemoryModeBlockY, improve enum class docstrings

* fix typo in python/extensions/pybind_isce3/polsar/symmetrize.cpp

* Update cxx/isce3/core/Constants.h

Co-authored-by: Geoffrey M Gunter <[email protected]>

* add band validation

* remove unused line

* improve symmetrizeCrossPolChannels() docstrings to explain kAutoBlocksY

* Update share/nisar/defaults/gcov.yaml

Co-authored-by: Geoffrey M Gunter <[email protected]>

* remove k prefix from enum class MemoryModeBlockY

* run "git-clang-format develop"

Co-authored-by: Gustavo H Shiroma <[email protected]>
Co-authored-by: Geoffrey M Gunter <[email protected]>

Author: 3 people

cxx/isce3/Headers.cmake

@@ -135,6 +135,7 @@ math/Bessel.h
 math/ComplexMultiply.h
 math/Sinc.h
 math/Sinc.icc
+polsar/symmetrize.h
 product/forward.h
 product/GeoGridParameters.h
 product/Metadata.h

cxx/isce3/Sources.cmake

@@ -71,6 +71,7 @@ matchtemplate/ampcor/correlators/tgtStats.cpp
 matchtemplate/ampcor/dom/Raster.cc
 matchtemplate/ampcor/dom/SLC.cc
 math/Bessel.cpp
+polsar/symmetrize.cpp
 product/GeoGridParameters.cpp
 product/Product.cpp
 product/RadarGridParameters.cpp

cxx/isce3/core/Constants.h

@@ -34,6 +34,16 @@ const int SINC_LEN = 8;
 const int SINC_ONE = 9;
 const int SINC_SUB = 8192;
 
+/** Enumeration type to indicate memory management for processes
+ * that require block processing in the Y direction
+ */
+enum class MemoryModeBlockY {
+    AutoBlocksY,     /**< auto mode (default value is defined by the
+                         module that is being executed) */
+    SingleBlockY,    /**< use a single block (disable block mode) */
+    MultipleBlocksY, /**< use multiple blocks (enable block mode) */
+};
+
 /** Convert string to dataInterpMethod */
 dataInterpMethod parseDataInterpMethod(const std::string & method);
 

cxx/isce3/geocode/GeocodeCov.cpp

@@ -1618,7 +1618,7 @@ void Geocode<T>::geocodeAreaProj(
 
     if (!std::isnan(min_nlooks))
         info << "nlooks min: " << min_nlooks << pyre::journal::newline;
-  
+
     // create projection based on epsg code
     std::unique_ptr<isce3::core::ProjectionBase> proj(
             isce3::core::createProj(_epsgOut));
@@ -1632,13 +1632,13 @@ void Geocode<T>::geocodeAreaProj(
             geogrid_upsampling, flag_upsample_radar_grid, dem_interp_method,
             &offset_y, &offset_x, &grid_size_y, &grid_size_x);
 
-    isce3::product::RadarGridParameters radar_grid_cropped =  
-        radar_grid.offsetAndResize(offset_y, offset_x, grid_size_y,
-                                   grid_size_x);
+    isce3::product::RadarGridParameters radar_grid_cropped =
+            radar_grid.offsetAndResize(
+                    offset_y, offset_x, grid_size_y, grid_size_x);
 
     bool is_radar_grid_single_block =
             (geocode_memory_mode !=
-             geocodeMemoryMode::BLOCKS_GEOGRID_AND_RADARGRID);
+                    geocodeMemoryMode::BLOCKS_GEOGRID_AND_RADARGRID);
 
     // RTC
     isce3::io::Raster* rtc_raster = nullptr;
@@ -1688,34 +1688,33 @@ void Geocode<T>::geocodeAreaProj(
             else
                 rtc_memory_mode = isce3::geometry::RTC_BLOCKS_GEOGRID;
 
-            computeRtc(dem_raster, *rtc_raster, radar_grid_cropped, _orbit, _doppler,
-                    _geoGridStartY, _geoGridSpacingY, _geoGridStartX,
+            computeRtc(dem_raster, *rtc_raster, radar_grid_cropped, _orbit,
+                    _doppler, _geoGridStartY, _geoGridSpacingY, _geoGridStartX,
                     _geoGridSpacingX, _geoGridLength, _geoGridWidth, _epsgOut,
                     input_terrain_radiometry, output_terrain_radiometry,
                     rtc_area_mode, rtc_algorithm, rtc_geogrid_upsampling,
                     rtc_min_value_db, radar_grid_nlooks, nullptr, nullptr,
-                    nullptr, rtc_memory_mode, dem_interp_method, _threshold, 
+                    nullptr, rtc_memory_mode, dem_interp_method, _threshold,
                     _numiter, 1.0e-8);
         } else {
             info << "reading pre-computed RTC..." << pyre::journal::newline;
             rtc_raster = input_rtc;
         }
 
         if (is_radar_grid_single_block) {
-            rtc_area.resize(radar_grid_cropped.length(), radar_grid_cropped.width());
-            rtc_raster->getBlock(rtc_area.data(), 0, 0, radar_grid_cropped.width(),
-                                 radar_grid_cropped.length(), 1);
+            rtc_area.resize(
+                    radar_grid_cropped.length(), radar_grid_cropped.width());
+            rtc_raster->getBlock(rtc_area.data(), 0, 0,
+                    radar_grid_cropped.width(), radar_grid_cropped.length(), 1);
         }
     }
 
-
-
-
     // number of bands in the input raster
     info << "nbands: " << nbands << pyre::journal::newline;
 
     info << "radar grid width: " << radar_grid_cropped.width()
-         << ", length: " << radar_grid_cropped.length() << pyre::journal::newline;
+         << ", length: " << radar_grid_cropped.length()
+         << pyre::journal::newline;
 
     info << "geogrid upsampling: " << geogrid_upsampling << pyre::journal::newline;
 
@@ -1765,15 +1764,14 @@ void Geocode<T>::geocodeAreaProj(
 
         // read slant-range image
         if (std::is_same<T, T_out>::value || nbands_off_diag_terms > 0) {
-            _getUpsampledBlock<T, T>(
-                rdrDataT, input_raster, offset_x, offset_y,
-                radar_grid_cropped.width(), radar_grid_cropped.length(),
-                flag_upsample_radar_grid, geocode_memory_mode, info);
+            _getUpsampledBlock<T, T>(rdrDataT, input_raster, offset_x, offset_y,
+                    radar_grid_cropped.width(), radar_grid_cropped.length(),
+                    flag_upsample_radar_grid, geocode_memory_mode, info);
         } else {
-            _getUpsampledBlock<T, T_out>(
-                rdrData, input_raster, offset_x, offset_y, 
-                radar_grid_cropped.width(), radar_grid_cropped.length(), 
-                flag_upsample_radar_grid, geocode_memory_mode, info);
+            _getUpsampledBlock<T, T_out>(rdrData, input_raster, offset_x,
+                    offset_y, radar_grid_cropped.width(),
+                    radar_grid_cropped.length(), flag_upsample_radar_grid,
+                    geocode_memory_mode, info);
         }
     }
     int block_size_x, nblocks_x, block_size_with_upsampling_x;
@@ -1809,41 +1807,43 @@ void Geocode<T>::geocodeAreaProj(
 
     info << "starting geocoding" << pyre::journal::endl;
     if (!std::is_same<T, T_out>::value && nbands_off_diag_terms == 0) {
-        _Pragma("omp parallel for schedule(dynamic)")
-        for (int block_y = 0; block_y < nblocks_y; ++block_y) {
+        _Pragma("omp parallel for schedule(dynamic)") for (int block_y = 0;
+                                                           block_y < nblocks_y;
+                                                           ++block_y)
+        {
             for (int block_x = 0; block_x < nblocks_x; ++block_x) {
-                _runBlock<T_out, T_out>(
-                        radar_grid_cropped, is_radar_grid_single_block, rdrData,
-                        block_size_y, block_size_with_upsampling_y, block_y,
-                        block_size_x, block_size_with_upsampling_x, block_x,
-                        numdone, progress_block, geogrid_upsampling, nbands,
+                _runBlock<T_out, T_out>(radar_grid_cropped,
+                        is_radar_grid_single_block, rdrData, block_size_y,
+                        block_size_with_upsampling_y, block_y, block_size_x,
+                        block_size_with_upsampling_x, block_x, numdone,
+                        progress_block, geogrid_upsampling, nbands,
                         nbands_off_diag_terms, dem_interp_method, dem_raster,
                         out_off_diag_terms, out_geo_rdr, out_geo_dem,
                         out_geo_nlooks, out_geo_rtc, proj.get(), flag_apply_rtc,
                         rtc_raster, input_raster, offset_y, offset_x,
-                        output_raster, rtc_area, 
-                        rtc_min_value, abs_cal_factor, clip_min, clip_max,
-                        min_nlooks, radar_grid_nlooks, flag_upsample_radar_grid,
-                        geocode_memory_mode, info);
+                        output_raster, rtc_area, rtc_min_value, abs_cal_factor,
+                        clip_min, clip_max, min_nlooks, radar_grid_nlooks,
+                        flag_upsample_radar_grid, geocode_memory_mode, info);
             }
         }
     } else {
-        _Pragma("omp parallel for schedule(dynamic)")
-        for (int block_y = 0; block_y < nblocks_y; ++block_y) {
+        _Pragma("omp parallel for schedule(dynamic)") for (int block_y = 0;
+                                                           block_y < nblocks_y;
+                                                           ++block_y)
+        {
             for (int block_x = 0; block_x < nblocks_x; ++block_x) {
-                _runBlock<T, T_out>(
-                        radar_grid_cropped, is_radar_grid_single_block, rdrDataT,
-                        block_size_y, block_size_with_upsampling_y, block_y,
-                        block_size_x, block_size_with_upsampling_x, block_x,
-                        numdone, progress_block, geogrid_upsampling, nbands,
+                _runBlock<T, T_out>(radar_grid_cropped,
+                        is_radar_grid_single_block, rdrDataT, block_size_y,
+                        block_size_with_upsampling_y, block_y, block_size_x,
+                        block_size_with_upsampling_x, block_x, numdone,
+                        progress_block, geogrid_upsampling, nbands,
                         nbands_off_diag_terms, dem_interp_method, dem_raster,
                         out_off_diag_terms, out_geo_rdr, out_geo_dem,
                         out_geo_nlooks, out_geo_rtc, proj.get(), flag_apply_rtc,
                         rtc_raster, input_raster, offset_y, offset_x,
-                        output_raster, rtc_area, 
-                        rtc_min_value, abs_cal_factor, clip_min, clip_max,
-                        min_nlooks, radar_grid_nlooks, flag_upsample_radar_grid,
-                        geocode_memory_mode, info);
+                        output_raster, rtc_area, rtc_min_value, abs_cal_factor,
+                        clip_min, clip_max, min_nlooks, radar_grid_nlooks,
+                        flag_upsample_radar_grid, geocode_memory_mode, info);
             }
         }
     }
@@ -2333,31 +2333,28 @@ void Geocode<T>::_runBlock(
                     out_geo_dem_array, out_geo_nlooks, out_geo_nlooks_array,
                     out_geo_rtc, out_geo_rtc_array);
 
-            isce3::core::Matrix<T_out> geoDataBlock(this_block_size_y,
-                                                    this_block_size_x);
+            isce3::core::Matrix<T_out> geoDataBlock(
+                    this_block_size_y, this_block_size_x);
 
             // fill both matrices with NaN
             geoDataBlock.fill(nan_t_out);
 
             for (int band = 0; band < nbands; ++band) {
                 _Pragma("omp critical")
                 {
-                    output_raster.setBlock(
-                            geoDataBlock.data(), block_x * block_size_x,
-                            block_y * block_size_y, this_block_size_x,
-                            this_block_size_y, band + 1);
+                    output_raster.setBlock(geoDataBlock.data(),
+                            block_x * block_size_x, block_y * block_size_y,
+                            this_block_size_x, this_block_size_y, band + 1);
                 }
             }
 
-
             if (nbands_off_diag_terms > 0) {
                 for (int band = 0; band < nbands_off_diag_terms; ++band) {
                     _Pragma("omp critical")
                     {
-                        out_off_diag_terms->setBlock(
-                                geoDataBlock.data(), block_x * block_size_x,
-                                block_y * block_size_y, this_block_size_x,
-                                this_block_size_y, band + 1);
+                        out_off_diag_terms->setBlock(geoDataBlock.data(),
+                                block_x * block_size_x, block_y * block_size_y,
+                                this_block_size_x, this_block_size_y, band + 1);
                     }
                 }
             }
@@ -2371,16 +2368,14 @@ void Geocode<T>::_runBlock(
 
         if (flag_apply_rtc) {
 
-            rtc_area_block.resize(radar_grid_block.length(),
-                                  radar_grid_block.width());
+            rtc_area_block.resize(
+                    radar_grid_block.length(), radar_grid_block.width());
             rtc_raster->getBlock(rtc_area_block.data(), offset_x, offset_y,
-                                 radar_grid_block.width(),
-                                 radar_grid_block.length(), 1);
+                    radar_grid_block.width(), radar_grid_block.length(), 1);
         }
 
-        _getUpsampledBlock<T, T2>(
-                rdrDataBlock, input_raster, offset_x + raster_offset_x,
-                offset_y + raster_offset_y,
+        _getUpsampledBlock<T, T2>(rdrDataBlock, input_raster,
+                offset_x + raster_offset_x, offset_y + raster_offset_y,
                 radar_grid_block.width(), radar_grid_block.length(),
                 flag_upsample_radar_grid, geocode_memory_mode, info);
     }
@@ -2452,8 +2447,7 @@ void Geocode<T>::_runBlock(
 
             const int jj = block_x * block_size_with_upsampling_x + j;
 
-            _Pragma("omp atomic")
-            numdone++;
+            _Pragma("omp atomic") numdone++;
             if (numdone % progress_block == 0)
                 _Pragma("omp critical")
                 {

cxx/isce3/geocode/GeocodeCov.h

@@ -37,6 +37,26 @@ enum geocodeMemoryMode {
     BLOCKS_GEOGRID_AND_RADARGRID = 3
 };
 
+constexpr static long long DEFAULT_MAX_BLOCK_SIZE = 1 << 29; // 512MB
+/** Compute the number of blocks and their length for block processing
+ * in the length direction
+ *
+ * @param[in]  array_length        Length of the data to be processed
+ * @param[in]  array_width         Width of the data to be processed
+ * @param[in]  nbands              Number of the bands to be processed
+ * @param[in]  type_size           Type size of the data to be processed
+ * @param[in]  channel             Pyre info channel
+ * @param[out] block_length        Block length
+ * @param[out] nblock_y            Number of blocks in the Y direction
+ * @param[in]  max_block_size      Maximum block size in Bytes (per thread)
+ */
+void getBlocksNumberAndLength(const int array_length, const int array_width,
+        const int nbands = 1,
+        const size_t type_size = 4, // Float32
+        pyre::journal::info_t* channel = nullptr, int* block_length = nullptr,
+        int* nblock_y = nullptr,
+        const long long max_block_size = DEFAULT_MAX_BLOCK_SIZE);
+
 template<class T>
 class Geocode {
 public:
@@ -309,26 +329,6 @@ class Geocode {
     void updateGeoGrid(const isce3::product::RadarGridParameters& radar_grid,
                        isce3::io::Raster& dem_raster);
 
-    constexpr static long long DEFAULT_MAX_BLOCK_SIZE = 1 << 29; // 512MB
-    /** Compute the number of blocks and their length for block processing
-     * in the length direction
-     *
-     * @param[in]  array_length        Length of the data to be processed
-     * @param[in]  array_width         Width of the data to be processed
-     * @param[in]  nbands              Number of the bands to be processed
-     * @param[in]  type_size           Type size of the data to be processed
-     * @param[in]  channel             Pyre info channel
-     * @param[out] block_length        Block length
-     * @param[out] nblock_y            Number of blocks in the Y direction
-     * @param[in]  max_block_size      Maximum block size in Bytes (per thread)
-     */
-    void getBlocksNumberAndLength(const int array_length, const int array_width,
-            const int nbands = 1,
-            const size_t type_size = 4, // Float32
-            pyre::journal::info_t* channel = nullptr,
-            int* block_length = nullptr, int* nblock_y = nullptr,
-            const long long max_block_size = DEFAULT_MAX_BLOCK_SIZE);
-
     // Get/set data interpolator
     isce3::core::dataInterpMethod dataInterpolator() const 
     { 
@@ -450,33 +450,27 @@ class Geocode {
             isce3::geometry::DEMInterpolator& dem_interp);
 
     template<class T2, class T_out>
-    void _runBlock(
-            const isce3::product::RadarGridParameters& radar_grid,
+    void _runBlock(const isce3::product::RadarGridParameters& radar_grid,
             bool is_radar_grid_single_block,
             std::vector<std::unique_ptr<isce3::core::Matrix<T2>>>& rdrData,
             int block_size_y, int block_size_with_upsampling_y, int block_y,
             int block_size_x, int block_size_with_upsampling_x, int block_x,
-            long long& numdone, const long long&  progress_block, 
-            double geogrid_upsampling,
-            int nbands, int nbands_off_diag_terms,
+            long long& numdone, const long long& progress_block,
+            double geogrid_upsampling, int nbands, int nbands_off_diag_terms,
             isce3::core::dataInterpMethod dem_interp_method,
             isce3::io::Raster& dem_raster,
             isce3::io::Raster* out_off_diag_terms,
-            isce3::io::Raster* out_geo_rdr,
-            isce3::io::Raster* out_geo_dem,
+            isce3::io::Raster* out_geo_rdr, isce3::io::Raster* out_geo_dem,
             isce3::io::Raster* out_geo_nlooks, isce3::io::Raster* out_geo_rtc,
-            isce3::core::ProjectionBase* proj, 
-            bool flag_apply_rtc, isce3::io::Raster* rtc_raster,
-            isce3::io::Raster& input_raster,
+            isce3::core::ProjectionBase* proj, bool flag_apply_rtc,
+            isce3::io::Raster* rtc_raster, isce3::io::Raster& input_raster,
             int raster_offset_y, int raster_offset_x,
             isce3::io::Raster& output_raster,
-            isce3::core::Matrix<float>& rtc_area,
-            float rtc_min_value,
+            isce3::core::Matrix<float>& rtc_area, float rtc_min_value,
             double abs_cal_factor, float clip_min, float clip_max,
             float min_nlooks, float radar_grid_nlooks,
             bool flag_upsample_radar_grid,
-            geocodeMemoryMode geocode_memory_mode, 
-            pyre::journal::info_t& info);
+            geocodeMemoryMode geocode_memory_mode, pyre::journal::info_t& info);
 
     void _loadDEM(isce3::io::Raster& demRaster,
                   isce3::geometry::DEMInterpolator& demInterp,

cxx/isce3/math/ComplexMultiply.h

@@ -2,15 +2,20 @@
 
 #include <complex>
 
-/** Promotion rules so T * complex<U> works for {T,U} in {float,double}. */
 namespace isce3 { namespace math { namespace complex_multiply {
 
-inline auto operator*(float t, std::complex<double> u) {
-    return double(t) * u;
+template<typename T1, typename T2>
+auto operator*(const std::complex<T1>& lhs, const T2& rhs)
+{
+    using U = typename std::common_type_t<T1, T2>;
+    return std::complex<U>(lhs) * U(rhs);
 }
 
-inline auto operator*(std::complex<double> t, float u) {
-    return t * double(u);
+template<typename T1, typename T2>
+auto operator*(const T1& lhs, const std::complex<T2>& rhs)
+{
+    using U = typename std::common_type_t<T1, T2>;
+    return U(lhs) * std::complex<U>(rhs);
 }
 
 }}}