Update geocode polygon (#518)

* moving files geocode/Geocode.* to geocode/GeocodeCov.*

* moving cython files geocode/Geocode.* to geocode/GeocodeCov.*

* updating c++ modules, cython- and pybind11 bindinds to new location for the Geocode module

* fixing RTC memory mode in the RTC cython binding

* fix bindings and tests

* Substitute function getGeoAreaElementMean by class GeocodePolygon and update respective bindings

* update branch after Geoff's reviews in #516

* fix unit tests

* update search for min. and max. height (to define the geogrid bbox)

* update GeocodePolygon geo2rdr threshold defaults

* substitute indexes by indices

* add pybind11 bindings, remove empty constructor, and refactor c++ code

* update docstrings

* refactor GeocodePolygon

* substitute pairs of std::min/std::max by std::minmax_element

* fix GeocodePolygon docstrings

* default flag_apply_rtc to false, substitute rtc_min_val_db by rtc_min_value_db

* simplify multiplication of abs_cal_factor / nlooks

* remove second (unnecessary) validation of GeocodePolygon attributes

* revert changes. _ValidatePolygon needs to run inside getPolygonMean because radar_grid may change

* separate GeocodePolygon from GeocodeCov compilation unit

* remove unnecessary imports

* remove unnecessary imports (pybinds)

* remove commented code

* move geocoded shared function from GeocodeCov.icc to GeocodeHelpers.h

* do not remove trailing newline in GeocodeCov.icc

* Update cxx/isce3/geocode/GeocodeHelpers.h

Co-Authored-By: Ryan Burns <[email protected]>

Co-authored-by: Ryan Burns <[email protected]>

Author: 2 people

cxx/isce3/Headers.cmake

@@ -89,6 +89,7 @@ geometry/Geo2rdr.h
 geometry/Geo2rdr.icc
 geocode/GeocodeCov.h
 geocode/GeocodeCov.icc
+geocode/GeocodePolygon.h
 geometry/geometry.h
 geometry/RTC.h
 geometry/Topo.h

cxx/isce3/Sources.cmake

@@ -43,6 +43,7 @@ geocode/loadDem.cpp
 geometry/DEMInterpolator.cpp
 geometry/Geo2rdr.cpp
 geocode/GeocodeCov.cpp
+geocode/GeocodePolygon.cpp
 geometry/geometry.cpp
 geometry/RTC.cpp
 geometry/Topo.cpp

cxx/isce3/geocode/GeocodeCov.cpp

@@ -1,4 +1,5 @@
 #include "GeocodeCov.h"
+#include "GeocodeHelpers.h"
 
 #include <algorithm>
 #include <cmath>
@@ -22,48 +23,6 @@ using isce3::core::Vec3;
 
 namespace isce3 { namespace geocode {
 
-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);
-}
-
-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);
-}
-
-template<typename T, typename T_out>
-void _convertToOutputType(T a, T_out& b)
-{
-    b = a;
-}
-
-template<typename T, typename T_out>
-void _convertToOutputType(std::complex<T> a, T_out& b)
-{
-    b = std::norm(a);
-}
-
-template<typename T, typename T_out>
-void _convertToOutputType(std::complex<T> a, std::complex<T_out>& b)
-{
-    b = a;
-}
-
-template<typename T, typename T_out>
-void _accumulate(T_out& band_value, T a, double b)
-{
-    if (b == 0)
-        return;
-    T_out a2;
-    _convertToOutputType(a, a2);
-    band_value += a2 * b;
-}
-
 template<class T>
 void Geocode<T>::updateGeoGrid(
         const isce3::product::RadarGridParameters& radar_grid,
@@ -1903,312 +1862,6 @@ template class Geocode<double>;
 template class Geocode<std::complex<float>>;
 template class Geocode<std::complex<double>>;
 
-// template <typename T>
-std::vector<float> getGeoAreaElementMean(
-        const std::vector<double>& x_vect, const std::vector<double>& y_vect,
-        const isce3::product::RadarGridParameters& radar_grid,
-        const isce3::core::Orbit& orbit,
-        const isce3::core::LUT2d<double>& input_dop,
-        isce3::io::Raster& input_raster, isce3::io::Raster& dem_raster,
-        bool flag_apply_rtc,
-        isce3::geometry::rtcInputTerrainRadiometry input_terrain_radiometry, 
-        isce3::geometry::rtcOutputTerrainRadiometry output_terrain_radiometry,  
-        int exponent, geocodeOutputMode output_mode, double geogrid_upsampling,
-        float rtc_min_value_db, double abs_cal_factor, float radar_grid_nlooks,
-        float* out_nlooks, isce3::core::dataInterpMethod dem_interp_method,
-        double threshold, int num_iter, double delta_range)
-{
-
-    pyre::journal::info_t info("isce.geometry.getGeoAreaElementMean");
-
-    if (std::isnan(geogrid_upsampling))
-        geogrid_upsampling = 1;
-    assert(geogrid_upsampling > 0);
-
-    const double x0 = std::min_element(x_vect.begin(), x_vect.end())[0];
-    const double xf = std::max_element(x_vect.begin(), x_vect.end())[0];
-
-    const double y0 = std::min_element(y_vect.begin(), y_vect.end())[0];
-    const double yf = std::max_element(y_vect.begin(), y_vect.end())[0];
-
-    const double dx = dem_raster.dx();
-    const double dy = dem_raster.dy();
-
-    std::string input_terrain_radiometry_str =
-            get_input_terrain_radiometry_str(input_terrain_radiometry);
-
-    info << "input radiometry: " << input_terrain_radiometry_str
-         << pyre::journal::newline << "look side: " << radar_grid.lookSide()
-         << pyre::journal::newline
-         << "radar_grid length: " << radar_grid.length()
-         << ", width: " << radar_grid.width() << pyre::journal::newline
-         << "RTC min value [dB]: " << rtc_min_value_db << pyre::journal::endl;
-
-    const double margin_x = std::abs(dx) * 10;
-    const double margin_y = std::abs(dy) * 10;
-    isce3::geometry::DEMInterpolator dem_interp;
-
-    dem_interp.loadDEM(dem_raster, x0 - margin_x, xf + margin_x,
-                       std::min(y0, yf) - margin_y,
-                       std::max(y0, yf) + margin_y);
-
-    isce3::core::Ellipsoid ellipsoid =
-            isce3::core::Ellipsoid(isce3::core::EarthSemiMajorAxis,
-                                   isce3::core::EarthEccentricitySquared);
-
-    int epsg = dem_raster.getEPSG();
-    std::unique_ptr<isce3::core::ProjectionBase> proj(
-            isce3::core::createProj(epsg));
-    const double pixazm = radar_grid.azimuthTimeInterval();
-    const double start = radar_grid.sensingStart() - 0.5 * pixazm;
-    const double dr = radar_grid.rangePixelSpacing();
-    double r0 = radar_grid.startingRange() - 0.5 * dr;
-
-    double a = radar_grid.sensingMid();
-    double r = radar_grid.midRange();
-
-    if (x_vect.size() != y_vect.size()) {
-        std::string error_message =
-                "ERROR x and y vectors have a different number of elements.";
-        throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_message);
-    }
-    int n_elements = x_vect.size();
-    std::vector<double> a_vect, r_vect;
-    a_vect.reserve(n_elements);
-    r_vect.reserve(n_elements);
-
-    info << "polygon indexes (a, r): ";
-
-    for (int i = 0; i < n_elements; ++i) {
-
-        const double x = x_vect[i];
-        const double y = y_vect[i];
-
-        const Vec3 dem11 = {x, y, dem_interp.interpolateXY(x, y)};
-        int converged = isce3::geometry::geo2rdr(
-                proj->inverse(dem11), ellipsoid, orbit, input_dop, a, r,
-                radar_grid.wavelength(), radar_grid.lookSide(), threshold,
-                num_iter, delta_range);
-        if (!converged) {
-            info << "WARNING convergence not found for vertex (x, y): " << x
-                 << ", " << y << pyre::journal::endl;
-            continue;
-        }
-        double a_index = (a - start) / pixazm;
-        double r_index = (r - r0) / dr;
-
-        info << "(" << a_index << ", " << r_index << "), ";
-
-        a_vect.emplace_back(a_index);
-        r_vect.emplace_back(r_index);
-    }
-
-    info << pyre::journal::endl;
-
-    const double a_min = std::min_element(a_vect.begin(), a_vect.end())[0];
-    const double a_max = std::max_element(a_vect.begin(), a_vect.end())[0];
-
-    const double r_min = std::min_element(r_vect.begin(), r_vect.end())[0];
-    const double r_max = std::max_element(r_vect.begin(), r_vect.end())[0];
-
-    const int y_min = std::max(0, (int) std::floor(a_min));
-    const int x_min = std::max(0, (int) std::floor(r_min));
-    const int ysize =
-            std::min((int) radar_grid.length(), (int) std::ceil(a_max)) - y_min;
-    const int xsize =
-            std::min((int) radar_grid.width(), (int) std::ceil(r_max)) - x_min;
-
-    info << "cropping radar grid from index (a0: " << y_min;
-    info << ", r0: " << x_min << ") to index (af: " << y_min + ysize;
-    info << ", rf: " << x_min + xsize << ")" << pyre::journal::endl;
-
-    isce3::product::RadarGridParameters radar_grid_cropped =
-            radar_grid.offsetAndResize(y_min, x_min, ysize, xsize);
-
-    info << "cropped radar_grid length: " << radar_grid_cropped.length()
-         << ", width: " << radar_grid_cropped.width() << pyre::journal::newline;
-
-    if (output_mode == geocodeOutputMode::INTERP) {
-        std::string error_msg = "invalid option";
-        throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg);
-    }
-
-    if (abs_cal_factor != 1)
-        info << "absolute calibration factor: " << abs_cal_factor
-             << pyre::journal::endl;
-
-    if (flag_apply_rtc) {
-
-        std::string input_terrain_radiometry_str =
-                get_input_terrain_radiometry_str(input_terrain_radiometry);
-        info << "input radiometry: " << input_terrain_radiometry_str
-             << pyre::journal::endl;
-    }
-
-    isce3::core::Matrix<float> rtc_area;
-    std::unique_ptr<isce3::io::Raster> rtc_raster_unique_ptr;
-    if (flag_apply_rtc) {
-
-        info << "computing RTC area factor..." << pyre::journal::endl;
-        rtc_raster_unique_ptr = std::make_unique<isce3::io::Raster>(
-                "/vsimem/dummy", radar_grid_cropped.width(),
-                radar_grid_cropped.length(), 1, GDT_Float32, "ENVI");
-        isce3::geometry::rtcAreaMode rtc_area_mode =
-                isce3::geometry::rtcAreaMode::AREA_FACTOR;
-        isce3::geometry::rtcAlgorithm rtc_algorithm =
-                isce3::geometry::rtcAlgorithm::RTC_AREA_PROJECTION;
-
-        isce3::geometry::rtcMemoryMode rtc_memory_mode =
-                isce3::geometry::rtcMemoryMode::RTC_SINGLE_BLOCK;
-
-        computeRtc(radar_grid_cropped, orbit, input_dop, dem_raster,
-                *rtc_raster_unique_ptr.get(), input_terrain_radiometry,
-                output_terrain_radiometry, rtc_area_mode, rtc_algorithm,
-                geogrid_upsampling * 2, rtc_min_value_db, radar_grid_nlooks,
-                nullptr, rtc_memory_mode, dem_interp_method, threshold, num_iter,
-                delta_range);
-
-        rtc_area.resize(radar_grid_cropped.length(),
-                        radar_grid_cropped.width());
-
-        rtc_raster_unique_ptr->getBlock(rtc_area.data(), 0, 0,
-                                              radar_grid_cropped.width(),
-                                              radar_grid_cropped.length(), 1);
-
-        info << "... done (RTC) " << pyre::journal::endl;
-    }
-
-    double rtc_min_value = 0;
-    if (!std::isnan(rtc_min_value_db) && flag_apply_rtc) {
-        rtc_min_value = std::pow(10., (rtc_min_value_db / 10.));
-        info << "RTC min. value: " << rtc_min_value_db
-             << " [dB] = " << rtc_min_value << pyre::journal::endl;
-    }
-
-    GDALDataType input_dtype = input_raster.dtype();
-    if (exponent == 0 && GDALDataTypeIsComplex(input_dtype))
-        exponent = 2;
-
-    if (input_raster.dtype() == GDT_Float32) {
-        info << "dtype: GDT_Float32" << pyre::journal::endl;
-        return _getGeoAreaElementMean<float>(
-                r_vect, a_vect, x_min, y_min, flag_apply_rtc,
-                rtc_area, radar_grid_cropped,
-                input_raster, rtc_min_value, out_nlooks,
-                abs_cal_factor, radar_grid_nlooks);
-    } else if (input_raster.dtype() == GDT_CFloat32) {
-        info << "dtype: GDT_CFloat32" << pyre::journal::endl;
-        return _getGeoAreaElementMean<std::complex<float>>(
-                r_vect, a_vect, x_min, y_min, flag_apply_rtc,
-                rtc_area, radar_grid_cropped,
-                input_raster, rtc_min_value, out_nlooks,
-                abs_cal_factor, radar_grid_nlooks);
-    } else
-        info << "ERROR not implemented for datatype: " << input_raster.dtype()
-             << pyre::journal::endl;
-
-    std::vector<float> empty_vector;
-    return empty_vector;
-}
-
-template<typename T>
-std::vector<float> _getGeoAreaElementMean(
-        const std::vector<double>& r_vect, const std::vector<double>& a_vect,
-        int x_min, int y_min, bool flag_apply_rtc,
-        isce3::core::Matrix<float>& rtc_area,
-        const isce3::product::RadarGridParameters& radar_grid,
-        isce3::io::Raster& input_raster, 
-        float rtc_min_value, float* out_nlooks, double abs_cal_factor,
-        float radar_grid_nlooks)
-{
-
-    pyre::journal::info_t info("isce.geometry._getGeoAreaElementMean");
-
-    // number of bands in the input raster
-    const int nbands = input_raster.numBands();
-    const int size_y = radar_grid.length();
-    const int size_x = radar_grid.width();
-    info << "nbands: " << nbands << pyre::journal::endl;
-
-    std::vector<std::unique_ptr<isce3::core::Matrix<T>>> rdrDataBlock;
-    rdrDataBlock.reserve(nbands);
-
-    for (int band = 0; band < nbands; ++band) {
-        if (nbands == 1)
-            info << "loading slant-range image..." << pyre::journal::endl;
-        else
-            info << "loading slant-range band: " << band << pyre::journal::endl;
-        rdrDataBlock.emplace_back(
-                std::make_unique<isce3::core::Matrix<T>>(size_y, size_x));
-
-        input_raster.getBlock(rdrDataBlock[band]->data(), x_min, y_min,
-                              size_x, size_y, band + 1);
-    }
-
-    isce3::core::Matrix<double> w_arr(size_y, size_x);
-    w_arr.fill(0);
-
-    int plane_orientation;
-    if (radar_grid.lookSide() == isce3::core::LookSide::Left)
-        plane_orientation = -1;
-    else
-        plane_orientation = 1;
-
-    double w_total = 0;
-    for (int i = 0; i < r_vect.size(); ++i) {
-
-        double x00, y00, x01, y01;
-        y00 = a_vect[i];
-        x00 = r_vect[i];
-
-        if (i < r_vect.size() - 1) {
-            y01 = a_vect[i + 1];
-            x01 = r_vect[i + 1];
-        } else {
-            y01 = a_vect[0];
-            x01 = r_vect[0];
-        }
-        isce3::geometry::areaProjIntegrateSegment(
-                y00 - y_min, y01 - y_min, x00 - x_min, x01 - x_min, size_y,
-                size_x, w_arr, w_total, plane_orientation);
-    }
-    std::vector<float> cumulative_sum(nbands);
-    double nlooks = 0;
-
-    Geocode<T> geo_obj;
-    for (int y = 0; y < size_y; ++y)
-        for (int x = 0; x < size_x; ++x) {
-            double w = w_arr(y, x);
-            if (w == 0)
-                continue;
-            if (flag_apply_rtc) {
-                const float rtc_value = rtc_area(y, x);
-                if (std::isnan(rtc_value) || rtc_value < rtc_min_value)
-                    continue;
-                nlooks += w;
-                w /= rtc_value;
-            } else {
-                nlooks += w;
-            }
-
-            for (int band = 0; band < nbands; ++band)
-                _accumulate(cumulative_sum[band],
-                            rdrDataBlock[band]->operator()(y, x), w);
-        }
-
-    info << "nlooks: " << radar_grid_nlooks * std::abs(nlooks)
-         << pyre::journal::endl;
-    for (int band = 0; band < nbands; ++band) {
-        cumulative_sum[band] = (cumulative_sum[band] * abs_cal_factor / nlooks);
-        info << "mean value (band = " << band << "): " << cumulative_sum[band]
-             << pyre::journal::endl;
-    }
-
-    if (out_nlooks != nullptr) {
-        *out_nlooks = radar_grid_nlooks * std::abs(nlooks);
-    }
-
-    return cumulative_sum;
-}
 
-}} // namespace isce3::geocode
+} // namespace geocode
+} // namespace isce3