Resample v2 flattening

cxx/isce3/Headers.cmake

@@ -114,6 +114,7 @@ geometry/metadataCubes.h
 geogrid/getRadarGrid.h
 geogrid/relocateRaster.h
 image/forward.h
+image/Flatten.h
 image/Resample.h
 image/ResampSlc.h
 image/ResampSlc.icc

cxx/isce3/Sources.cmake

@@ -59,6 +59,7 @@ geometry/TopoLayers.cpp
 geometry/metadataCubes.cpp
 geogrid/getRadarGrid.cpp
 geogrid/relocateRaster.cpp
+image/Flatten.cpp
 image/Resample.cpp
 image/ResampSlc.cpp
 io/gdal/Dataset.cpp

cxx/isce3/image/Flatten.cpp

@@ -0,0 +1,142 @@
+
+#include "Flatten.h"
+
+#include <isce3/product/RadarGridProduct.h>
+
+
+namespace isce3::image::flatten {
+
+
+/** Calculate the phase flattening parameter for a given pixel
+ * 
+ * This code assumes that the output and input grids have the same wavelength
+ * and pixel spacing.
+ * 
+ * @returns     double      The flattening phase for this pixel.
+ *
+ * RadarGridParameters of output radar data
+ * @param[in] radarGridOut
+ * RadarGridParameters of input radar data
+ * @param[in] radarGridIn
+ * The difference between the input range index and its point on the output grid;
+ * unit: range column indices (int)
+ * @param[in] rangeOffset
+*/
+double _pixelFlatteningPhase(
+    const isce3::product::RadarGridParameters& radarGridOut,
+    const isce3::product::RadarGridParameters& radarGridIn,
+    const double rangeOffset
+)
+{
+    // Values from the output radar grid
+    // unit: distance (meters)
+    const auto startingRange = radarGridOut.startingRange();
+    // unit: distance (meters) per range index position
+    const auto rangePixelSpacing = radarGridOut.rangePixelSpacing();
+    // unit: distance (meters) per cycle
+    const auto wavelength = radarGridOut.wavelength();
+
+    // Values from the input radar grid
+    // unit: distance (meters)
+    const auto inStartingRange = radarGridIn.startingRange();
+    // The number of cycles completed by this wavelength of light in a given unit
+    // distance.
+    // unit: cycles per unit distance
+    const auto spatialFrequency = 1. / wavelength;
+
+    // The difference between the starting range and the input starting range.
+    // unit: distance (meters)
+    const auto startingRangeDifference = startingRange - inStartingRange;
+    // The distance from the point referenced at the range pixel grid and the point that
+    // its offset points to.
+    // unit: distance (meters)
+    const auto offsetDistance = rangeOffset * rangePixelSpacing;
+    // The total distance difference between the indexed point on the output radar
+    // grid and its offset point on the input radar grid.
+    // unit: distance (meters)
+    const auto distanceDiff = startingRangeDifference + offsetDistance;
+
+    // The number of cycles covered over that distance.
+    // unit: cycles
+    const auto cycleDiff = distanceDiff * spatialFrequency;
+
+    // The flattening phase parameter for this pixel is the cycles in a two-way trip,
+    // multiplied by 2 * pi.
+    // unit: radians
+    return 4. * M_PI * cycleDiff;
+}
+
+
+void flattenAtCoords(
+    ArrayRef2D<std::complex<float>> dataBlock,
+    const ArrayRefConst2D<double> rangeIndices,
+    const isce3::product::RadarGridParameters& radarGridOut,
+    const isce3::product::RadarGridParameters& radarGridIn,
+    const size_t outRgFirstPixel
+)
+{
+    const size_t outWidth = dataBlock.cols();
+    const size_t outLength = dataBlock.rows();
+
+#pragma omp parallel for collapse(2)
+    for (size_t azIndexIn = 0; azIndexIn < outLength; ++azIndexIn){
+        for (size_t rgIndexIn = 0; rgIndexIn < outWidth; ++rgIndexIn){
+
+            // Get the range indices on the output grid corresponding to these indices
+            // on the input grid overall.
+            const double rgIndexOut = rangeIndices(azIndexIn, rgIndexIn);
+
+            double rangeOffset = rgIndexOut - static_cast<double>(rgIndexIn) -
+                static_cast<double>(outRgFirstPixel);
+
+            const double flattenPhase = _pixelFlatteningPhase(
+                radarGridOut, radarGridIn, rangeOffset
+            );
+
+            // Update dataBlock column and row from index
+            const std::complex<float> cpxVal(
+                std::cos(flattenPhase), std::sin(flattenPhase)
+            );
+            dataBlock(azIndexIn, rgIndexIn) *= cpxVal;
+        }
+    } // end multithreaded block
+}
+
+
+void getFlatteningPhase(
+    ArrayRef2D<std::complex<float>> dataBlock,
+    const ArrayRefConst2D<double> rangeIndices,
+    const isce3::product::RadarGridParameters& radarGridOut,
+    const isce3::product::RadarGridParameters& radarGridIn,
+    const size_t outRgFirstPixel
+)
+{
+    const size_t outWidth = dataBlock.cols();
+    const size_t outLength = dataBlock.rows();
+
+#pragma omp parallel for collapse(2)
+    for (size_t azIndexIn = 0; azIndexIn < outLength; ++azIndexIn){
+        for (size_t rgIndexIn = 0; rgIndexIn < outWidth; ++rgIndexIn){
+
+            // Get the range indices on the output grid corresponding to these indices
+            // on the input grid overall.
+            const double rgIndexOut = rangeIndices(azIndexIn, rgIndexIn);
+
+            double rangeOffset = rgIndexOut - static_cast<double>(rgIndexIn) -
+                static_cast<double>(outRgFirstPixel);
+
+            const double flattenPhase = _pixelFlatteningPhase(
+                radarGridOut, radarGridIn, rangeOffset
+            );
+
+            // Update dataBlock column and row from index
+            const std::complex<float> cpxVal(
+                std::cos(flattenPhase), std::sin(flattenPhase)
+            );
+            dataBlock(azIndexIn, rgIndexIn) = cpxVal;
+        }
+    } // end multithreaded block
+}
+
+
+}  // namespace isce3::image::flatten

cxx/isce3/image/Flatten.h

@@ -0,0 +1,59 @@
+
+
+#include <complex>
+#include <Eigen/Core>
+#include <limits>
+
+#include <isce3/product/RadarGridParameters.h>
+
+namespace isce3::image::flatten {
+
+template<class T, int Options = Eigen::RowMajor>
+using Array2D = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Options>;
+
+template<class T, int Options = Eigen::RowMajor>
+using ArrayRef2D = Eigen::Ref<Array2D<T, Options>>;
+
+template<class T, int Options = Eigen::RowMajor>
+using ArrayRefConst2D = Eigen::Ref<const Array2D<T, Options>>;
+
+/**
+ * Re-flatten a grid of SLC data from its' original grid parameters into a new set of
+ * grid parameters.
+ *
+ * @param[out] dataBlock        The SLC data to flatten
+ * @param[in] rangeIndices      range index of each coordinate pixel in the data block
+ *                              in the coordinate system of the alternate radar grid
+ * @param[in] radarGridOut      radar grid parameters of the alternate grid
+ * @param[in] radarGridIn       radar grid parameters of the original grid
+ * @param[in] outRgFirstPixel   range index of the first sample of the alternate grid
+ */
+void flattenAtCoords(
+    ArrayRef2D<std::complex<float>> dataBlock,
+    const ArrayRefConst2D<double> rangeIndices,
+    const isce3::product::RadarGridParameters& radarGridOut,
+    const isce3::product::RadarGridParameters& radarGridIn,
+    const size_t outRgFirstPixel
+);
+
+/**
+ * Get the relative flattening phase, in complex format, necessary to re-flatten a radar
+ * block from one radar grid to another.
+ *
+ * @param[out] dataBlock        The data block to write the flattening phase to.
+ *                              Any values in this block may be overwritten.
+ * @param[in] rangeIndices      range index of each coordinate pixel in the data block
+ *                              in the coordinate system of the alternate radar grid
+ * @param[in] radarGridOut      radar grid parameters of the alternate grid
+ * @param[in] radarGridIn       radar grid parameters of the original grid
+ * @param[in] outRgFirstPixel   range index of the first sample of the alternate grid
+ */
+void getFlatteningPhase(
+    ArrayRef2D<std::complex<float>> dataBlock,
+    const ArrayRefConst2D<double> rangeIndices,
+    const isce3::product::RadarGridParameters& radarGridOut,
+    const isce3::product::RadarGridParameters& radarGridIn,
+    const size_t outRgFirstPixel
+);
+
+} // namespace isce3::image::flatten

python/extensions/pybind_isce3/Sources.cmake

@@ -58,6 +58,7 @@ geogrid/getRadarGrid.cpp
 geogrid/relocateRaster.cpp
 geogrid/geogrid.cpp
 geometry/lookIncFromSr.cpp
+image/Flatten.cpp
 image/image.cpp
 image/Resample.cpp
 image/ResampSlc.cpp

python/extensions/pybind_isce3/image/Flatten.cpp

@@ -0,0 +1,79 @@
+#include "Flatten.h"
+
+#include <isce3/image/Flatten.h>
+#include <isce3/product/RadarGridParameters.h>
+
+#include <pybind11/complex.h>
+#include <pybind11/eigen.h>
+#include <pybind11/numpy.h>
+
+namespace py = pybind11;
+
+
+void addbindings_modulate(py::module & m)
+{
+    m.def(
+        "_flatten_at_coords",
+        &isce3::image::flatten::flattenAtCoords,
+        py::arg("data_block"),
+        py::arg("range_indices"),
+        py::arg("radar_grid_out"),
+        py::arg("radar_grid_in"),
+        py::arg("out_rg_first_pixel"),
+        R"(
+        Re-flatten a grid of SLC data from its' original grid parameters into a new set
+        of grid parameters.
+
+        This function does not check the sizes of the input grids against each other,
+        and should be considered an implementation. See the Python function
+        isce3.image.flatten.flatten_at_coords for the full implementation.
+
+        Parameters
+        ----------
+        data_block : np.ndarray (complex64)
+            The SLC data to flatten
+        range_indices : np.ndarray (float64)
+            range index of each coordinate pixel in the data block in the coordinate
+            system of the alternate radar grid
+        radar_grid_out : isce3.product.RadarGridParameters
+            radar grid parameters of the alternate grid
+        radar_grid_in : isce3.product.RadarGridParameters
+            radar grid parameters of the original grid
+        out_rg_first_pixel : int
+            range index of the first sample of the alternate grid
+        )"
+    );
+
+    m.def(
+        "_get_flattening_phase_at_coords",
+        &isce3::image::flatten::getFlatteningPhase,
+        py::arg("data_block"),
+        py::arg("range_indices"),
+        py::arg("radar_grid_out"),
+        py::arg("radar_grid_in"),
+        py::arg("out_rg_first_pixel"),
+        R"(
+        Re-flatten a grid of SLC data from its' original grid parameters into a new set
+        of grid parameters.
+
+        This function does not check the sizes of the input grids against each other,
+        and should be considered an implementation. See the Python function
+        isce3.image.flatten.get_flattening_phase_at_coords for the full implementation.
+
+        Parameters
+        ----------
+        data_block : np.ndarray (complex64)
+            The data block to write the flattening phase to. Any values in this block
+            may be overwritten.
+        range_indices : np.ndarray (float64)
+            range index of each coordinate pixel in the data block in the coordinate
+            system of the alternate radar grid
+        radar_grid_out : isce3.product.RadarGridParameters
+            radar grid parameters of the alternate grid
+        radar_grid_in : isce3.product.RadarGridParameters
+            radar grid parameters of the original grid
+        out_rg_first_pixel : int
+            range index of the first sample of the alternate grid
+        )"
+    );
+}

python/extensions/pybind_isce3/image/Flatten.h

@@ -0,0 +1,4 @@
+#pragma once
+#include <pybind11/pybind11.h>
+
+void addbindings_flatten(pybind11::module&);

python/extensions/pybind_isce3/image/image.cpp

@@ -2,17 +2,22 @@
 
 #include "Resample.h"
 #include "ResampSlc.h"
+#include "Flatten.h"
 
 namespace py = pybind11;
 
 void addsubmodule_image(py::module & m)
 {
     py::module m_image = m.def_submodule("image");
     py::module m_image_v2 = m_image.def_submodule("v2");
+    py::module m_image_flatten = m_image.def_submodule("flatten");
 
     // Add the resample v2 functionality to the v2 module.
     addbindings_resamp(m_image_v2); 
 
+    // Add the flattening functionality to the image.flatten module
+    addbindings_flatten(m_image_flatten);
+
     // forward declare bound classes for v1
     py::class_<isce3::image::ResampSlc> pyResampSlc(m_image, "ResampSlc");