diff --git a/cxx/isce3/Headers.cmake b/cxx/isce3/Headers.cmake index 5d2543079..e1d21bf6d 100644 --- a/cxx/isce3/Headers.cmake +++ b/cxx/isce3/Headers.cmake @@ -114,6 +114,7 @@ geometry/metadataCubes.h geogrid/getRadarGrid.h geogrid/relocateRaster.h image/forward.h +image/Modulate.h image/Resample.h image/ResampSlc.h image/ResampSlc.icc diff --git a/cxx/isce3/Sources.cmake b/cxx/isce3/Sources.cmake index 171dbe5c2..30005a8cc 100644 --- a/cxx/isce3/Sources.cmake +++ b/cxx/isce3/Sources.cmake @@ -59,6 +59,7 @@ geometry/TopoLayers.cpp geometry/metadataCubes.cpp geogrid/getRadarGrid.cpp geogrid/relocateRaster.cpp +image/Modulate.cpp image/Resample.cpp image/ResampSlc.cpp io/gdal/Dataset.cpp diff --git a/cxx/isce3/core/Poly2d.h b/cxx/isce3/core/Poly2d.h index 5d96c3875..539b1e540 100644 --- a/cxx/isce3/core/Poly2d.h +++ b/cxx/isce3/core/Poly2d.h @@ -72,6 +72,17 @@ class isce3::core::Poly2d { /**Get coefficient by indices*/ inline double getCoeff(int row, int col) const; + /** + * Check if point resides in domain of Poly2d. Added for interoperability with + * LUT2d, but always assumed to be true. + */ + inline bool contains(double y, double x) const { + if (std::isnan(y) || std::isnan(x)){ + return false; + } + return true; + } + /**Evaluate polynomial at given y/azimuth/row ,x/range/col*/ double eval(double y, double x) const; diff --git a/cxx/isce3/image/Modulate.cpp b/cxx/isce3/image/Modulate.cpp new file mode 100644 index 000000000..2743612d0 --- /dev/null +++ b/cxx/isce3/image/Modulate.cpp @@ -0,0 +1,250 @@ +#include "Modulate.h" + +#include +#include +#include + +namespace isce3::image::modulate { + + +template +std::complex _getPixelCarrierPhase( + const double azimuth, + const double range, + const AzRgFunc& carrier_phase, + bool conjugate = false +) { + // Evaluate the pixel's carrier phase + // unit: phase (radians) + const double phase = carrier_phase.eval(azimuth, range); + + // Convert the carrier phase into a unit phasor (i.e. an angle on the unit + // circle in the complex plane). + // unit: unitless (complex) + if(conjugate){ + return std::complex(std::cos(phase), -std::sin(phase)); + } else { + return std::complex(std::cos(phase), std::sin(phase)); + } +} + + +template +void getCarrierPhase( + ArrayRef2D> out, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const bool conjugate, + const std::complex fill_value +) +{ + // unit: azimuth row indices (int) + const size_t phase_block_length = out.rows(); + // unit: range column indices (int) + const size_t phase_block_width = out.cols(); + + // remove carrier from radar data +#pragma omp parallel for collapse(2) + for (size_t az_index = 0; az_index < phase_block_length; ++az_index) { + for (size_t rg_index = 0; rg_index < phase_block_width; ++rg_index) { + // unit: time (seconds) + const double azimuth = + radar_grid.sensingStart() + az_index / radar_grid.prf(); + + // unit: distance (meters) + const double range = radar_grid.startingRange() + + rg_index * radar_grid.rangePixelSpacing(); + + if(not carrier_phase.contains(azimuth, range)){ + out(az_index, rg_index) = fill_value; + continue; + } + + // Get the carrier unit phasor for this pixel + const auto phasor = _getPixelCarrierPhase( + azimuth, range, carrier_phase, conjugate + ); + + // Write the phasor into the output data block + out(az_index, rg_index) = phasor; + } + } // end multithreaded block +} + + +template +void _getCarrierPhaseAtCoords( + ArrayRef2D> out, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const ConstArrayRef2D azimuth_indices, + const ConstArrayRef2D range_indices, + const bool conjugate, + const std::complex fill_value +) +{ + const size_t outWidth = out.cols(); + const size_t outLength = out.rows(); + +#pragma omp parallel for collapse(2) + for (size_t az_index = 0; az_index < outLength; ++az_index){ + for (size_t rg_index = 0; rg_index < outWidth; ++rg_index){ + + // Get the indices on the output grid corresponding to these indices + // on the input grid overall. + const double az_carrier_index = azimuth_indices(az_index, rg_index); + const double rg_carrier_index = range_indices(az_index, rg_index); + + // Azimuth time at the current output pixel + const double azimuth = radar_grid.sensingStart() + az_carrier_index / + radar_grid.prf(); + + // Slant Range at the current output pixel + const double range = radar_grid.startingRange() + rg_carrier_index * + radar_grid.rangePixelSpacing(); + + if(not carrier_phase.contains(azimuth, range)){ + out(az_index, rg_index) = fill_value; + continue; + } + + // Get the carrier phasor for this pixel + const auto phasor = _getPixelCarrierPhase( + azimuth, range, carrier_phase, conjugate + ); + + // Write the phasor into the output data block + out(az_index, rg_index) = phasor; + } + } // end multithreaded block +} + + +template +void modulateCarrierPhase( + ArrayRef2D> slc_data_block, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const bool conjugate, + const std::complex fill_value +) +{ + // unit: azimuth row indices (int) + const size_t phase_block_length = slc_data_block.rows(); + // unit: range column indices (int) + const size_t phase_block_width = slc_data_block.cols(); + + // remove carrier from radar data +#pragma omp parallel for collapse(2) + for ( size_t az_index = 0; az_index < phase_block_length; ++az_index) { + for (size_t rg_index = 0; rg_index < phase_block_width; ++rg_index) { + // unit: time (seconds) + const double azimuth = + radar_grid.sensingStart() + az_index / radar_grid.prf(); + + // unit: distance (meters) + const double range = radar_grid.startingRange() + + rg_index * radar_grid.rangePixelSpacing(); + + if(not carrier_phase.contains(azimuth, range)){ + slc_data_block(az_index, rg_index) = fill_value; + continue; + } + + // Get the carrier phasor for this pixel + const auto phasor = _getPixelCarrierPhase( + azimuth, range, carrier_phase, conjugate + ); + + // Modulate the phasor into the output data block + slc_data_block(az_index, rg_index) *= phasor; + } + } // end multithreaded block +} + + +template +void _modulateCarrierPhaseAtCoords( + ArrayRef2D> slc_data_block, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const ConstArrayRef2D azimuth_indices, + const ConstArrayRef2D range_indices, + const bool conjugate, + const std::complex fill_value +) +{ + const size_t outWidth = slc_data_block.cols(); + const size_t outLength = slc_data_block.rows(); + +#pragma omp parallel for collapse(2) + for (size_t az_index = 0; az_index < outLength; ++az_index){ + for (size_t rg_index = 0; rg_index < outWidth; ++rg_index){ + + // Get the indices on the output grid corresponding to these indices + // on the input grid overall. + const double az_carrier_index = azimuth_indices(az_index, rg_index); + const double rg_carrier_index = range_indices(az_index, rg_index); + + // Azimuth time at the current output pixel + const double azimuth = radar_grid.sensingStart() + az_carrier_index / + radar_grid.prf(); + + // Slant Range at the current output pixel + const double range = radar_grid.startingRange() + rg_carrier_index * + radar_grid.rangePixelSpacing(); + + if(not carrier_phase.contains(azimuth, range)){ + slc_data_block(az_index, rg_index) = fill_value; + continue; + } + + // Get the carrier phasor for this pixel + const auto phasor = _getPixelCarrierPhase( + azimuth, range, carrier_phase, conjugate + ); + + // Modulate the phasor into the output data block + slc_data_block(az_index, rg_index) *= phasor; + } + } // end multithreaded block +} + + +#define EXPLICIT_INSTANTIATION(AzRgFunc) \ +template void getCarrierPhase( \ + ArrayRef2D> out, \ + const AzRgFunc& carrier_phase, \ + const isce3::product::RadarGridParameters& radar_grid, \ + const bool conjugate, \ + const std::complex fill_value \ +); \ +template void _getCarrierPhaseAtCoords( \ + ArrayRef2D> out, \ + const AzRgFunc& carrier_phase, \ + const isce3::product::RadarGridParameters& radar_grid, \ + const ConstArrayRef2D azimuth_indices, \ + const ConstArrayRef2D range_indices, \ + const bool conjugate, \ + const std::complex fill_value \ +); \ +template void modulateCarrierPhase( \ + ArrayRef2D> slc_data_block, \ + const AzRgFunc& carrier_phase, \ + const isce3::product::RadarGridParameters& radar_grid, \ + const bool conjugate, \ + const std::complex fill_value \ +); \ +template void _modulateCarrierPhaseAtCoords( \ + ArrayRef2D> slc_data_block, \ + const AzRgFunc& carrier_phase, \ + const isce3::product::RadarGridParameters& radar_grid, \ + const ConstArrayRef2D azimuth_indices, \ + const ConstArrayRef2D range_indices, \ + const bool conjugate, \ + const std::complex fill_value \ +) +EXPLICIT_INSTANTIATION(isce3::core::LUT2d); +EXPLICIT_INSTANTIATION(isce3::core::Poly2d); + +} // end namespace isce3::image::modulate diff --git a/cxx/isce3/image/Modulate.h b/cxx/isce3/image/Modulate.h new file mode 100644 index 000000000..a7cc80835 --- /dev/null +++ b/cxx/isce3/image/Modulate.h @@ -0,0 +1,155 @@ +#include +#include +#include + +#include +#include +#include + +namespace isce3::image::modulate { + +template +using Array2D = Eigen::Array; + +template +using ArrayRef2D = Eigen::Ref>; + +template +using ConstArrayRef2D = Eigen::Ref>; + +/** + * Acquire the phase of the given carrier of a radar scene. + * + * @param[out] out + * Block of data to be written to. All data in block will be overwritten. + * unit: array2D of complex + * @tparam[in] carrier_phase + * azimuth carrier phase of the SLC data, ::modulate, as a function of azimuth and range. + * @param[in] radar_grid + * parameters for the given radar grid corresponding to the output block. + * @param[in] conjugate + * if true, get the conjugate of the phase. + * @param[in] fill_value + * The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined here as + * any pixels that cannot be evaluated by the carrier_phase function. Poly2d functions + * do not have out-of-bounds pixels, but LUT2d functions may. Defaults to NaN + j*NaN. + */ +template +void getCarrierPhase( + ArrayRef2D> out, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const bool conjugate, + const std::complex fill_value +); + + +/** + * Acquire the phase of the given carrier at each given index of a radar scene. + * + * Note that this function does not perform size checks between `out` and the indices + * blocks and radar grid. + * + * @param[out] out + * Block of data to be written to. All data in block will be overwritten. + * unit: array2D of complex + * @tparam[in] carrier_phase + * azimuth carrier phase of the SLC data, in radians, as a function of azimuth and range. + * @param[in] radar_grid + * parameters for the given radar grid corresponding to the output block. + * @param[in] azimuth_indices + * azimuth index of each pixel in the output block w.r.t the radar grid. Must be the + * same shape as `out`. + * unit: azimuth row indices (int) + * @param[in] range_indices + * range index of each pixel in the output block w.r.t the radar grid. Must be the + * same shape as `out`. + * unit: range column indices (int) + * @param[in] conjugate + * if true, modulate the conjugate of the phase. + * @param[in] fill_value + * The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined here as + * any pixels that cannot be evaluated by the carrier_phase function. Poly2d functions + * do not have out-of-bounds pixels, but LUT2d functions may. Defaults to NaN + j*NaN. + */ +template +void _getCarrierPhaseAtCoords( + ArrayRef2D> out, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const ConstArrayRef2D azimuth_indices, + const ConstArrayRef2D range_indices, + const bool conjugate, + const std::complex fill_value +); + + +/** + * Remove range and azimuth phase carrier from a block of input radar SLC data + * + * @param[out] slc_data_block + * the block of data SLC to be modulated. + * unit: array2D of complex + * @tparam[in] carrier_phase + * carrier phase of the SLC data, in radians, as a function of azimuth and range. + * This phase will be modulated to or demodulated from the image. + * @param[in] radar_grid + * parameters for the given radar grid corresponding to `slc_data_block`. + * @param[in] conjugate + * if true, modulate the conjugate of the phase. + * @param[in] fill_value + * The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined here as + * any pixels that cannot be evaluated by the carrier_phase function. Poly2d functions + * do not have out-of-bounds pixels, but LUT2d functions may. Defaults to NaN + j*NaN. + */ +template +void modulateCarrierPhase( + ArrayRef2D> slc_data_block, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const bool conjugate, + const std::complex fill_value +); + + +/** + * Evaluate and modulate or demodulate the phase carrier onto the given SLC data block. + * + * Note that this function does not perform size checks between `slc_data_block` and the + * indices blocks and radar grid. + * + * @param[out] slc_data_block + * the block of data SLC to be modulated. + * unit: array2D of complex + * @tparam[in] carrier_phase + * carrier phase of the SLC data, in radians, as a function of azimuth and range. + * This phase will be modulated to or demodulated from the image. + * @param[in] radar_grid + * parameters for the given radar grid corresponding to `slc_data_block`. + * @param[in] azimuth_indices + * azimuth index of each pixel in the output block w.r.t the radar grid. Must be the + * same shape as `slc_data_block`. + * unit: azimuth row indices (int) + * @param[in] range_indices + * range index of each pixel in the output block w.r.t the radar grid. Must be the + * same shape as `slc_data_block`. + * unit: range column indices (int) + * @param[in] conjugate + * if true, modulate the conjugate of the phase. + * @param[in] fill_value + * The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined here as + * any pixels that cannot be evaluated by the carrier_phase function. Poly2d functions + * do not have out-of-bounds pixels, but LUT2d functions may. Defaults to NaN + j*NaN. + */ +template +void _modulateCarrierPhaseAtCoords( + ArrayRef2D> slc_data_block, + const AzRgFunc& carrier_phase, + const isce3::product::RadarGridParameters& radar_grid, + const ConstArrayRef2D azimuth_indices, + const ConstArrayRef2D range_indices, + const bool conjugate, + const std::complex fill_value +); + +} // namespace isce3::image::modulate diff --git a/cxx/isce3/image/Resample.cpp b/cxx/isce3/image/Resample.cpp index 0beacc7de..577636e4a 100644 --- a/cxx/isce3/image/Resample.cpp +++ b/cxx/isce3/image/Resample.cpp @@ -3,10 +3,10 @@ #include #include #include +#include namespace isce3::image::v2 { - void resampleToCoords( ArrayRef2D> resampled_data_block, const ConstArrayRef2D> input_data_block, diff --git a/cxx/isce3/image/Resample.h b/cxx/isce3/image/Resample.h index f05ff5fee..05647b132 100644 --- a/cxx/isce3/image/Resample.h +++ b/cxx/isce3/image/Resample.h @@ -27,7 +27,7 @@ using ConstArrayRef2D = Eigen::Ref>; * @param[in] azimuth_input_indices * azimuth (radar-coordinates y) index of the pixels in input grid * @param[in] radar_grid - * RadarGridParameters of radar data + * parameters for the given radar grid * @param[in] native_doppler_lut * native doppler of SLC image * @param[in] fill_value diff --git a/python/extensions/pybind_isce3/Sources.cmake b/python/extensions/pybind_isce3/Sources.cmake index 6732d3dcc..5b23c924b 100644 --- a/python/extensions/pybind_isce3/Sources.cmake +++ b/python/extensions/pybind_isce3/Sources.cmake @@ -59,6 +59,7 @@ geogrid/relocateRaster.cpp geogrid/geogrid.cpp geometry/lookIncFromSr.cpp image/image.cpp +image/Modulate.cpp image/Resample.cpp image/ResampSlc.cpp io/gdal/Dataset.cpp diff --git a/python/extensions/pybind_isce3/image/Modulate.cpp b/python/extensions/pybind_isce3/image/Modulate.cpp new file mode 100644 index 000000000..92f789d43 --- /dev/null +++ b/python/extensions/pybind_isce3/image/Modulate.cpp @@ -0,0 +1,198 @@ +#include "Modulate.h" + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace py = pybind11; + + +template +void addbindings_modulate(py::module & m) +{ + m.def( + "_get_carrier_phase", + py::overload_cast< + isce3::image::modulate::ArrayRef2D>, + const AzRgFunc&, + const isce3::product::RadarGridParameters&, + const bool, + const std::complex + >(&isce3::image::modulate::getCarrierPhase), + py::arg("out"), + py::arg("carrier_phase"), + py::arg("radar_grid"), + py::arg("conjugate"), + py::arg("fill_value") = + std::complex(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()), + R"( + Acquire the phase of the given carrier of a radar scene. + + Parameters + ---------- + out: numpy.ndarray (complex64) + The output phase array to modify. Anything in this array will be + overwritten. + carrier_phase: isce3.core.LUT2d or isce3.core.Poly2d + Carrier phase, in radian, as a function of azimuth and range. + radar_grid: isce3.product.RadarGridParameters + Parameters for the given radar grid. + conjugate: bool + If True, get the conjugate of the phase. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are + defined here as any pixels that cannot be evaluated by the carrier_phase + function. Poly2d functions do not have out-of-bounds pixels, but LUT2d + functions may. Defaults to NaN + j*NaN. + )" + ); + + + m.def( + "_get_carrier_phase_at_coords", + py::overload_cast< + isce3::image::modulate::ArrayRef2D>, + const AzRgFunc&, + const isce3::product::RadarGridParameters&, + const isce3::image::modulate::ConstArrayRef2D, + const isce3::image::modulate::ConstArrayRef2D, + const bool, + const std::complex + >(&isce3::image::modulate::_getCarrierPhaseAtCoords), + py::arg("out"), + py::arg("carrier_phase"), + py::arg("radar_grid"), + py::arg("azimuth_indices"), + py::arg("range_indices"), + py::arg("conjugate"), + py::arg("fill_value") = + std::complex(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()), + R"( + Acquire the phase of the given carrier at each given index of a radar scene. + + Parameters + ---------- + out: numpy.ndarray (complex64) + The output phase array to modify. Anything in this array will be + overwritten. + carrier_phase: isce3.core.LUT2d or isce3.core.Poly2d + Carrier phase, in radian, as a function of azimuth and range. + radar_grid: isce3.product.RadarGridParameters + Parameters for the given radar grid. + azimuth_indices: numpy.ndarray (float64) + Azimuth index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + range_indices: numpy.ndarray (float64) + Range index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + conjugate: bool + If True, get the conjugate of the phase. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are + defined here as any pixels that cannot be evaluated by the carrier_phase + function. Poly2d functions do not have out-of-bounds pixels, but LUT2d + functions may. Defaults to NaN + j*NaN. + )" + ); + + + m.def( + "_modulate_carrier_phase", + py::overload_cast< + isce3::image::modulate::ArrayRef2D>, + const AzRgFunc&, + const isce3::product::RadarGridParameters&, + const bool, + const std::complex + >(&isce3::image::modulate::modulateCarrierPhase), + py::arg("slc_data_block"), + py::arg("carrier_phase"), + py::arg("radar_grid"), + py::arg("conjugate"), + py::arg("fill_value") = + std::complex(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()), + R"( + Evaluate and modulate or demodulate the phase carrier onto the given SLC data + block. + + Parameters + ---------- + slc_data_block: numpy.ndarray (complex64) + The block of SLC data to modulate. + carrier_phase: isce3.core.LUT2d or isce3.core.Poly2d + Carrier phase, in radian, as a function of azimuth and range. This phase + will be modulated to or demodulated from the image. + radar_grid: isce3.product.RadarGridParameters + Parameters for the given radar grid. + conjugate: bool, optional + If True, modulate the conjugate of the phase. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are + defined here as any pixels that cannot be evaluated by the carrier_phase + function. Poly2d functions do not have out-of-bounds pixels, but LUT2d + functions may. Defaults to NaN + j*NaN. + )" + ); + + + m.def( + "_modulate_carrier_phase_at_coords", + py::overload_cast< + isce3::image::modulate::ArrayRef2D>, + const AzRgFunc&, + const isce3::product::RadarGridParameters&, + const isce3::image::modulate::ConstArrayRef2D, + const isce3::image::modulate::ConstArrayRef2D, + const bool, + const std::complex + >(&isce3::image::modulate::_modulateCarrierPhaseAtCoords), + py::arg("slc_data_block"), + py::arg("carrier_phase"), + py::arg("radar_grid"), + py::arg("azimuth_indices"), + py::arg("range_indices"), + py::arg("conjugate"), + py::arg("fill_value") = + std::complex(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()), + R"( + Evaluate and modulate or demodulate the phase carrier onto the given SLC data + block at the given indices. + + Parameters + ---------- + slc_data_block: numpy.ndarray (complex64) + The output phase array to modulate. + carrier_phase: isce3.core.LUT2d or isce3.core.Poly2d + Carrier phase, in radian, as a function of azimuth and range. This phase + will be modulated to or demodulated from the image. + radar_grid: isce3.product.RadarGridParameters + Parameters for the given radar grid. + azimuth_indices: numpy.ndarray (float64) + Azimuth index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + range_indices: numpy.ndarray (float64) + Range index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + conjugate: bool, optional + If True, modulate the conjugate of the phase. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are + defined here as any pixels that cannot be evaluated by the carrier_phase + function. Poly2d functions do not have out-of-bounds pixels, but LUT2d + functions may. Defaults to NaN + j*NaN. + )" + ); +} + +template void addbindings_modulate>(py::module & m); +template void addbindings_modulate(py::module & m); diff --git a/python/extensions/pybind_isce3/image/Modulate.h b/python/extensions/pybind_isce3/image/Modulate.h new file mode 100644 index 000000000..4555ad681 --- /dev/null +++ b/python/extensions/pybind_isce3/image/Modulate.h @@ -0,0 +1,5 @@ +#pragma once +#include + +template +void addbindings_modulate(pybind11::module&); diff --git a/python/extensions/pybind_isce3/image/Resample.cpp b/python/extensions/pybind_isce3/image/Resample.cpp index edd2f9bfe..b915c34c6 100644 --- a/python/extensions/pybind_isce3/image/Resample.cpp +++ b/python/extensions/pybind_isce3/image/Resample.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include diff --git a/python/extensions/pybind_isce3/image/image.cpp b/python/extensions/pybind_isce3/image/image.cpp index c71d72910..339d831ca 100644 --- a/python/extensions/pybind_isce3/image/image.cpp +++ b/python/extensions/pybind_isce3/image/image.cpp @@ -1,5 +1,6 @@ #include "image.h" +#include "Modulate.h" #include "Resample.h" #include "ResampSlc.h" @@ -9,9 +10,14 @@ 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_modulate = m_image.def_submodule("modulate"); // Add the resample v2 functionality to the v2 module. - addbindings_resamp(m_image_v2); + addbindings_resamp(m_image_v2); + + // Add modulation functionality. + addbindings_modulate>(m_image_modulate); + addbindings_modulate(m_image_modulate); // forward declare bound classes for v1 py::class_ pyResampSlc(m_image, "ResampSlc"); diff --git a/python/packages/isce3/image/modulate.py b/python/packages/isce3/image/modulate.py new file mode 100644 index 000000000..3dbce6ed4 --- /dev/null +++ b/python/packages/isce3/image/modulate.py @@ -0,0 +1,255 @@ +from __future__ import annotations + +import journal +import numpy as np + +from isce3.core import LUT2d +from isce3.core.poly2d import Poly2d +from isce3.ext.isce3.image.modulate import ( + _get_carrier_phase, + _get_carrier_phase_at_coords, + _modulate_carrier_phase, + _modulate_carrier_phase_at_coords, +) +from isce3.product import RadarGridParameters + + +def modulate_carrier_phase( + slc_data_block: np.ndarray[np.complex64], + carrier_phase: LUT2d | Poly2d, + radar_grid: RadarGridParameters, + conjugate: bool = False, + fill_value: np.complex64 = np.nan + 1.j * np.nan, +) -> np.ndarray[np.complex64]: + """ + Evaluate and modulate or demodulate the phase carrier onto the given SLC data block. + + Parameters + ---------- + slc_data_block : np.ndarray of np.complex64 + The block of SLC data to modulate. + carrier_phase : LUT2d or Poly2d + Carrier phase of the SLC data, in radian, as a function of azimuth and range. + This phase will be modulated to or demodulated from the image. + radar_grid : RadarGridParameters + Parameters for the given radar grid corresponding to `slc_data_block`. + conjugate : bool, optional + If True, modulate the conjugate of the phase. Defaults to False. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined + here as any pixels that cannot be evaluated by the carrier_phase function. + Poly2d functions do not have out-of-bounds pixels, but LUT2d functions may. + Defaults to NaN + j*NaN. + + Returns + ------- + np.ndarray of np.complex64 + The modulated SLC block. + """ + out_arr = slc_data_block.copy() + out_arr = np.require(out_arr, dtype=np.complex64, requirements=["C", "W"]) + + _modulate_carrier_phase( + slc_data_block=out_arr, + carrier_phase=carrier_phase, + radar_grid=radar_grid, + conjugate=conjugate, + fill_value=fill_value, + ) + + return out_arr + + +def modulate_carrier_phase_at_coords( + slc_data_block: np.ndarray[np.complex64], + carrier_phase: LUT2d | Poly2d, + radar_grid: RadarGridParameters, + azimuth_indices: np.ndarray[np.float64], + range_indices: np.ndarray[np.float64], + conjugate: bool = False, + fill_value: np.complex64 = np.nan + 1.j * np.nan, +) -> np.ndarray[np.complex64]: + """ + Evaluate and modulate or demodulate the phase carrier onto the given SLC data block + at the given indices. + + Parameters + ---------- + slc_data_block : np.ndarray of np.complex64 + The block of SLC data to modulate. + carrier_phase : LUT2d or Poly2d + Carrier phase of the SLC data, in radian, as a function of azimuth and range. + This phase will be modulated to or demodulated from the image. + radar_grid : RadarGridParameters + Parameters for the given radar grid corresponding to `slc_data_block`. + azimuth_indices : np.ndarray of np.float64 + Azimuth index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + range_indices : np.ndarray of np.float64 + Range index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + conjugate : bool, optional + If True, modulate the conjugate of the phase. Defaults to False. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined + here as any pixels that cannot be evaluated by the carrier_phase function. + Poly2d functions do not have out-of-bounds pixels, but LUT2d functions may. + Defaults to NaN + j*NaN. + + Returns + ------- + np.ndarray of np.complex64 + The modulated SLC block. + """ + error_channel = journal.error("modulate.modulate_carrier_phase_at_coords") + + if azimuth_indices.shape != range_indices.shape: + err_log = ( + f"Azimuth indices block shape {azimuth_indices.shape} and range indices " + f"block shape {range_indices.shape} are unequal." + ) + error_channel.log(err_log) + raise ValueError(err_log) + + if azimuth_indices.shape != slc_data_block.shape: + err_log = ( + f"Indices block shapes {azimuth_indices.shape} and SLC data block shape " + f"{slc_data_block.shape} are unequal." + ) + error_channel.log(err_log) + raise ValueError(err_log) + + out_arr = slc_data_block.copy() + out_arr = np.require(out_arr, dtype=np.complex64, requirements=["C", "W"]) + + # Ensure that all of the index data blocks meet the requirements of the + # _modulate_carrier_phase_at_coords pybind (correct dtype, with flag C_CONTIGUOUS) + # These function calls will return conforming copies of the data blocks if they + # are not already conforming. + range_indices = np.require(range_indices, dtype=np.float64, requirements=["C"]) + azimuth_indices = np.require(azimuth_indices, dtype=np.float64, requirements=["C"]) + + _modulate_carrier_phase_at_coords( + slc_data_block=out_arr, + carrier_phase=carrier_phase, + radar_grid=radar_grid, + azimuth_indices=azimuth_indices, + range_indices=range_indices, + conjugate=conjugate, + fill_value=fill_value, + ) + + return out_arr + + +def get_carrier_phase( + carrier_phase: LUT2d | Poly2d, + radar_grid: RadarGridParameters, + conjugate: bool = False, + fill_value: np.complex64 = np.nan + 1.j * np.nan, +) -> np.ndarray[np.complex64]: + """ + Acquire the phase of the given carrier of a radar scene. + + Parameters + ---------- + carrier_phase : LUT2d or Poly2d + Carrier phase, in radian, as a function of azimuth and range. + radar_grid : RadarGridParameters + Parameters for the given radar grid corresponding to the output block. + conjugate : bool, optional + If True, get the conjugate of the phase. Defaults to False. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined + here as any pixels that cannot be evaluated by the carrier_phase function. + Poly2d functions do not have out-of-bounds pixels, but LUT2d functions may. + Defaults to NaN + j*NaN. + + Returns + ------- + np.ndarray of np.complex64 + The carrier phase, in the form of complex unit vectors. + """ + out_arr = np.full( + (radar_grid.length, radar_grid.width), + fill_value=fill_value, + dtype=np.complex64, + ) + + _get_carrier_phase( + out=out_arr, + carrier_phase=carrier_phase, + radar_grid=radar_grid, + conjugate=conjugate, + fill_value=fill_value, + ) + + return out_arr + + +def get_carrier_phase_at_coords( + carrier_phase: LUT2d | Poly2d, + radar_grid: RadarGridParameters, + azimuth_indices: np.ndarray[np.float64], + range_indices: np.ndarray[np.float64], + conjugate: bool = False, + fill_value: np.complex64 = np.nan + 1.j * np.nan, +) -> np.ndarray[np.complex64]: + """ + Acquire the phase of the given carrier at each given index of a radar scene. + + Parameters + ---------- + carrier_phase : LUT2d or Poly2d + Carrier phase, in radian, as a function of azimuth and range. + radar_grid : RadarGridParameters + Parameters for the radar grid corresponding to the output block. + azimuth_indices : np.ndarray of np.float64 + Azimuth index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + range_indices : np.ndarray of np.float64 + Range index of each output coordinate pixel in the given radar coordinate + system. Must be the same shape as phase_data_block. + conjugate : bool, optional + If True, get the conjugate of the phase. Defaults to False. + fill_value: complex + The value to fill out-of-bounds pixels with. Out-of-bounds pixels are defined + here as any pixels that cannot be evaluated by the carrier_phase function. + Poly2d functions do not have out-of-bounds pixels, but LUT2d functions may. + Defaults to NaN + j*NaN. + + Returns + ------- + np.ndarray of np.complex64 + The carrier phase, in the form of complex unit vectors. + """ + error_channel = journal.error("modulate.get_carrier_phase_at_coords") + + if azimuth_indices.shape != range_indices.shape: + err_log = ( + f"Azimuth indices block shape {azimuth_indices.shape} and range indices " + f"block shape {range_indices.shape} are unequal." + ) + error_channel.log(err_log) + raise ValueError(err_log) + + out_arr = np.full(azimuth_indices.shape, fill_value=fill_value, dtype=np.complex64) + + # Ensure that all of the index data blocks meet the requirements of the + # _get_carrier_phase_at_coords pybind (correct dtype, with flag C_CONTIGUOUS) + # These function calls will return conforming copies of the data blocks if they + # are not already conforming. + range_indices = np.require(range_indices, dtype=np.float64, requirements=["C"]) + azimuth_indices = np.require(azimuth_indices, dtype=np.float64, requirements=["C"]) + + _get_carrier_phase_at_coords( + out=out_arr, + carrier_phase=carrier_phase, + radar_grid=radar_grid, + azimuth_indices=azimuth_indices, + range_indices=range_indices, + conjugate=conjugate, + fill_value=fill_value, + ) + + return out_arr diff --git a/python/packages/isce3/image/v2/resample_slc.py b/python/packages/isce3/image/v2/resample_slc.py index dd3159ce1..8e4be4ad0 100644 --- a/python/packages/isce3/image/v2/resample_slc.py +++ b/python/packages/isce3/image/v2/resample_slc.py @@ -1,13 +1,18 @@ from __future__ import annotations -from collections.abc import Sequence +from collections.abc import Iterable, Sequence from time import perf_counter import journal import numpy as np from isce3.core import SINC_HALF, LUT2d +from isce3.core.poly2d import Poly2d from isce3.core.resample_block_generators import get_blocks, get_blocks_by_offsets from isce3.ext.isce3.image.v2 import _resample_to_coords +from isce3.image.modulate import ( + modulate_carrier_phase, + modulate_carrier_phase_at_coords, +) from isce3.io.dataset import DatasetReader, DatasetWriter from isce3.product import RadarGridParameters @@ -24,6 +29,9 @@ def resample_slc_blocks( quiet: bool = False, fill_value: np.complex64 = np.nan + 1.0j * np.nan, with_gpu: bool = False, + *, + phase_carriers: Iterable[LUT2d | Poly2d] | None = None, + remodulate: bool = False, ) -> None: """ Resamples one or more SLCs onto a geometry described by given offsets datasets. @@ -63,6 +71,12 @@ def resample_slc_blocks( with_gpu : bool, optional If True, run the GPU resample workflow. If False, run the CPU resample workflow. Defaults to False. + phase_carriers: Iterable of LUT2d or Poly2d or None, optional + Carrier phase, in radians, to remove prior to resampling. If None, the carrier + phase is assumed to be zero. Defaults to None. + remodulate: bool, optional + If True, remodulate the phase_carriers phase into the output data. Use only if + phase_carriers is also given. Defaults to False. """ info_channel = journal.info("resample_slc.resample_slc_blocks") warning_channel = journal.warning("resample_slc.resample_slc_blocks") @@ -81,6 +95,11 @@ def resample_slc_blocks( error_channel.log(err_log) raise ValueError(err_log) + if remodulate and (phase_carriers is None): + err_log = "If remodulate is True, phase_carriers must also be given." + error_channel.log(err_log) + raise ValueError(err_log) + for in_dataset in input_slcs: in_dataset_dtype = in_dataset.dtype if in_dataset_dtype not in [np.complex64, np.complex128]: @@ -212,6 +231,22 @@ def resample_slc_blocks( # Run the resampling algorithm on the given blocks. for i in range(len(input_blocks)): input_block = input_blocks[i] + block_grid = input_radar_grid[in_slices] + + if phase_carriers is not None: + if not quiet: + info_channel.log( + f"demodulating input SLC for block {out_block_slice}..." + ) + + for carrier in phase_carriers: + input_block = modulate_carrier_phase( + slc_data_block=input_block, + carrier_phase=carrier, + radar_grid=block_grid, + conjugate=True, + ) + if not quiet: info_channel.log( f"interpolating to output SLC for block {out_block_slice}..." @@ -219,15 +254,34 @@ def resample_slc_blocks( # Reporting input block shape for debugging info_channel.log(f"Input block: {in_slices}") - output_blocks[i] = resample( + output_block = resample( input_block, range_index_grid, azimuth_index_grid, - input_radar_grid[in_slices], + block_grid, doppler, fill_value, ) + + if remodulate and (phase_carriers is not None): + if not quiet: + info_channel.log( + f"remodulating output SLC for block {out_block_slice}..." + ) + + for carrier in phase_carriers: + output_block = modulate_carrier_phase_at_coords( + slc_data_block=output_block, + carrier_phase=carrier, + radar_grid=block_grid, + azimuth_indices=azimuth_index_grid, + range_indices=range_index_grid, + conjugate=False, + ) + + output_blocks[i] = output_block + block_processing_timer += perf_counter() # The resampling blocks have now been filled. For each output dataset, write @@ -405,7 +459,7 @@ def resample_to_coords( raise ValueError(err_log) # Ensure that all of the input data blocks meet the requirements of the - # _resample_to_coords pybind (correct dtype, with flags C_CONTIGUOUS and WRITABLE) + # _resample_to_coords pybind (correct dtype, with flag C_CONTIGUOUS) # These function calls will return conforming copies of the data blocks if they # are not already conforming. input_data_block = np.require( diff --git a/tests/python/packages/CMakeLists.txt b/tests/python/packages/CMakeLists.txt index d37588fe1..7b92997f9 100644 --- a/tests/python/packages/CMakeLists.txt +++ b/tests/python/packages/CMakeLists.txt @@ -19,6 +19,7 @@ isce3/geometry/bounding_radar_grid.py isce3/geometry/polygons.py isce3/geometry/radar_grid_spacing.py isce3/image/resample_slc.py +isce3/image/modulate.py isce3/io/gdal/gdal_raster.py isce3/io/background.py isce3/matchtemplate/import_ampcor.py diff --git a/tests/python/packages/isce3/image/modulate.py b/tests/python/packages/isce3/image/modulate.py new file mode 100644 index 000000000..e409438b8 --- /dev/null +++ b/tests/python/packages/isce3/image/modulate.py @@ -0,0 +1,287 @@ +"""Test resample_slc ver. 2""" +from __future__ import annotations + +import pytest +import numpy as np + +from isce3.image.modulate import ( + get_carrier_phase, + get_carrier_phase_at_coords, + modulate_carrier_phase, + modulate_carrier_phase_at_coords, +) +from isce3.core import DateTime, LUT2d +from isce3.product import RadarGridParameters + +from ..resample_slc_utils import validate_test_results + + +def generate_carrier_ramp_complex( + grid_params: RadarGridParameters, + az_indices: np.ndarray, + rg_indices: np.ndarray, + az_frequency: float, + rg_frequency: float, +) -> np.ndarray: + """ + Evaluate the carrier phase ramp at the given set of indices. + + Parameters + ---------- + grid_params : RadarGridParameters + The radar grid parameters object for the doppler ramp. + az_indices : np.ndarray + The indices to evaluate the carrier at. + rg_indices : np.ndarray + The indices to evaluate the carrier at. + az_frequency : float + The azimuth carrier frequency. + rg_frequency : float + The range carrier frequency. + + Returns + ------- + np.ndarray + An array of the carrier phase ramp evaluated at each index passed in. + """ + # Trivially, for zero-doppler, just return an array of 1+0j + if az_frequency == 0 or rg_frequency == 0: + return np.ones(az_indices.shape, dtype=np.complex64) + + # Get the absolute azimuth time at each index. + azimuth = grid_params.sensing_start + az_indices / grid_params.prf + range = grid_params.starting_range + rg_indices * grid_params.range_pixel_spacing + + # Evaluate and return the ramp. + az_ramp = np.exp(1.0j * 2 * np.pi * az_frequency * azimuth) + rg_ramp = np.exp(1.0j * 2 * np.pi * rg_frequency * range) + + return az_ramp * rg_ramp + + +def generate_carrier_lut( + grid_params: RadarGridParameters, + az_frequency: float, + rg_frequency: float, +) -> LUT2d: + """ + Create a constant carrier LUT. + + Parameters + ---------- + grid_params : RadarGridParameters + The radar grid parameters to evaluate azimuth and range with. + az_frequency : float + The azimuth carrier frequency of the LUT to be created. + rg_frequency : float + The range carrier frequency of the LUT to be created. + + Returns + ------- + LUT2d + The generated LUT. + """ + + # Get the dimensions and indices of the array + array_length = grid_params.length + array_width = grid_params.width + az_indices = np.arange(array_length) + rg_indices = np.arange(array_width) + + # Get the azimuth times and range distances + azimuth = grid_params.sensing_start + az_indices / grid_params.prf + range = grid_params.starting_range + rg_indices * grid_params.range_pixel_spacing + + # Get the LUT array by evaluating the array at each index + az_array = np.full( + shape=(array_length, array_width), + fill_value=2 * np.pi * az_frequency, + dtype=np.float64, + ) * azimuth[:, np.newaxis] + rg_array = np.full( + shape=(array_length, array_width), + fill_value=2 * np.pi * rg_frequency, + dtype=np.float64, + ) * range[np.newaxis, :] + lut_array = az_array + rg_array + + # Generate the LUT + return LUT2d(range, azimuth, lut_array) + + +@pytest.mark.parametrize( + "function", + [ + "get_carrier_phase", + "get_carrier_phase_at_coords", + "modulate_carrier_phase", + "modulate_carrier_phase_at_coords" + ] +) +class TestModulate: + """Tests for the image.modulate code.""" + + @pytest.mark.parametrize("az_freq", [0.1, 0.25, 0.5]) + @pytest.mark.parametrize("rg_freq", [0.1, 0.25, 0.5]) + @pytest.mark.parametrize("conjugate", [True, False]) + def test_modulation_constant_frequency( + self, + az_freq: float, + rg_freq: float, + function: str, + conjugate: bool, + ) -> tuple[np.ndarray[np.complex64], np.ndarray[np.complex64]]: + """ + Tests the four carrier phase acquisition/modulation functions. + + Fixtures + -------- + az_freq : float + The azimuth carrier frequency for this test, in Hz. + rg_freq : float + The azimuth carrier frequency for this test, in Hz. + function : str + The name of the function to test. + conjugate : bool + If True, conjugate the output signal; else do not. + """ + # The size of the generated image. + az_length = 100 + rg_width = 100 + out_shape = (az_length, rg_width) + + # Set up a dummy radar grid + radar_grid: RadarGridParameters = RadarGridParameters( + sensing_start=1, + wavelength=1, + prf=1, + starting_range=1, + range_pixel_spacing=1, + look_side="right", + length=az_length, + width=rg_width, + ref_epoch=DateTime(), + ) + + # Generate a doppler centroid that has a ramp + lut: LUT2d = generate_carrier_lut( + grid_params=radar_grid, + az_frequency=az_freq, + rg_frequency=rg_freq, + ) + + # All of the functions call for a radar grid, carrier phase, and conjugate bool. + # Begin putting together a set of keyword arguments, since the function + # signatures are all very similar. + signal = np.full( + (az_length, rg_width), + fill_value=1. + 0.j, + dtype=np.complex64, + ) + + carrier_ramp_complex: np.ndarray[np.complex64] + + kwargs = { + "radar_grid": radar_grid, + "carrier_phase": lut, + "conjugate": conjugate + } + + # Perform the interpolation. + if function in ["get_carrier_phase", "modulate_carrier_phase"]: + az_indices, rg_indices = np.indices((az_length, rg_width), dtype=np.float64) + + # Create the expected output signal. For this test, a simple phase ramp in + # complex phasor format is the output. + carrier_ramp_complex = generate_carrier_ramp_complex( + grid_params=radar_grid, + az_indices=az_indices, + rg_indices=rg_indices, + az_frequency=az_freq, + rg_frequency=rg_freq, + ) + + if function == "get_carrier_phase": + signal = get_carrier_phase(**kwargs) + elif function == "modulate_carrier_phase": + # A dummy SLC of 1 + 0j to modulate - this will give the carrier + # phase. + input_slc = np.full( + (az_length, rg_width), + fill_value=1. + 0.j, + dtype=np.complex64, + ) + kwargs["slc_data_block"] = input_slc + signal = modulate_carrier_phase(**kwargs) + + elif function in [ + "get_carrier_phase_at_coords", + "modulate_carrier_phase_at_coords" + ]: + # Set the offsets at random positions with a range of -1.5 to 1.5 with a + # flat probability distribution. This ensures that the difference in phase + # and potential edge effects near the ends of an image are detectable. + mag_offset = 3 + az_offsets = np.random.random(out_shape) * mag_offset - mag_offset/2 + rg_offsets = np.random.random(out_shape) * mag_offset - mag_offset/2 + + # Add the offsets to these indices to get the indices in the ground truth + # grid. + rows, cols = np.indices(out_shape) + azimuth_indices = np.array(az_offsets + rows, dtype=np.float64) + range_indices = np.array(rg_offsets + cols, dtype=np.float64) + + # Create the expected output signal. For this test, a simple phase ramp in + # complex phasor format is the output. + carrier_ramp_complex = generate_carrier_ramp_complex( + grid_params=radar_grid, + az_indices=azimuth_indices, + rg_indices=range_indices, + az_frequency=az_freq, + rg_frequency=rg_freq, + ) + + kwargs["azimuth_indices"] = azimuth_indices + kwargs["range_indices"] = range_indices + + if function == "get_carrier_phase_at_coords": + signal = get_carrier_phase_at_coords(**kwargs) + elif function == "modulate_carrier_phase_at_coords": + # A dummy SLC of 1 + 0j to modulate - this will give the carrier + # phase. + input_slc = np.full( + (az_length, rg_width), + fill_value=1. + 0.j, + dtype=np.complex64, + ) + kwargs["slc_data_block"] = input_slc + signal = modulate_carrier_phase_at_coords(**kwargs) + + # If conjugate is True, the ground-truth array is currently the conjugate + # of the output array (we hope) and must be conjugated for validation. + if conjugate: + carrier_ramp_complex = np.conjugate(carrier_ramp_complex) + try: + # Validate the generated data against the true data. + # The correlation is expected to be very high and the standard deviation + # very low. The percentage of NaN values will be variable depending on + # offset distance for the "at_coords" functions and zero for the others. + validate_test_results( + test_arr=signal, + true_arr=carrier_ramp_complex, + correlation_min=0.99999, + phase_stdev_max=1e-6, + nan_percent_max=3, + buffer_size_az=2, + buffer_size_rg=2, + az_offset=0, + rg_offset=0, + ) + except AssertionError as err: + # If an error is caught in the validation function, add some clarifying + # notes for readability. + err.add_note(f"function: {function}") + err.add_note(f"az_freq: {az_freq}") + err.add_note(f"rg_freq: {rg_freq}") + err.add_note(f"conjugate: {conjugate}") + raise err diff --git a/tests/python/packages/isce3/resample_slc_utils.py b/tests/python/packages/isce3/resample_slc_utils.py index 7728de2c6..30e9a822b 100644 --- a/tests/python/packages/isce3/resample_slc_utils.py +++ b/tests/python/packages/isce3/resample_slc_utils.py @@ -133,20 +133,20 @@ def validate_test_results( if not corr > correlation_min: fail = True fail_strings.append( - f"Correlation {corr} of resample output is less than minimum acceptable " + f"Correlation {corr} of test output is less than minimum acceptable " f"value {correlation_min}" ) if not phase_stdev_degrees < phase_stdev_max: fail = True fail_strings.append( - f"Phase standard deviation between ground-truth and resample output " + f"Phase standard deviation between ground-truth and test output " f"{phase_stdev_degrees} is greater than maximum acceptable value " f"{phase_stdev_max}" ) if not nan_percent < nan_percent_max: fail = True fail_strings.append( - f"percentage of NaNs in resample output: {nan_percent}% is greater than " + f"percentage of NaNs in test output: {nan_percent}% is greater than " f"maximum acceptable value: {nan_percent_max}%" )