diff --git a/CMakeLists.txt b/CMakeLists.txt index 827dae3cd..030448a26 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -90,6 +90,11 @@ if(ISCE3_WITH_CYTHON) message(ERROR "isce3's cython extension has been removed!") endif() +option(ISCE3_ENABLE_FBP_TIMING "Log factorized backprojection timing" OFF) +if(ISCE3_ENABLE_FBP_TIMING) + add_definitions(-DISCE3_ENABLE_FBP_TIMING) +endif() + ###Layout same install directory structure as pyre include(GNUInstallDirs) InitInstallDirLayout() diff --git a/cxx/isce3/Headers.cmake b/cxx/isce3/Headers.cmake index 5d2543079..26eac57fd 100644 --- a/cxx/isce3/Headers.cmake +++ b/cxx/isce3/Headers.cmake @@ -175,6 +175,7 @@ signal/Looks.h signal/Looks.icc signal/multilook.h signal/NFFT.h +signal/NFFT2d.h signal/shiftSignal.h signal/signalUtils.h signal/Signal.h diff --git a/cxx/isce3/Sources.cmake b/cxx/isce3/Sources.cmake index 171dbe5c2..88cb37559 100644 --- a/cxx/isce3/Sources.cmake +++ b/cxx/isce3/Sources.cmake @@ -103,6 +103,7 @@ signal/Filter.cpp signal/flatten.cpp signal/Looks.cpp signal/NFFT.cpp +signal/NFFT2d.cpp signal/shiftSignal.cpp signal/signalUtils.cpp signal/Signal.cpp diff --git a/cxx/isce3/core/Interp2d.icc b/cxx/isce3/core/Interp2d.icc index 441683b36..09feb2f59 100644 --- a/cxx/isce3/core/Interp2d.icc +++ b/cxx/isce3/core/Interp2d.icc @@ -33,8 +33,13 @@ DataType interp2d(const Kernel& kernelx, // Do X interp at each Y index. for (int i_kernely = 0; i_kernely < widthy; ++i_kernely) { long i_datay = i_kernely + lowy; - if (periodic) - i_datay %= ny; + if (periodic) { + // XXX need both operands signed for correct result + const auto lny = static_cast(ny); + i_datay %= lny; + // Careful that C++ modulo retains sign of dividend. + if (i_datay < 0) i_datay += lny; + } if ((i_datay >= 0) and (i_datay < ny)) { const DataType* zi = &z[i_datay * stridey]; const DataType* px = detail::get_contiguous_view_or_copy( diff --git a/cxx/isce3/core/Kernels.h b/cxx/isce3/core/Kernels.h index 17e1af44c..6fe428cdd 100644 --- a/cxx/isce3/core/Kernels.h +++ b/cxx/isce3/core/Kernels.h @@ -103,6 +103,10 @@ class NFFTKernel : public Kernel { T operator()(double x) const override; + int kernel_radius() const { return _m; } + int data_size() const { return _n; } + int fft_size() const { return _fft_size; } + private: int _m; int _n; diff --git a/cxx/isce3/core/Linspace.h b/cxx/isce3/core/Linspace.h index cc390d21b..07ce68fdf 100644 --- a/cxx/isce3/core/Linspace.h +++ b/cxx/isce3/core/Linspace.h @@ -1,6 +1,7 @@ #pragma once #include "Common.h" +#include namespace isce3 { namespace core { @@ -120,6 +121,17 @@ class Linspace { constexpr int search(U) const; + /** + * Get the boundaries implied by a Linspace of bin center coordinates. + * + * \returns The array [leading_edge, trailing_edge] where leading_edge is a + * half step ahead of first() and trailing_edge is a half step beyond + * last() + */ + CUDA_HOSTDEV + constexpr + std::array bounds() const; + private: T _first = {}; T _spacing = {}; diff --git a/cxx/isce3/core/Linspace.icc b/cxx/isce3/core/Linspace.icc index 172f98432..0ec2d0c54 100644 --- a/cxx/isce3/core/Linspace.icc +++ b/cxx/isce3/core/Linspace.icc @@ -106,4 +106,13 @@ bool operator!=(const Linspace & lhs, const Linspace & rhs) return !(lhs == rhs); } +template +CUDA_HOSTDEV +constexpr +std::array Linspace::bounds() const +{ + const T half = spacing() / 2; + return std::array{first() - half, last() + half}; +} + }} diff --git a/cxx/isce3/core/detail/Interp1d.h b/cxx/isce3/core/detail/Interp1d.h index 6df956bfa..0e9ac0e0a 100644 --- a/cxx/isce3/core/detail/Interp1d.h +++ b/cxx/isce3/core/detail/Interp1d.h @@ -62,8 +62,16 @@ const DataType* get_contiguous_view_or_copy(DataType block[], int width, } // else if (periodic) { + // Careful that C++ modulo retains sign of dividend. + if (low < 0) { + // XXX need both operands signed for correct result + const auto lsize = static_cast(size); + low %= lsize; + low += lsize; + } for (int i = 0; i < width; ++i) { - long j = ((low + i) % size) * stride; + // Already guaranteed positive from above. + auto j = ((low + i) % size) * stride; block[i] = data[j]; } } else { diff --git a/cxx/isce3/cuda/CMakeLists.txt b/cxx/isce3/cuda/CMakeLists.txt index 7198fc0ea..52add3b04 100644 --- a/cxx/isce3/cuda/CMakeLists.txt +++ b/cxx/isce3/cuda/CMakeLists.txt @@ -33,6 +33,7 @@ target_link_libraries(${LISCECUDA} PUBLIC ${LISCE} ${CUDART_LIBRARY} ${CUDAFFT_LIBRARY} + OpenMP::OpenMP_CUDA_Optional ) # Specify API version and build version (used to generate name and soname diff --git a/cxx/isce3/cuda/Headers.cmake b/cxx/isce3/cuda/Headers.cmake index 6c4ea178a..e97d0541d 100644 --- a/cxx/isce3/cuda/Headers.cmake +++ b/cxx/isce3/cuda/Headers.cmake @@ -13,6 +13,8 @@ core/gpuPoly2d.h core/gpuProjections.h core/Interp1d.h core/Interp1d.icc +core/Interp2d.h +core/Interp2d.icc core/Kernels.h core/Kernels.icc core/Orbit.h @@ -48,4 +50,5 @@ signal/gpuCrossMul.h signal/gpuFilter.h signal/gpuLooks.h signal/gpuSignal.h +signal/NFFT2d.h ) diff --git a/cxx/isce3/cuda/Sources.cmake b/cxx/isce3/cuda/Sources.cmake index b5ab0817f..b602afa71 100644 --- a/cxx/isce3/cuda/Sources.cmake +++ b/cxx/isce3/cuda/Sources.cmake @@ -53,4 +53,5 @@ signal/gpuFilter.cu signal/gpuLooks.cu signal/gpuRangeFilter.cu signal/gpuSignal.cu +signal/NFFT2d.cu ) diff --git a/cxx/isce3/cuda/core/Interp2d.h b/cxx/isce3/cuda/core/Interp2d.h new file mode 100644 index 000000000..341911f2e --- /dev/null +++ b/cxx/isce3/cuda/core/Interp2d.h @@ -0,0 +1,37 @@ +#pragma once + +#include "forward.h" +#include + +namespace isce3::cuda::core { + +/** Interpolate Matrix z at point (x,y) + * + * @tparam KernelType kernel element type + * @tparam DataType data element type + * + * @param[in] kernelx Kernel function to use for interpolation in x direction + * @param[in] kernely Kernel function to use for interpolation in y direction + * @param[in] z Matrix to interpolate. + * @param[in] nx Number of x samples. + * @param[in] stridex Stride between x samples. + * @param[in] ny Number of y samples. + * @param[in] stridex Stride between y samples. + * @param[in] x Desired sample (0 <= x < nx). + * @param[in] y Desired sample (0 <= y < ny). + * @param[in] periodic Use periodic boundary condition. Default = false. + * @returns Interpolated value or 0 if kernel would run off array. + * + * Matrix z will be addressed as z[ix * stridex + iy * stridey] for + * 0 <= ix < nx and 0 <= iy < ny. + */ +template +CUDA_HOSTDEV +DataType interp2d(const KernelX& kernelx, + const KernelY& kernely, const DataType* z, size_t nx, + size_t stridex, size_t ny, size_t stridey, double x, double y, + bool periodic = false); + +} + +#include "Interp2d.icc" diff --git a/cxx/isce3/cuda/core/Interp2d.icc b/cxx/isce3/cuda/core/Interp2d.icc new file mode 100644 index 000000000..2debcbeff --- /dev/null +++ b/cxx/isce3/cuda/core/Interp2d.icc @@ -0,0 +1,139 @@ +#include + +#include "Interp1d.h" + +namespace isce3::cuda::core { + +namespace detail { + template + CUDA_HOSTDEV + void interp1d_coeffs(const Kernel& kernel, const double t, + long* low, KT coeffs[]) + { + int width = int(ceil(kernel.width())); + long i0 = 0; + if (width % 2 == 0) { + i0 = static_cast(ceil(t)); + } else { + i0 = static_cast(round(t)); + } + *low = i0 - width / 2; // integer division implicit floor() + for (int i = 0; i < width; ++i) { + double ti = i + (*low) - t; + coeffs[i] = kernel(ti); + } + } + + template + CUDA_HOSTDEV + const DataType* get_contiguous_view_or_copy(DataType block[], int width, + long low, const DataType* data, size_t size, size_t stride, + bool periodic) + { + const long high = low + width; + if ((stride == 1) and (low >= 0) and (high <= size)) { + return &data[low]; + } + // else + if (periodic) { + // Careful that C++ modulo retains sign of dividend. + if (low < 0) { + // XXX need both operands signed for correct result + const auto lsize = static_cast(size); + low %= lsize; + low += lsize; + } + for (int i = 0; i < width; ++i) { + // Already guaranteed positive from above. + auto j = ((low + i) % size) * stride; + block[i] = data[j]; + } + } else { + for (int i = 0; i < width; ++i) { + long j = (low + i) * stride; + if ((j >= 0) and (j < size)) { + block[i] = data[j]; + } else { + block[i] = static_cast(0); + } + } + } + return block; + } + + template + CUDA_HOSTDEV + auto inner_product(const int width, const TX x[], const TY y[]) + { + using namespace isce3::math::complex_operations; + using TO = typename std::common_type::type; + TO sum = 0; + + for (int i = 0; i < width; ++i) { + sum += x[i] * y[i]; + } + return sum; + } +} + +template +CUDA_HOSTDEV +DataType interp2d(const KernelX& kernelx, + const KernelY& kernely, const DataType* z, size_t nx, + size_t stridex, size_t ny, size_t stridey, double x, double y, + bool periodic) +{ + using namespace isce3::math::complex_operations; + + const int MAX_WIDTH = 16; + + // Small-size optimization to avoid heap allocation. The fixed-size stack + // buffers below hold widthx/widthy elements, so an over-wide kernel would + // overflow them. Guard against that (device-side assert traps the kernel; + // host-side aborts) rather than silently corrupting memory. + const int widthx = static_cast(ceil(kernelx.width())); + const int widthy = static_cast(ceil(kernely.width())); + assert(widthx > 0 and widthx <= MAX_WIDTH); + assert(widthy > 0 and widthy <= MAX_WIDTH); + + using TX = typename KernelX::value_type; + using TY = typename KernelY::value_type; + static_assert(std::is_same::value); + TX coeffsx[MAX_WIDTH]; + TY coeffsy[MAX_WIDTH]; + DataType datax[MAX_WIDTH], datay[MAX_WIDTH]; + + // Pre-compute (widthx + widthy) coefficients rather than calculating them + // (widthx * widthy) times inside the loop. + // In principle we could delay calculating coeffsy until after the loop, in + // which case we'd only need a single coefficient buffer. But we need lowy + // which comes from the same API call, and we're probably okay on stack + // space anyhow. + long lowx = 0, lowy = 0; + detail::interp1d_coeffs(kernelx, x, &lowx, coeffsx); + detail::interp1d_coeffs(kernely, y, &lowy, coeffsy); + + // Do X interp at each Y index. + for (int i_kernely = 0; i_kernely < widthy; ++i_kernely) { + long i_datay = i_kernely + lowy; + if (periodic) { + // XXX need both operands signed for correct result + const auto lny = static_cast(ny); + i_datay %= lny; + // Careful that C++ modulo retains sign of dividend. + if (i_datay < 0) i_datay += lny; + } + if ((i_datay >= 0) and (i_datay < ny)) { + const DataType* zi = &z[i_datay * stridey]; + const DataType* px = detail::get_contiguous_view_or_copy( + datax, widthx, lowx, zi, nx, stridex, periodic); + datay[i_kernely] = detail::inner_product(widthx, coeffsx, px); + } else { + datay[i_kernely] = 0; + } + } + // Do Y interp. + return detail::inner_product(widthy, coeffsy, datay); +} + +} diff --git a/cxx/isce3/cuda/core/Kernels.h b/cxx/isce3/cuda/core/Kernels.h index 0e62448c5..dd82a4d62 100644 --- a/cxx/isce3/cuda/core/Kernels.h +++ b/cxx/isce3/cuda/core/Kernels.h @@ -125,6 +125,60 @@ class KnabKernel : public Kernel> { double _bandwidth; }; +/** + * NFFT time-domain kernel + * + * This is called \f$ \phi(x) \f$ in the NFFT papers @cite keiner2009, + * specifically the Kaiser-Bessel window function. + * The domain is scaled so that usage is the same as other ISCE kernels, e.g., + * for x in [0,n) instead of [-0.5,0.5). + */ +template +class NFFTKernel : public Kernel> { + using Base = Kernel>; + friend Base; + +public: + /** A non-owning kernel view type that can be passed to device code */ + using view_type = NFFTKernel; + + /** + * Construct a new NFFTKernel object. + * + * \param[in] m Half kernel size (width = 2*m+1) + * \param[in] n Length of input signal + * \param[in] fft_size FFT transform size (> n) + */ + NFFTKernel(int m, int n, int fft_size); + + /** Construct from corresponding host kernel object */ + NFFTKernel(const isce3::core::NFFTKernel& other); + + /** Get half kernel size. */ + int kernel_radius() const { return m_; } + + /** Get length of input signal. */ + int data_size() const { return n_; } + + /** Get FFT transform size. */ + int fft_size() const { return fft_size_; } + + explicit operator isce3::core::NFFTKernel() const { + return {kernel_radius(), data_size(), fft_size()}; + } + +protected: + /** \internal Implementation of \p operator() */ + CUDA_HOSTDEV T eval(double t) const; + +private: + int m_; + int n_; + int fft_size_; + T scale_; + T b_; +}; + /** A non-owning reference to a TabulatedKernel object */ template class TabulatedKernelView : public Kernel> { diff --git a/cxx/isce3/cuda/core/Kernels.icc b/cxx/isce3/cuda/core/Kernels.icc index 3db5ce99d..fe3fb2eb9 100644 --- a/cxx/isce3/cuda/core/Kernels.icc +++ b/cxx/isce3/cuda/core/Kernels.icc @@ -49,6 +49,54 @@ CUDA_HOSTDEV inline T KnabKernel::eval(double t) const return sinc(x) * detail::samplingWindow(t, Base::halfwidth(), bandwidth()); } +namespace detail { +template +inline void calc_nfft_scale_factors(T& b, T& scale, int m, int n, int fft_size) +{ + b = M_PI * (2.0 - 1.0 * n / fft_size); + scale = 1.0 / (M_PI * isce3::math::bessel_i0(m * b)); +} +} // namespace detail + +template +NFFTKernel::NFFTKernel(int m, int n, int fft_size) + : Base(2 * m + 1), m_(m), n_(n), fft_size_(fft_size) +{ + if ((m < 1) || (n < 1) || (fft_size < 1)) { + throw isce3::except::LengthError( + ISCE_SRCINFO(), "NFFT parameters must be positive."); + } + detail::calc_nfft_scale_factors(b_, scale_, m, n, fft_size); +} + +template +NFFTKernel::NFFTKernel(const isce3::core::NFFTKernel& other) + : Base(2 * other.kernel_radius() + 1), m_(other.kernel_radius()), + n_(other.data_size()), fft_size_(other.fft_size()) +{ + detail::calc_nfft_scale_factors(b_, scale_, other.kernel_radius(), + other.data_size(), other.fft_size()); +} + +template +CUDA_HOSTDEV inline T NFFTKernel::eval(double t) const +{ + T x2 = t * t - m_ * m_; + // x=0 + if (std::abs(x2) < std::numeric_limits::epsilon()) { + return b_ * scale_; + } + T out = 1.0; + if (x2 < 0.0) { + T x = std::sqrt(std::abs(x2)); + out = std::sinh(b_ * x) / x; + } else { + T x = std::sqrt(x2); + out = std::sin(b_ * x) / x; + } + return scale_ * out; +} + template TabulatedKernelView::TabulatedKernelView(const TabulatedKernel& kernel) : Base(kernel.width()), _table(kernel._table.data().get()), diff --git a/cxx/isce3/cuda/focus/Backproject.cu b/cxx/isce3/cuda/focus/Backproject.cu index 45117dcfd..6acdc02b2 100644 --- a/cxx/isce3/cuda/focus/Backproject.cu +++ b/cxx/isce3/cuda/focus/Backproject.cu @@ -3,7 +3,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -23,9 +25,17 @@ #include #include #include +#include #include #include +#include +#include +#include #include +#include +#include +#include +#include using namespace isce3::core; using namespace isce3::cuda::geometry; @@ -34,6 +44,10 @@ using isce3::cuda::core::interp1d; using isce3::error::ErrorCode; using isce3::focus::bistaticDelay; using isce3::focus::dryTropoDelayTSX; +using isce3::focus::PolarGrid; +using isce3::focus::setupPolarGridForPulses; +using isce3::cuda::signal::NFFT2dResult; +using isce3::cuda::signal::NFFT2dResultView; using HostDEMInterpolator = isce3::geometry::DEMInterpolator; using HostRadarGeometry = isce3::container::RadarGeometry; @@ -61,8 +75,6 @@ template using DeviceTabulatedKernel = isce3::cuda::core::TabulatedK namespace isce3 { namespace cuda { namespace focus { -namespace { - /** * \internal * Interpolate platform position and velocity at a range of uniformly-spaced @@ -370,6 +382,57 @@ __global__ void getCPIBounds(int* kstart_out, int* kstop_out, kstop_out[tid] = std::min(kstop, azimuth_time.size()); } +/** + * \internal + * Estimate coherent processing window bounds for one or more targets. + * + * Returns the indices of the first pulse and one past the last pulse to + * coherently integrate for each target. + * + * \param[out] tstart_out Processing window start time (inclusive) + * \param[out] tstop_out Processing window end time (exclusive) + * \param[in] t_in Azim. time of each target w.r.t. reference epoch (s) + * \param[in] r_in Slant range of each target (m) + * \param[in] x_in Position of each target in ECEF coords (m) + * \param[in] p_in Platform position at each target's azimuth time (m) + * \param[in] v_in Platform velocity at each target's azimuth time (m) + * \param[in] n Number of targets + * \param[in] wvl Radar wavelength (m) + * \param[in] ds Desired azimuth resolution (m) + */ +__global__ void getCPITimeBounds(double* tstart_out, double* tstop_out, + const double* t_in, const double* r_in, + const Vec3* x_in, const Vec3* p_in, + const Vec3* v_in, const size_t n, + const double wvl, const double ds) +{ + // thread index (1d grid of 1d blocks) + const auto tid = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + + // bounds check + if (tid >= n) { + return; + } + + // load inputs + const double t = t_in[tid]; + const double r = r_in[tid]; + const Vec3 p = p_in[tid]; + const Vec3 v = v_in[tid]; + const Vec3 x = x_in[tid]; + + // estimate synthetic aperture length required to achieve the desired + // azimuth resolution + const double l = wvl * r * (p.norm() / x.norm()) / (2. * ds); + + // approximate CPI duration (assuming constant platform velocity) + const double cpi = l / v.norm(); + + // get coherent processing window start & end time + tstart_out[tid] = t - 0.5 * cpi; + tstop_out[tid] = t + 0.5 * cpi; +} + /** * \internal * Backprojection core processing loop @@ -463,7 +526,41 @@ __global__ void sumCoherentBatch( out[tid] += thrust::complex(batch_sum); } -} // namespace +template +__global__ void broadcastMultiply(const T* row, size_t ncol, T* image, size_t npix) +{ + // thread index (1d grid of 1d blocks) + const auto tid = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + + // bounds check + if (tid >= npix) { + return; + } + + auto i_col = tid % ncol; + + image[tid] *= row[i_col]; +} + +__global__ void +makeSubApertureMask( + const double subaperture_start, const double subaperture_end, + const size_t n, + const double* pixel_start, + const double* pixel_end, + bool* mask) +{ + const auto tid = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + + // bounds check + if (tid >= n) { + return; + } + + mask[tid] = (subaperture_end > pixel_start[tid]) + and (subaperture_start < pixel_end[tid]); +} + template ErrorCode backproject(std::complex* out, @@ -753,4 +850,820 @@ ErrorCode backproject(std::complex* out, return ec; } + +// FBP junk + +/** + * \internal + * Transform a 2D radar grid from polar coordinates (cos_squint, range) to + * ECEF XYZ coordinates. + * + * The global error code is set if any thread encounters an error. + * + * \param[out] xyz_out ECEF XYZ of each target (m) + * \param[in] grid Polar grid + * \param[in] dem DEM sampling interface + * \param[in] ellipsoid Reference ellipsoid + * \param[in] side Radar look side + * \param[in] params Root-finding algorithm parameters + * \param[out] errc Error flag + */ +__global__ void runPolar2Geo(Vec3* xyz_out, const PolarGrid grid, + DeviceDEMInterpolator dem, const Ellipsoid ellipsoid, + const LookSide side, + const Rdr2GeoBracketParams params, + ErrorCode* errc) +{ + using isce3::geometry::detail::polar2geo_bracket; + + // thread index (1d grid of 1d blocks) + const auto tid = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + + // bounds check + const auto lines = static_cast(grid.length()); + const auto samples = static_cast(grid.width()); + if (tid >= lines * samples) { + return; + } + + // convert flat index to 2D array indices + const auto j = static_cast(tid / samples); + const auto i = static_cast(tid % samples); + + const double r = grid.range[i]; + const double q = grid.sin_squint[j]; + const double c = sqrt(1.0 - q * q); + + Vec3 xyz; + double look_angle; + + const auto status = polar2geo_bracket(&xyz, &look_angle, + grid.origin, grid.axis, r, q, c, dem, ellipsoid, + side, params); + + // check convergence + if (status == isce3::error::ErrorCode::Success) { + xyz_out[tid] = xyz; + } else { + // set output to NaN + constexpr static auto nan = std::numeric_limits::quiet_NaN(); + xyz_out[tid] = {nan, nan, nan}; + + // set global error flag + *errc = ErrorCode::FailedToConverge; + } +} + + +template +class TimingReporter { +public: + TimingReporter(const std::string& channel_id) : + prev_time_{Clock::now()}, log_{channel_id} {} + + void report(const std::string& step) { + using namespace std::chrono; + auto cur_time = Clock::now(); + auto duration = cur_time - prev_time_; + auto msec = duration_cast(duration).count(); + log_ << step << " took " << msec << " ms" << pyre::journal::endl; + prev_time_ = cur_time; + } +private: + std::chrono::time_point prev_time_; + pyre::journal::info_t log_; +}; + +// macro to enable logging of timing info +#ifdef ISCE3_ENABLE_FBP_TIMING +#define ISCE3_FBP_TIMING(x) x +#else +#define ISCE3_FBP_TIMING(x) // no-op +#endif + +template +std::tuple[]>, std::unique_ptr> +backprojectToPolarGrid( + const std::complex* in, + const Linspace& in_slant_range, + const std::vector& pos, + const std::vector& vel, + const PolarGrid& out_grid, + DeviceDEMInterpolator& dem, double fc, + const Kernel& kernel, DryTroposphereModel dry_tropo_model, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params) +{ + const auto nt = static_cast(pos.size()); + if (nt < 2) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "require at least two pulses in FBP stage"); + } + if (vel.size() != nt) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "require same number of position and velocity vectors"); + } + + static constexpr double c = isce3::core::speed_of_light; + + // check that dry_tropo_model is supported internally + if (not(dry_tropo_model == DryTroposphereModel::NoDelay or + dry_tropo_model == DryTroposphereModel::TSX)) { + + std::string errmsg = "unexpected dry troposphere model"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + ISCE3_FBP_TIMING( + auto timing = TimingReporter("isce3.cuda.focus.backprojectToPolarGrid"); + ) + + const auto npix = static_cast(out_grid.length()) * out_grid.width(); + auto height = std::make_unique(npix); + auto out = std::make_unique[]>(npix); + + // range sampling window + double swst = 2. * in_slant_range.first() / c; + double dtau = 2. * in_slant_range.spacing() / c; + int nr = in_slant_range.size(); + Linspace sampling_window(swst, dtau, nr); + + // reference ellipsoid + int epsg = dem.epsgCode(); + const Ellipsoid ellipsoid = makeProjection(epsg)->ellipsoid(); + + // init device variable to return error codes from device code + thrust::device_vector errc(1, ErrorCode::Success); + + thrust::device_vector x(npix); + + { + const unsigned block = 256; + const unsigned grid = (npix + block - 1) / block; + + runPolar2Geo<<>>(x.data().get(), out_grid, + dem, ellipsoid, + out_grid.look_side, r2g_params, + errc.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + ISCE3_FBP_TIMING(timing.report("polar2geo");) + + // transform each target position from ECEF to LLH coordinates + // NOTE only really needed if dumping height layer or doing TSX atmosphere + // correction, but just compute it unconditionally. + thrust::device_vector llh(npix); + + { + const unsigned block = 256; + const unsigned grid = (npix + block - 1) / block; + + ecef2llh<<>>(llh.data().get(), x.data().get(), + npix, ellipsoid); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + if (height != nullptr) { + thrust::device_vector d_height(npix); + thrust::transform(llh.begin(), llh.end(), d_height.begin(), + [] __device__ (const Vec3& x) { return (float)x[2]; }); + checkCudaErrors(cudaMemcpy(height.get(), d_height.data().get(), + npix * sizeof(float), cudaMemcpyDeviceToHost)); + } + ISCE3_FBP_TIMING(timing.report("ecef2llh");) + + // estimate dry troposphere delay + thrust::device_vector tau_atm(npix); + + if (dry_tropo_model == DryTroposphereModel::NoDelay) { + checkCudaErrors(cudaMemset(tau_atm.data().get(), 0, + npix * sizeof(double))); + } else if (dry_tropo_model == DryTroposphereModel::TSX) { + const unsigned block = 256; + const unsigned grid = (npix + block - 1) / block; + + // TODO new interface for constant aperture center + thrust::device_vector p(npix); + thrust::fill(p.begin(), p.end(), out_grid.origin); + + estimateDryTropoDelayTSX<<>>( + tau_atm.data().get(), p.data().get(), llh.data().get(), + npix, ellipsoid); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } else { + std::string errmsg = "unexpected dry troposphere model"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), errmsg); + } + ISCE3_FBP_TIMING(timing.report("dry troposphere");) + + // Assume we can fit all pulses for a subimage in device memory. + const auto npix_in = static_cast(in_slant_range.size()) * nt; + thrust::device_vector> rc(npix_in); + checkCudaErrors(cudaMemcpy(rc.data().get(), in, npix_in * sizeof(*in), + cudaMemcpyHostToDevice)); + + thrust::device_vector> img(npix, 0.0); + { + // integrate pulses + const unsigned block = 256; + const unsigned grid = (npix + block - 1) / block; + + using KV = typename Kernel::view_type; + + thrust::device_vector d_pos(pos); + thrust::device_vector d_vel(vel); + + // TODO interface with scalar kstart & kstop + thrust::device_vector kstart(npix, 0); + thrust::device_vector kstop(npix, nt); + + sumCoherentBatch<<>>( + img.data().get(), rc.data().get(), d_pos.data().get(), + d_vel.data().get(), sampling_window, x.data().get(), + tau_atm.data().get(), kstart.data().get(), kstop.data().get(), + npix, fc, kernel, 0, nt); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + ISCE3_FBP_TIMING(timing.report("sum coherent");) + + // copy output back to the host + checkCudaErrors(cudaMemcpy(out.get(), img.data().get(), + npix * sizeof(std::complex), + cudaMemcpyDeviceToHost)); + + + // baseband + const double kw = 4 * M_PI / (c / fc); + #pragma omp parallel for + for (int i = 0; i < out_grid.range.size(); ++i) { + const double phi = -kw * out_grid.range[i]; + const auto phasor = std::complex(std::cos(phi), std::sin(phi)); + for (int j = 0; j < out_grid.sin_squint.size(); ++j) { + out[j * out_grid.width() + i] *= phasor; + } + } + ISCE3_FBP_TIMING(timing.report("baseband");) + + return std::make_tuple(errc[0], std::move(out), std::move(height)); +} + + +std::tuple< + ErrorCode, + std::unique_ptr[]>, + std::unique_ptr> +backprojectToPolarGrid( + const std::complex* in, const Linspace& in_slant_range, + const std::vector& pos, const std::vector& vel, + const PolarGrid& out_grid, + const HostDEMInterpolator& dem, double fc, + const Kernel& kernel, DryTroposphereModel dry_tropo_model, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params) +{ + DeviceDEMInterpolator d_dem(dem); + + if (typeid(kernel) == typeid(HostBartlettKernel)) { + const DeviceBartlettKernel d_kernel( + dynamic_cast&>(kernel)); + return backprojectToPolarGrid(in, in_slant_range, pos, vel, + out_grid, d_dem, fc, d_kernel, dry_tropo_model, r2g_params); + } + else if (typeid(kernel) == typeid(HostLinearKernel)) { + const DeviceLinearKernel d_kernel( + dynamic_cast&>(kernel)); + return backprojectToPolarGrid(in, in_slant_range, pos, vel, + out_grid, d_dem, fc, d_kernel, dry_tropo_model, r2g_params); + } + else if (typeid(kernel) == typeid(HostKnabKernel)) { + const DeviceKnabKernel d_kernel( + dynamic_cast&>(kernel)); + return backprojectToPolarGrid(in, in_slant_range, pos, vel, + out_grid, d_dem, fc, d_kernel, dry_tropo_model, r2g_params); + } + else if (typeid(kernel) == typeid(HostTabulatedKernel)) { + const DeviceTabulatedKernel d_kernel( + dynamic_cast&>(kernel)); + return backprojectToPolarGrid(in, in_slant_range, pos, vel, + out_grid, d_dem, fc, d_kernel, dry_tropo_model, r2g_params); + } + else if (typeid(kernel) == typeid(HostChebyKernel)) { + const DeviceChebyKernel d_kernel( + dynamic_cast&>(kernel)); + return backprojectToPolarGrid(in, in_slant_range, pos, vel, + out_grid, d_dem, fc, d_kernel, dry_tropo_model, r2g_params); + } + throw isce3::except::RuntimeError(ISCE_SRCINFO(), "not implemented"); +} + +__global__ void +interpPolar(thrust::complex* geo_image, const Vec3* geo_points, + size_t n, const PolarGrid grid, + const NFFT2dResultView nfft, const double kw, + std::optional mask = std::nullopt, + std::optional extra_range_delays = std::nullopt) +{ + // thread index (1d grid of 1d blocks) + const auto i = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + + // bounds check + if (i >= n) { + return; + } + // mask check + if (mask.has_value() and not mask.value()[i]) { + return; + } + + // compute target location in polar grid + // TODO not sure why I get a linker error when I try using the API (with CUDA_HOSTDEV added) + // double sin_squint, range; + // isce3::geometry::geo2polar(&sin_squint, &range, geo_points[i], grid.origin, grid.axis); + const Vec3 lookvec = geo_points[i] - grid.origin; + double range = lookvec.norm(); + const double sin_squint = lookvec.dot(grid.axis) / range; + + if (extra_range_delays.has_value()) { + range += extra_range_delays.value()[i]; + } + + // convert to image index + const double ix = (range - grid.range.first()) / grid.range.spacing(), + iy = (sin_squint - grid.sin_squint.first()) / grid.sin_squint.spacing(); + + // interpolate baseband data + const auto z = nfft.interp({iy, ix}, /* periodic */ false); + + // compensate phase and sum contribution + double sin_phi, cos_phi; + ::sincos(kw * range, &sin_phi, &cos_phi); + geo_image[i] += z * thrust::complex(cos_phi, sin_phi); +} + +ErrorCode +projectPolarToGeo( + std::complex* geo_image, + const Vec3* geo_points, + const size_t n, + const PolarGrid& grid, + const std::complex* polar_image, + const double wavelength, + const isce3::signal::NFFT2dParams& params) +{ + using isce3::fft::nextFastPower; + using dims_t = isce3::cuda::signal::NFFT2d::dims_t; + using std::lround; + + const dims_t m = {params.rows.m, params.cols.m}; + + // NFFTKernel width = 2*m+1; interp2d uses fixed-size stack arrays of + // MAX_WIDTH=16, so reject m >= 8 here at the API entry point instead of + // corrupting memory in the device kernel. + constexpr int MAX_NFFT_M = 7; + if (m[0] > MAX_NFFT_M or m[1] > MAX_NFFT_M) { + return ErrorCode::InvalidKernelSize; + } + + const dims_t dims_in = {grid.length(), grid.width()}; + const dims_t dims_out = { + nextFastPower(static_cast(lround(params.rows.s * dims_in[0]))), + nextFastPower(static_cast(lround(params.cols.s * dims_in[1]))) + }; + + auto nfft = isce3::cuda::signal::NFFT2d(m, dims_in, dims_out); + const size_t nin = static_cast(grid.length()) * grid.width(); + std::vector> spectrum(nin); + // Okay to discard const because fft is planned with FFTW_EXECUTE which + // doesn't modify input. + ISCE3_FBP_TIMING( + auto timing = TimingReporter("isce3.cuda.focus.backprojectToPolarGrid");) + // TODO do this FFT on GPU + isce3::fft::fft2d(spectrum.data(), + const_cast*>(polar_image), + {dims_in[0], dims_in[1]}); + ISCE3_FBP_TIMING(timing.report("FFT");) + + // zero-pad and filter + auto result = nfft.transform_host(dims_in, /* strides = */ {dims_in[1], 1}, + spectrum.data()); + auto nfft_view = NFFT2dResultView(result); + ISCE3_FBP_TIMING(timing.report("IFFT");) + + const double kw = 4 * M_PI / wavelength; + + // TODO add interface that just accepts device pointers as argument? + // NOTE You can construct device_vector using iterators, but that's super + // slow, hence the manual copy. + thrust::device_vector> d_geo_image(n); + thrust::device_vector d_geo_points(n); + checkCudaErrors(cudaMemcpy(d_geo_image.data().get(), geo_image, + n * sizeof(*geo_image), cudaMemcpyHostToDevice)); + checkCudaErrors(cudaMemcpy(d_geo_points.data().get(), geo_points, + n * sizeof(*geo_points), cudaMemcpyHostToDevice)); + ISCE3_FBP_TIMING(timing.report("copy to device");) + + { + const unsigned block = 256; + const unsigned cugrid = (n + block - 1) / block; + + interpPolar<<>>(d_geo_image.data().get(), + d_geo_points.data().get(), n, grid, nfft_view, kw); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + ISCE3_FBP_TIMING(timing.report("interp");) + + checkCudaErrors(cudaMemcpy(geo_image, d_geo_image.data().get(), + n * sizeof(*geo_image), cudaMemcpyDeviceToHost)); + ISCE3_FBP_TIMING(timing.report("copy to host");) + + return ErrorCode::Success; +} + + +isce3::error::ErrorCode +accumulatePolarImagesToRadarGrid(std::complex* out, + const isce3::container::RadarGeometry& out_geometry, + const isce3::core::Orbit& in_orbit, + const isce3::core::LUT2d& in_doppler, + const std::vector& grids, + const std::vector*>& image_interpolators, + const isce3::geometry::DEMInterpolator& dem, double fc, double ds, + DryTroposphereModel dry_tropo_model, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params, + float* height) +{ + using namespace isce3::core; + using namespace isce3::cuda::geometry; + using isce3::focus::dryTropoDelayTSX; + using isce3::error::ErrorCode; + using isce3::focus::PolarGrid; + using DeviceDEMInterpolator = isce3::cuda::geometry::gpuDEMInterpolator; + using DeviceRadarGeometry = isce3::cuda::container::RadarGeometry; + + static constexpr double c = isce3::core::speed_of_light; + static constexpr auto nan = std::numeric_limits::quiet_NaN(); + + // will search sorted intervals to figure out active sub images per target + auto starts = std::vector(grids.size()); + std::transform(grids.begin(), grids.end(), starts.begin(), + [](const PolarGrid& grid) { return grid.aztime_start; }); + auto ends = std::vector(grids.size()); + std::transform(grids.begin(), grids.end(), ends.begin(), + [](const PolarGrid& grid) { return grid.aztime_end; }); + + // get input & output radar grid azimuth time & slant range + Linspace out_azimuth_time = out_geometry.sensingTime(); + Linspace out_slant_range = out_geometry.slantRange(); + + // reference ellipsoid + int epsg = dem.epsgCode(); + Ellipsoid ellipsoid = makeProjection(epsg)->ellipsoid(); + + // carrier wavelength + const double wvl = c / fc; + const double kw = 4 * M_PI / wvl; + + // copy inputs to device + const DeviceRadarGeometry d_out_geometry(out_geometry); + DeviceDEMInterpolator d_dem(dem); + + const size_t nout = out_geometry.gridLength() * out_geometry.gridWidth(); + + thrust::device_vector d_x(nout); + thrust::device_vector errc(1, ErrorCode::Success); + { + const unsigned block = 256; + const unsigned grid = (nout + block - 1) / block; + + runRdr2Geo<<>>(d_x.data().get(), out_azimuth_time, + out_slant_range, d_out_geometry.doppler(), + d_out_geometry.orbit(), d_dem, ellipsoid, wvl, + d_out_geometry.lookSide(), r2g_params, errc.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + // transform each target position from ECEF to LLH coordinates + // NOTE only really needed if dumping height layer or doing TSX atmosphere + // correction, but just compute it unconditionally. + thrust::device_vector d_llh(nout); + + { + const unsigned block = 256; + const unsigned grid = (nout + block - 1) / block; + + ecef2llh<<>>(d_llh.data().get(), d_x.data().get(), + nout, ellipsoid); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + if (height != nullptr) { + thrust::device_vector d_height(nout); + thrust::transform(d_llh.begin(), d_llh.end(), d_height.begin(), + [] __device__ (const Vec3& x) { return (float)x[2]; }); + checkCudaErrors(cudaMemcpy(height, d_height.data().get(), + nout * sizeof(float), cudaMemcpyDeviceToHost)); + } + + // Running geo2rdr to get integration bounds seems like overkill. + // TODO Maybe mask on Doppler instead? + thrust::device_vector d_t(nout); + thrust::device_vector d_r(nout); + auto d_in_orbit = isce3::cuda::core::Orbit(in_orbit); + auto d_in_orbit_view = isce3::cuda::core::OrbitView(d_in_orbit); + auto d_in_doppler = isce3::cuda::core::gpuLUT2d(in_doppler); + + { + const unsigned block = 256; + const unsigned grid = (nout + block - 1) / block; + + runGeo2Rdr<<>>( + d_t.data().get(), d_r.data().get(), d_x.data().get(), nout, + d_in_orbit_view, d_in_doppler, wvl, + d_out_geometry.lookSide(), g2r_params, errc.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + // get platform position & velocity at center of CPI for each target + thrust::device_vector d_p(nout); + thrust::device_vector d_v(nout); + + { + const unsigned block = 256; + const unsigned grid = (nout + block - 1) / block; + + interpolateOrbit<<>>(d_p.data().get(), d_v.data().get(), + d_in_orbit_view, d_t.data().get(), + nout, errc.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + // Calculate dry troposphere delay. To re-use code and memory, first we'll + // calculate in time units and then convert to spatial units. + thrust::device_vector dr_atm(nout); + + if (dry_tropo_model == DryTroposphereModel::NoDelay) { + checkCudaErrors(cudaMemset(dr_atm.data().get(), 0, + nout * sizeof(double))); + } else if (dry_tropo_model == DryTroposphereModel::TSX) { + const unsigned block = 256; + const unsigned grid = (nout + block - 1) / block; + + estimateDryTropoDelayTSX<<>>( + dr_atm.data().get() /* time units */, d_p.data().get(), + d_llh.data().get(), nout, ellipsoid); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } else { + std::string errmsg = "unexpected dry troposphere model"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + // two-way time -> one-way range + constexpr double halfspeed = isce3::core::speed_of_light / 2.0; + thrust::transform(dr_atm.begin(), dr_atm.end(), dr_atm.begin(), + [halfspeed] __host__ __device__ (double dt) { return halfspeed * dt; }); + + d_llh.clear(); d_llh.shrink_to_fit(); + + // get coherent integration bounds (pulse indices) for each target + thrust::device_vector d_tstart(nout); + thrust::device_vector d_tstop(nout); + + { + const unsigned block = 256; + const unsigned grid = (nout + block - 1) / block; + + getCPITimeBounds<<>>(d_tstart.data().get(), + d_tstop.data().get(), d_t.data().get(), d_r.data().get(), + d_x.data().get(), d_p.data().get(), d_v.data().get(), nout, + wvl, ds); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + d_p.clear(); d_p.shrink_to_fit(); + d_v.clear(); d_v.shrink_to_fit(); + d_t.clear(); d_t.shrink_to_fit(); + d_r.clear(); d_r.shrink_to_fit(); + + // NOTE Thrust does not use the bit-packing strategy as STL. + thrust::device_vector d_mask(nout); + + // TODO reduce tstart & tend + // TODO check this O(log(n)) algorithm + //const auto kstart = std::distance(ends.begin(), + // std::lower_bound(ends.begin(), ends.end(), tstart)); + //const auto kstop = std::distance(starts.begin(), + // std::upper_bound(starts.start(), starts.end(), tstart + cpi)); + const auto num_images = image_interpolators.size(); + const decltype(num_images) kstart = 0, kstop = num_images; + + // Copy image to device since we accumulate (don't init to zero). + // thrust::device_vector> d_out(out, out + nout); + thrust::device_vector> d_out(nout); + checkCudaErrors(cudaMemcpy(d_out.data().get(), out, nout * sizeof(*out), + cudaMemcpyHostToDevice)); + + // Advance forward iterator as needed. + auto image_iter = image_interpolators.begin(); + for (int k = 0; k < kstart; ++k) ++image_iter; + + for (auto k = kstart; k < kstop; ++k, ++image_iter) { + const auto& image_grid = grids[k]; + const auto d_nfft = *image_iter; + + { + const unsigned block = 256; + const unsigned cuda_grid = (nout + block - 1) / block; + + makeSubApertureMask<<>>(image_grid.aztime_start, + image_grid.aztime_end, nout, d_tstart.data().get(), + d_tstop.data().get(), d_mask.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + const auto d_nfft_view = NFFT2dResultView(*d_nfft); + + { + const unsigned block = 256; + const unsigned cuda_grid = (nout + block - 1) / block; + + interpPolar<<>>(d_out.data().get(), + d_x.data().get(), nout, image_grid, d_nfft_view, kw, + d_mask.data().get(), dr_atm.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + } + + // Copy result back to host. + checkCudaErrors(cudaMemcpy(out, d_out.data().get(), nout * sizeof(*out), + cudaMemcpyDeviceToHost)); + + return errc[0]; +} + + +void mergePolarImages( + const std::vector& grids, + const std::vector*>& image_interpolators, + const isce3::focus::PolarGrid& output_grid, + Eigen::Ref>> output_image, + const double fc, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + int az_block_size) +{ + // check that output grid dimensions match buffer size + const auto m = output_grid.length(), n = output_grid.width(); + if ((m != output_image.rows()) or (n != output_image.cols())) { + std::string msg = "Dimensions of image grid (" + std::to_string(m) + + ", " + std::to_string(n) + ") do not match dimensions of image " + "buffer (" + std::to_string(output_image.rows()) + ", " + + std::to_string(output_image.cols()) + ")"; + throw isce3::except::LengthError(ISCE_SRCINFO(), msg); + } + + // check that we have a grid for each input image + const auto num_images = image_interpolators.size(); + if (grids.size() != num_images) { + std::string msg = "Size mismatch: got " + std::to_string(num_images) + + " sub images but " + std::to_string(grids.size()) + " grids"; + throw isce3::except::LengthError(ISCE_SRCINFO(), msg); + } + + // check look directions for consistency + const auto look_side = output_grid.look_side; + for (const auto& grid : grids) { + if (grid.look_side != look_side) { + std::string msg = "Output grid look direction does not match " + "input grid look direction"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), msg); + } + } + + // Check block size and allocate scratch space. + if (az_block_size <= 0) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "azimuth block size must be positive"); + } + az_block_size = std::min(az_block_size, output_grid.sin_squint.size()); + + // auto block_positions = isce3::core::EArray2D(); + const auto npix = static_cast(az_block_size) * output_grid.width(); + auto block_positions = thrust::device_vector(npix); + auto block_image = thrust::device_vector>(npix); + + // reference ellipsoid + isce3::core::Ellipsoid ellipsoid = isce3::core::makeProjection(dem.epsgCode())->ellipsoid(); + isce3::cuda::geometry::gpuDEMInterpolator d_dem(dem); + + // wavenumber + const double kw = 4 * M_PI * fc / isce3::core::speed_of_light; + + // Baseband. Note that we could do this at the same time as the + // reprojection but it'd require a fair bit of copy/paste. + thrust::host_vector> h_phasors(n); + #pragma omp parallel for + for (auto j = decltype(n){0}; j < n; ++j) { + const double arg = -kw * output_grid.range[j]; + h_phasors[j] = std::complex(std::cos(arg), std::sin(arg)); + } + const auto d_phasors = thrust::device_vector>(h_phasors); + + using isce3::error::ErrorCode; + thrust::device_vector errc(1, ErrorCode::Success); + + // loop over output blocks + auto n_blocks = (m + az_block_size - 1) / az_block_size; + for (auto i_block = decltype(n_blocks){0}; i_block < n_blocks; ++i_block) { + const auto i_row0 = i_block * az_block_size; + const auto i_row1 = std::min(i_row0 + az_block_size, m); + const auto block_npix = static_cast(n) * (i_row1 - i_row0); + + // Zero out block_image since interpPolar accumulates + checkCudaErrors(cudaMemset(block_image.data().get(), 0, + block_npix * sizeof(thrust::complex))); + + const auto output_grid_subset = output_grid.offsetAndResize(i_row0, 0, + i_row1 - i_row0, output_grid.width()); + + // Compute output pixel 3D locations + { + const unsigned cu_block = 256; + const unsigned cu_grid = (block_npix + cu_block - 1) / cu_block; + + runPolar2Geo<<>>(block_positions.data().get(), + output_grid_subset, d_dem, ellipsoid, output_grid.look_side, + r2g_params, errc.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + const auto ec = errc[0]; + if (ec != ErrorCode::Success) { + throw isce3::except::DomainError(ISCE_SRCINFO(), + "polar2geo failed with ErrorCode " + + isce3::error::getErrorString(ec)); + } + + // loop over input images + auto image_it = image_interpolators.begin(); + for (int i_img = 0; i_img < num_images; ++i_img, ++image_it) { + const auto& input_grid = grids[i_img]; + const auto d_nfft = *image_it; + const auto d_nfft_view = NFFT2dResultView(*d_nfft); + { + const unsigned cu_block = 256; + const unsigned cu_grid = (block_npix + cu_block - 1) / cu_block; + + interpPolar<<>>(block_image.data().get(), + block_positions.data().get(), block_npix, input_grid, + d_nfft_view, kw); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + } // images + + // baseband + { + const unsigned cu_block = 256; + const unsigned cu_grid = (block_npix + cu_block - 1) / cu_block; + + broadcastMultiply<<>>(d_phasors.data().get(), + n, block_image.data().get(), block_npix); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + checkCudaErrors(cudaMemcpy(output_image.row(i_row0).data(), + block_image.data().get(), block_npix * sizeof(std::complex), + cudaMemcpyDeviceToHost)); + } // blocks +} + }}} // namespace isce3::cuda::focus diff --git a/cxx/isce3/cuda/focus/Backproject.h b/cxx/isce3/cuda/focus/Backproject.h index 260d23846..29b37694e 100644 --- a/cxx/isce3/cuda/focus/Backproject.h +++ b/cxx/isce3/cuda/focus/Backproject.h @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -87,4 +88,154 @@ backproject(std::complex* out, const Geo2RdrBracketParams& geo2rdr_params = {}, int batch = 1024, float* height = nullptr); + +/** + * @brief Backproject range-compressed signal into a polar grid (GPU). + * + * Performs time-domain backprojection of range-compressed SAR signal + * data onto a polar coordinate grid. For each output pixel in the polar + * grid, the signal is resampled from the input data by accumulating + * contributions from all pulses according to the instantaneous slant + * range, with phase compensation for motion and the troposphere. + * + * @param[in] in Input range-compressed signal data + * @param[in] in_slant_range Slant range grid of the input data (m) + * @param[in] pos Platform position vectors at each pulse + * (ECEF, m) + * @param[in] vel Platform velocity vectors at each pulse + * (ECEF, m/s) + * @param[in] out_grid Target polar grid to backproject onto + * @param[in] dem Digital elevation model (DEM) + * @param[in] fc Center frequency (Hz) + * @param[in] kernel 1-D interpolation kernel + * @param[in] dry_tropo_model Dry troposphere path delay model + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * + * @returns A tuple containing: + * - error code (non-zero if geometry fails to converge for + * any pixel, in which case values for those pixels are NaN) + * - focused signal data on the polar grid (size = + * out_grid.width() * out_grid.length()) + * - per-pixel height above the ellipsoid (m) + */ +std::tuple< + isce3::error::ErrorCode, + std::unique_ptr[]>, // image + std::unique_ptr> // height +backprojectToPolarGrid( + const std::complex* in, + const isce3::core::Linspace& in_slant_range, + const std::vector& pos, + const std::vector& vel, + const isce3::focus::PolarGrid& out_grid, + const isce3::geometry::DEMInterpolator& dem, + double fc, + const isce3::core::Kernel& kernel, + DryTroposphereModel dry_tropo_model, + const Rdr2GeoBracketParams& r2g_params = {}); + +/** + * @brief Project a polar grid image to geographic points (GPU). + * + * Interpolates a polar grid image to a set of 3D geographic (XYZ) + * positions using NFFT-based interpolation. + * + * @param[out] geo_image Output interpolated signal at geo points + * @param[in] geo_points Target 3D positions (ECEF, m) + * @param[in] n Number of target positions + * @param[in] grid Polar grid containing the image data + * @param[in] polar_image Input polar grid image data + * @param[in] wavelength Radar wavelength (m) + * @param[in] params NFFT interpolation parameters + * + * @returns Error code indicating success or failure + */ +isce3::error::ErrorCode +projectPolarToGeo( + std::complex* geo_image, + const isce3::core::Vec3* geo_points, + const size_t n, + const isce3::focus::PolarGrid& grid, + const std::complex* polar_image, + const double wavelength, + const isce3::signal::NFFT2dParams& params); + +/** + * @brief Accumulate polar grid images onto an output radar grid (GPU). + * + * Combines multiple subaperture polar grid images back onto a + * Cartesian radar geometry grid. For each pixel in the output grid, + * the target position is computed via rdr2geo, the corresponding + * coherent processing interval is determined via geo2rdr, and the + * polar image data is accumulated using NFFT-based interpolation. + * + * The caller is responsible for allocating the output arrays to the + * appropriate size (out_geometry.gridLength() * + * out_geometry.gridWidth()). + * + * @param[out] out Accumulated focused signal data + * @param[in] out_geometry Target output grid, orbit, and Doppler + * @param[in] in_orbit Input data orbit + * @param[in] in_doppler Input data Doppler centroid LUT + * @param[in] grids List of subaperture polar grids + * @param[in] image_interpolators NFFT interpolators for each grid + * @param[in] dem Digital elevation model (DEM) + * @param[in] fc Center frequency (Hz) + * @param[in] ds Desired azimuth resolution (m) + * @param[in] dry_tropo_model Dry troposphere path delay model + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * @param[in] g2r_params geo2rdr_bracket configuration parameters + * @param[out] height Height of each pixel (m) above the + * ellipsoid (optional, may be nullptr) + * + * @returns Non-zero error code if rdr2geo or geo2rdr fails to + * converge for any pixel, and the values for these + * pixels are set to NaN. + */ +isce3::error::ErrorCode +accumulatePolarImagesToRadarGrid(std::complex* out, + const isce3::container::RadarGeometry& out_geometry, + const isce3::core::Orbit& in_orbit, + const isce3::core::LUT2d& in_doppler, + const std::vector& grids, + const std::vector*>& + image_interpolators, + const isce3::geometry::DEMInterpolator& dem, double fc, double ds, + DryTroposphereModel dry_tropo_model = DryTroposphereModel::TSX, + const Rdr2GeoBracketParams& r2g_params = {}, + const Geo2RdrBracketParams& g2r_params = {}, + float* height = nullptr); + +/** + * @brief Merge subaperture polar grid images into a single grid (GPU). + * + * Combines multiple subaperture polar grid images onto a merged + * output polar grid. For each pixel in the output grid, the 3D + * target position is computed via polar2geo, and the input image + * data is accumulated via NFFT-based interpolation. The output + * image is expected to be zero-initialized by the caller. + * + * @param[in] grids List of subaperture input polar grids + * @param[in] image_interpolators NFFT interpolators for each input grid + * @param[in] output_grid Merged output polar grid + * @param[out] output_image Accumulated output image (must be + * zero-initialized); dimensions must + * match output_grid + * @param[in] fc Center frequency (Hz) + * @param[in] dem Digital elevation model (DEM) + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * @param[in] az_block_size Number of azimuth rows to process + * at a time (defaults to 1024) + */ +void mergePolarImages( + const std::vector& grids, + const std::vector*>& + image_interpolators, + const isce3::focus::PolarGrid& output_grid, + Eigen::Ref>> output_image, + const double fc, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params = {}, + int az_block_size = 1024); + }}} // namespace isce3::cuda::focus diff --git a/cxx/isce3/cuda/signal/NFFT2d.cu b/cxx/isce3/cuda/signal/NFFT2d.cu new file mode 100644 index 000000000..4f908c908 --- /dev/null +++ b/cxx/isce3/cuda/signal/NFFT2d.cu @@ -0,0 +1,267 @@ +#include "NFFT2d.h" + +#include +#include +#include + +template +using Kernel = isce3::cuda::core::NFFTKernel; + +namespace isce3::cuda::signal { + +// constructor +template +NFFT2d::NFFT2d(const dims_t& m, const dims_t& sizes, const dims_t& fft_sizes) + : m_(m), sizes_(sizes), fft_sizes_(fft_sizes), kernels_( + {Kernel{m[0], sizes[0], fft_sizes[0]}, + Kernel{m[1], sizes[1], fft_sizes[1]}}) +{ + size_t nout = static_cast(fft_sizes[0]) * fft_sizes[1]; + xf_.resize(nout); + + // Just compute weights on CPU for now and then copy to GPU. + + // Pre-compute spectral weights (1/phi_hat in NFFT papers). + // Also include factor of n since FFTW does not normalize DFT. + for (int idim = 0; idim < ndims; ++idim) { + auto weight = std::vector(sizes[idim]); + T b = M_PI * (2.0 - 1.0 * sizes[idim] / fft_sizes[idim]); + T norm = isce3::math::bessel_i0(b * m[idim]) / sizes[idim]; + size_t n2 = (sizes[idim] - 1) / 2 + 1; + for (size_t i = 0; i < n2; ++i) { + double f = 2 * M_PI * i / fft_sizes_[idim]; + weight[i] = norm / + isce3::math::bessel_i0(m[idim] * std::sqrt(b * b - f * f)); + } + for (size_t i = n2; i < sizes[idim]; ++i) { + double f = 2 * M_PI * ((double)i - sizes[idim]) / fft_sizes[idim]; + weight[i] = norm / + isce3::math::bessel_i0(m[idim] * std::sqrt(b * b - f * f)); + } + weights_[idim].assign(weight.begin(), weight.end()); + } +} + +template +__global__ void setSpectrum2d(thrust::complex* xout, int rows_out, int cols_out, + const thrust::complex* xin, int rows_in, int cols_in, int row_stride_in, + int col_stride_in, T* weights_rows, T* weights_cols) +{ + int col = blockIdx.x * blockDim.x + threadIdx.x; + int row = blockIdx.y * blockDim.y + threadIdx.y; + + if ((col >= cols_in) || (row >= rows_in)) { + return; + } + + // NOTE For even lengths we're not splitting Nyquist bins. + const auto kin = static_cast(row) * row_stride_in + + col * col_stride_in; + + long row_out = 0, col_out = 0; + + int m2 = rows_in / 2; + int n2 = cols_in / 2; + + if (rows_in % 2 == 1) { + ++m2; + } + if (cols_in % 2 == 1) { + ++n2; + } + + if (row < m2) { + row_out = row; + } else { + row_out = rows_out - (rows_in - row); + } + if (col < n2) { + col_out = col; + } else { + col_out = cols_out - (cols_in - col); + } + + const auto kout = row_out * cols_out + col_out; + xout[kout] = weights_rows[row] * weights_cols[col] * xin[kin]; +} + +// Digest some data. +template +NFFT2dResult +NFFT2d::transform_host(const dims_t& sizes, const dims_t& strides, const std::complex *x) +{ + if ((strides[0] != sizes[1]) or (strides[1] != 1)) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Only implemented for C-ordered data"); + } + // Copy input data to device. + auto nin = static_cast(sizes[0]) * sizes[1]; + thrust::device_vector> d_x(nin); + checkCudaErrors(cudaMemcpy(d_x.data().get(), x, nin * sizeof(*x), + cudaMemcpyHostToDevice)); + + return NFFT2d::transform_device(sizes, strides, d_x.data().get()); +} + +template +NFFT2dResult +NFFT2d::transform_device(const dims_t& sizes, const dims_t& strides, const thrust::complex* x) +{ + for (int idim = 0; idim < ndims; ++idim) { + if (sizes[idim] != sizes_[idim]) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "Spectrum size != NFFT size."); + } + } + // Clear any old data. + auto nout = static_cast(fft_sizes_[0]) * fft_sizes_[1]; + xf_.assign(nout, thrust::complex(0, 0)); + + // Pad and weight + { + dim3 cu_block(16, 16); + dim3 cu_grid( + (sizes[1] + cu_block.x - 1) / cu_block.x, + (sizes[0] + cu_block.y - 1) / cu_block.y); + + setSpectrum2d<<>>( + xf_.data().get(), fft_sizes_[0], fft_sizes_[1], + x, sizes[0], sizes[1], strides[0], strides[1], + weights_[0].data().get(), weights_[1].data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaStreamSynchronize(cudaStreamDefault)); + } + + // Transform to time domain. + auto result = NFFT2dResult(m_, sizes_, fft_sizes_, kernels_); + int dims[] = {fft_sizes_[0], fft_sizes_[1]}; + isce3::cuda::fft::ifft2d(result.xt_.data().get(), xf_.data().get(), dims); + + return result; +} + +template +NFFT2dResult::operator isce3::signal::NFFT2dResult() const +{ + using CpuKernel = isce3::core::NFFTKernel; + auto result = isce3::signal::NFFT2dResult(m_, sizes_, fft_sizes_, + {CpuKernel {kernels_[0]}, CpuKernel {kernels_[1]}}); + const auto npix = static_cast(fft_sizes_[0]) * fft_sizes_[1]; + checkCudaErrors(cudaMemcpy(result.data(), xt_.data().get(), + npix * sizeof(std::complex), cudaMemcpyDeviceToHost)); + return result; +} + +template +NFFT2dResult::NFFT2dResult(const isce3::signal::NFFT2dResult& other) + : m_ {other.kernel_radii()}, sizes_ {other.sizes()}, + fft_sizes_ {other.fft_sizes()}, + kernels_ {Kernel {other.kernels()[0]}, Kernel {other.kernels()[1]}} +{ + const auto npix = static_cast(fft_sizes_[0]) * fft_sizes_[1]; + xt_.resize(npix); + checkCudaErrors(cudaMemcpy(xt_.data().get(), other.data(), + npix * sizeof(std::complex), cudaMemcpyHostToDevice)); +} + +template +NFFT2dResultView::NFFT2dResultView(const NFFT2dResult& result) : + sizes_{result.sizes_}, + fft_sizes_{result.fft_sizes_}, + pxt_{result.xt_.data().get()}, + kernels_{result.kernels_} + {} + +template +NFFT2dResult makeImageNFFT2d( + const Eigen::Ref>>& image, + const isce3::signal::NFFT2dParams& params, + bool pad_input) +{ + using isce3::fft::nextFastPower; + + auto rows_in = image.rows(); + auto cols_in = image.cols(); + using image_t = isce3::core::EArray2D>; + auto image_copy = image_t(0, 0); + + // Pointer to input image or padded/copied version so we can have fewer + // conditionals later. + // FIXME figure out how to do this with an Eigen type... + auto image_ptr = image.data(); + + // Need to copy if image is not contiguous row-major since we don't have + // high-level interface for strided FFTs. + bool need_copy = (image.innerStride() != 1) or (image.outerStride() != cols_in); + if (need_copy) { + image_copy.resize(rows_in, cols_in); + // assign later + } + + if (pad_input) { + auto padded_rows_in = nextFastPower(rows_in); + auto padded_cols_in = nextFastPower(cols_in); + if ((rows_in == padded_rows_in) and (cols_in == padded_cols_in)) { + // User asked for padding but we don't actually need it. + pad_input = false; + } else { + image_copy.resize(padded_rows_in, padded_cols_in); + image_copy.setZero(); + rows_in = padded_rows_in; + cols_in = padded_cols_in; + // assign later + } + } + + if (need_copy or pad_input) { + // This way NFFT2d::interp() coordinates are preserved, though user + // will be able to get some extra data. + image_copy.topLeftCorner(image.rows(), image.cols()) = image; + image_ptr = image_copy.data(); + } + + // now copy contiguous data to GPU + thrust::device_vector> d_image(rows_in * cols_in); + checkCudaErrors(cudaMemcpy(d_image.data().get(), image_ptr, + d_image.size() * sizeof(*image_ptr), cudaMemcpyHostToDevice)); + + // Use fft2 b/c planfft2d could modify inputs and we won't reuse it anyway. + using dims_t = typename NFFT2d::dims_t; + dims_t dims = { + static_cast(rows_in), + static_cast(cols_in)}; + // Since we've already made a copy, we can just FFT in-place. + isce3::cuda::fft::fft2d(d_image.data().get(), d_image.data().get(), + {dims[0], dims[1]}); + + // Calculate sizes for padded inverse transform. + dims_t dims_out = { + nextFastPower(static_cast(std::round(params.rows.s * dims[0]))), + nextFastPower(static_cast(std::round(params.cols.s * dims[1])))}; + + const dims_t m = {params.rows.m, params.cols.m}; + auto plan = NFFT2d(m, dims, dims_out); + return plan.transform_device(dims, {dims[1], 1}, d_image.data().get()); +} + +} + +template class isce3::cuda::signal::NFFT2d; +template class isce3::cuda::signal::NFFT2d; +template class isce3::cuda::signal::NFFT2dResult; +template class isce3::cuda::signal::NFFT2dResult; +template class isce3::cuda::signal::NFFT2dResultView; +template class isce3::cuda::signal::NFFT2dResultView; + +template isce3::cuda::signal::NFFT2dResult +isce3::cuda::signal::makeImageNFFT2d( + const Eigen::Ref>>& image, + const isce3::signal::NFFT2dParams& params, + bool pad_input); + +template isce3::cuda::signal::NFFT2dResult +isce3::cuda::signal::makeImageNFFT2d( + const Eigen::Ref>>& image, + const isce3::signal::NFFT2dParams& params, + bool pad_input); \ No newline at end of file diff --git a/cxx/isce3/cuda/signal/NFFT2d.h b/cxx/isce3/cuda/signal/NFFT2d.h new file mode 100644 index 000000000..024fce584 --- /dev/null +++ b/cxx/isce3/cuda/signal/NFFT2d.h @@ -0,0 +1,234 @@ +#pragma once + +#include "forward.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace isce3::cuda::signal { + +/** + * @brief GPU-accelerated Non-Uniform Fast Fourier Transform (NFFT) in 2D. + * + * This class performs a zero-padded, pre-filtered 2D inverse FFT to convert + * image-spectrum data into the time domain for interpolation. All transform + * and intermediate spectral data reside in GPU device memory. + */ +template +class NFFT2d { + public: + static constexpr int ndims = 2; + using dims_t = std::array; + + NFFT2d() = delete; + + /** + * @brief Construct a new NFFT2d object + * + * @param m Interpolator half-length along {rows, columns} + * @param sizes Image spectrum dimensions {rows, columns} + * @param fft_sizes Transform sizes along {rows, columns}. + * Usually larger than image size. + */ + NFFT2d(const dims_t& m, const dims_t& sizes, const dims_t& fft_sizes); + + /** + * @brief Transform image spectrum to time domain, input in host memory. + * + * @param sizes Image spectrum dimensions {rows, columns}. + * Must match dimensions provided in ctor. + * @param strides Strides (in pixels) along each dimension {rows, columns}. + * @param x Image spectrum. + * + * The input data is copied to device memory. The spectrum will be + * zero-padded, pre-filtered, and transformed to the time-domain. + * + * @return NFFT2dResult object containing the time-domain data. + */ + NFFT2dResult transform_host(const dims_t& sizes, + const dims_t& strides, const std::complex* x); + + /** + * @brief Transform image spectrum to time domain, input on device. + * + * @param sizes Image spectrum dimensions {rows, columns}. + * Must match dimensions provided in ctor. + * @param strides Strides (in pixels) along each dimension {rows, columns}. + * @param x Image spectrum (device pointer). + * + * The spectrum will be zero-padded, pre-filtered, and transformed + * to the time-domain entirely on the GPU. + * + * @return NFFT2dResult object containing the time-domain data. + */ + NFFT2dResult transform_device(const dims_t& sizes, + const dims_t& strides, const thrust::complex* x); + + /** Image spectrum dimensions */ + const dims_t& sizes() const { return sizes_; } + + /** Transform sizes */ + const dims_t& fft_sizes() const { return fft_sizes_; } + + /** Pointer to most recent spectral data (filtered and padded) */ + const thrust::complex* spectrum() const { return xf_.data().get(); } + + private: + dims_t m_, sizes_, fft_sizes_; + thrust::device_vector> xf_; + std::array, 2> weights_; + std::array, 2> kernels_; +}; + + +/** + * @brief Result of a GPU NFFT2d transform, holding time-domain data on the device. + * + * The time-domain data (xt_) is stored in GPU device memory as a + * thrust::device_vector. This object can be converted to the CPU + * counterpart (isce3::signal::NFFT2dResult) or wrapped in a + * NFFT2dResultView for GPU-side interpolation. + */ +template +class NFFT2dResult { + friend class NFFT2d; + friend class NFFT2dResultView; + +public: + using dims_t = typename NFFT2d::dims_t; + + NFFT2dResult() = delete; + + /** + * @brief Construct a NFFT2dResult. + * + * @param m Interpolator half-length along {rows, columns} + * @param sizes Image spectrum dimensions {rows, columns} + * @param fft_sizes Transform sizes along {rows, columns} + * @param kernels Interpolator kernels along {rows, columns} + * @param xt Optional device pointer to time-domain data of length + * fft_sizes[0] * fft_sizes[1]. If null, a zero-filled + * device buffer is allocated. + */ + NFFT2dResult(const dims_t& m, const dims_t& sizes, const dims_t& fft_sizes, + const std::array, 2>& kernels, + const thrust::complex* xt = nullptr) + : m_ {m}, sizes_ {sizes}, fft_sizes_ {fft_sizes}, kernels_ {kernels} + { + const auto n = static_cast(fft_sizes[0]) * fft_sizes[1]; + if (xt == nullptr) { + xt_.resize(n); + } else { + xt_.assign(xt, xt + n); + } + }; + + /** copy to host */ + explicit operator isce3::signal::NFFT2dResult() const; + + /** copy from host */ + NFFT2dResult(const isce3::signal::NFFT2dResult& other); + + /** Interpolator half-lengths along {rows, columns} */ + const dims_t& kernel_radii() const { return m_; } + /** Image spectrum dimensions */ + const dims_t& sizes() const { return sizes_; } + /** Transform sizes */ + const dims_t& fft_sizes() const { return fft_sizes_; } + /** Interpolator kernels along {rows, columns} */ + const auto& kernels() const { return kernels_; } + +private: + dims_t m_, sizes_, fft_sizes_; + std::array, 2> kernels_; + thrust::device_vector> xt_; +}; + + +/** + * @brief Lightweight view for GPU-side interpolation of NFFT2d results. + * + * This class wraps a NFFT2dResult and provides a CUDA device-callable + * interpolation method. It is intended for use within device kernels where + * multiple pixels need to be interpolated without transferring data back to + * the host. + */ +template +class NFFT2dResultView { +public: + static constexpr int ndims = 2; + using dims_t = std::array; + + NFFT2dResultView() = delete; + + /** Construct a lightweight view of an NFFT2dResult. */ + NFFT2dResultView(const NFFT2dResult& result); + + /** + * @brief Interpolate the image on the device. + * + * @param t Desired pixel location {row, column} + * Values should be in 0 <= t[i] < sizes()[i]. + * @param periodic Whether to use a periodic boundary condition. + * @return Interpolated value. + */ + CUDA_DEV inline + thrust::complex interp( + const std::array& t, bool periodic = true) const + { + constexpr int xdim = 1, ydim = 0; + + // scale time index to account for zero-padding of spectrum. + double x = t[xdim] * fft_sizes_[xdim] / sizes_[xdim]; + double y = t[ydim] * fft_sizes_[ydim] / sizes_[ydim]; + + return isce3::cuda::core::interp2d(kernels_[xdim], kernels_[ydim], pxt_, + fft_sizes_[xdim], /* stridex */ 1, fft_sizes_[ydim], + /* stridey */ fft_sizes_[xdim], x, y, periodic); + }; + + /** Transform sizes */ + CUDA_DEV + const dims_t& fft_sizes() const { return fft_sizes_; } + +private: + dims_t sizes_, fft_sizes_; + const thrust::complex* pxt_; + std::array, 2> kernels_; +}; + + +/** + * @brief Create an NFFT2dResult object for interpolating an image. + * + * @tparam T Format of real/imag pixel data, typically float or double + * @param image Input time-domain image. A temporary copy will be made if + * it is not row-major with a column stride of one. + * @param m Half-length of interpolator along {rows, columns} + * @param s Minimum factors (> 1) for frequency-domain zero-padding + * along {rows, columns}. Actual padding may be larger to + * achieve efficient inverse transform size. + * @param pad_input Whether to also zero-pad input data to an efficient + * forward transform size. Requires extra memory. + * + * @return NFFT2dResult object for interpolating the image. + */ +template +NFFT2dResult makeImageNFFT2d( + const Eigen::Ref>>& image, + const isce3::signal::NFFT2dParams& params = {}, + bool pad_input = false); + +} diff --git a/cxx/isce3/cuda/signal/forward.h b/cxx/isce3/cuda/signal/forward.h index 41f96e406..843e06346 100644 --- a/cxx/isce3/cuda/signal/forward.h +++ b/cxx/isce3/cuda/signal/forward.h @@ -9,5 +9,8 @@ namespace isce3 { namespace cuda { namespace signal { template class gpuLooks; template class gpuRangeFilter; template class gpuSignal; + template class NFFT2d; + template class NFFT2dResult; + template class NFFT2dResultView; }}} // namespace isce3::cuda::signal diff --git a/cxx/isce3/error/ErrorCode.cpp b/cxx/isce3/error/ErrorCode.cpp index d91ed4cb5..99efa65a4 100644 --- a/cxx/isce3/error/ErrorCode.cpp +++ b/cxx/isce3/error/ErrorCode.cpp @@ -30,6 +30,8 @@ std::string getErrorString(ErrorCode status) return "specified tolerance or number of iterations is invalid"; case ErrorCode::InvalidInterval: return "specified interval does not contain a solution"; + case ErrorCode::InvalidKernelSize: + return "kernel size exceeds compile-time limit"; } throw isce3::except::RuntimeError(ISCE_SRCINFO(), "unknown error code"); diff --git a/cxx/isce3/error/ErrorCode.h b/cxx/isce3/error/ErrorCode.h index 4c801610a..615fc03d4 100644 --- a/cxx/isce3/error/ErrorCode.h +++ b/cxx/isce3/error/ErrorCode.h @@ -18,6 +18,7 @@ enum class ErrorCode { NullDereference, InvalidTolerance, InvalidInterval, + InvalidKernelSize, }; /** Return a string describing the error code */ diff --git a/cxx/isce3/focus/Backproject.cpp b/cxx/isce3/focus/Backproject.cpp index fa94b2820..f6556b22c 100644 --- a/cxx/isce3/focus/Backproject.cpp +++ b/cxx/isce3/focus/Backproject.cpp @@ -5,14 +5,18 @@ #include #include #include +#include #include #include #include +#include +#include #include #include #include #include #include +#include #include #include @@ -23,6 +27,10 @@ using namespace isce3::geometry; using isce3::error::ErrorCode; using isce3::container::RadarGeometry; +using isce3::signal::NFFT2dResult; +using isce3::signal::NFFT2dParams; +using isce3::fft::planfft2d; +using isce3::fft::nextFastPower; namespace isce3 { namespace focus { @@ -211,5 +219,847 @@ backproject(std::complex* out, const RadarGeometry& out_geometry, return ErrorCode::Success; } + +static Vec3 vector_mean(const std::vector& vecs) +{ + Vec3 sum = {0, 0, 0}; + for (const auto& vec : vecs) { + sum += vec; + } + return sum * (1.0 / vecs.size()); +} + +double +getPolarAngleTimeConstant(const double fc, const double vs, + const double bandwidth, const double c) +{ + // Yegulalp, Eq. (11) + const auto fmax = fc + bandwidth / 2; + return c / (2 * fmax * vs); +} + + +std::tuple, std::vector> +setupPolarGridForPulses( + const RadarGeometry& in_geometry, + const Eigen::Ref& azimuth_time, + double range_bandwidth, + double azimuth_resolution, + double oversample_range, double oversample_azimuth, + int num_doppler_eval, std::optional pri) +{ + // Interpolate platform position & velocity at each pulse + const auto nt = azimuth_time.size(); + if (nt < 1) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Need at least one pulse to setup polar grid."); + } + std::vector pos(nt), vel(nt); + + for (auto i = decltype(nt){0}; i < nt; ++i) { + double t = azimuth_time[i]; + in_geometry.orbit().interpolate(&pos[i], &vel[i], t); + } + + // For the along-track axis we could fit a line to the positions, or use the + // dominant eigenvector of the position sample covariance. But the average + // velocity is probably about the same and simpler to compute. + Vec3 axis = vector_mean(vel); + const auto vs = axis.norm(); + axis *= 1.0 / vs; + + constexpr auto c = isce3::core::speed_of_light; + const auto fc = c / in_geometry.wavelength(); + const auto slant_range = in_geometry.slantRange(); + + // Our polar data structures use a constant Doppler centroid (DC) vs range. + // If we have some DC variation over the swath, we'll increase the Doppler + // bandwidth enough to accommodate it. Later we can mask out the pixels + // outside the desired azimuth band if desired. + // We will assume the DC is stable over the slow-time span of the pulses. + const auto + r0 = slant_range.first(), + r1 = slant_range.last(), + t0 = azimuth_time[0], + t1 = azimuth_time[nt - 1], + tmid = (t0 + t1) / 2, + dop2q = c / (fc * 2 * vs), + pri_ = pri.value_or((t1 - t0) / (nt - 1)); + + auto q0 = in_geometry.doppler().eval(tmid, r0) * dop2q; + auto q1 = q0; + for (int i = 1; i < num_doppler_eval; ++i) { + const auto ri = r0 + i * (r1 - r0) / (num_doppler_eval - 1); + const auto qi = in_geometry.doppler().eval(tmid, ri) * dop2q; + q0 = std::min(q0, qi); + q1 = std::max(q1, qi); + } + auto qmid = (q0 + q1) / 2; + auto qspan = (q1 - q0) + c / (fc * 2 * azimuth_resolution); + + // Use mean position as origin of polar grid. + Vec3 origin = vector_mean(pos); + + // Bistatic correction, roughly 22 m for NISAR-like geometry (many pulses). + // If neglected causes a noticeable spectral shift for short apertures + // that can mess up baseband interpolation. + const auto rmid = (r0 + r1) / 2; + const auto ds_dr = vs / (c - qmid * vs); + origin += rmid * ds_dr * axis; + + // Depends on range, so adjust aperture duration by variation in shift. + // This will cause a higher sample rate and hopefully avoid aliasing. + // Roughly 4 m for NISAR-like geometry (one pulse, almost negligible). + const auto duration = (t1 - t0 + pri_) + (r1 - r0) * ds_dr / vs; + + // Yegulalp, Eq. (11) and (12) + const auto tq = getPolarAngleTimeConstant(fc, vs, range_bandwidth, c); + auto dq = tq / (duration * oversample_azimuth); + auto dr = c / (2 * range_bandwidth * oversample_range); + + // Though inefficient, user might try to combine more pulses than are + // needed to achieve the desired azimuth resolution. For example, they + // might try to backproject all pulses from a stripmap radar in one shot. + const auto dq_min = azimuth_resolution / + (slant_range.last() * oversample_azimuth); + if (dq < dq_min) { + // TODO emit a warning? + dq = dq_min; + } + + int nr = 1 + static_cast(std::ceil((r1 - r0) / dr)); + int nq = 1 + static_cast(std::ceil(qspan / dq)); + + auto pgrid = PolarGrid{t0, t1 + pri_, + origin, axis, Linspace(r0, dr, nr), + Linspace(qmid - dq * (nq - 1) / 2, dq, nq), + in_geometry.lookSide()}; + + return {pgrid, pos, vel}; +} + + +std::tuple[]>, std::unique_ptr> +backprojectToPolarGrid( + const std::complex* in, const Linspace& in_slant_range, + const std::vector& pos, + const std::vector& vel, + const PolarGrid& out_grid, + const DEMInterpolator& dem, double fc, + const Kernel& kernel, DryTroposphereModel dry_tropo_model, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params) +{ + using isce3::geometry::detail::polar2geo_bracket; + + const auto nt = pos.size(); + if (vel.size() != nt) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "require same number of position and velocity vectors"); + } + + static constexpr double c = isce3::core::speed_of_light; + static constexpr auto nan = std::numeric_limits::quiet_NaN(); + + // check that dry_tropo_model is supported internally + if (not(dry_tropo_model == DryTroposphereModel::NoDelay or + dry_tropo_model == DryTroposphereModel::TSX)) { + + std::string errmsg = "unexpected dry troposphere model"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + const auto npix = static_cast(out_grid.length()) * out_grid.width(); + auto height = std::make_unique(npix); + auto out = std::make_unique[]>(npix); + + // range sampling window + double swst = 2. * in_slant_range.first() / c; + double dtau = 2. * in_slant_range.spacing() / c; + int nr = in_slant_range.size(); + Linspace sampling_window(swst, dtau, nr); + + // reference ellipsoid + int epsg = dem.epsgCode(); + const Ellipsoid ellipsoid = makeProjection(epsg)->ellipsoid(); + + // loop over targets in output grid + bool all_converged = true; +#pragma omp parallel for + for (int j = 0; j < out_grid.sin_squint.size(); ++j) { + const double + q = out_grid.sin_squint[j], + c = std::sqrt(1.0 - q * q); + for (int i = 0; i < out_grid.range.size(); ++i) { + + // Run polar2geo to get target position. + // Only need LLH if dumping height or using TSX atmosphere model, + // but just compute it unconditionally. + Vec3 x, llh; + { + const double r = out_grid.range[i]; + double look_angle; + + const auto status = polar2geo_bracket(&x, &look_angle, + out_grid.origin, out_grid.axis, r, q, c, dem, ellipsoid, + out_grid.look_side, r2g_params); + + llh = ellipsoid.xyzToLonLat(x); + height[j * out_grid.width() + i] = llh[2]; + + if (status != isce3::error::ErrorCode::Success) { + all_converged = false; + out[j * out_grid.width() + i] = {nan, nan}; + height[j * out_grid.width() + i] = nan; + continue; + } + } + + // estimate dry troposphere delay + double tau_atm = 0.; + if (dry_tropo_model == DryTroposphereModel::TSX) { + tau_atm = dryTropoDelayTSX(out_grid.origin, llh, ellipsoid); + } + + // TODO range-dependent Doppler mask? + int kstart = 0, kstop = static_cast(nt); + + // integrate pulses + out[j * out_grid.width() + i] = + sumCoherent(in, sampling_window, pos, vel, x, fc, tau_atm, + kernel, kstart, kstop); + } + } + + // baseband + const double kw = 4 * M_PI / (c / fc); + #pragma omp parallel for + for (int i = 0; i < out_grid.range.size(); ++i) { + const double phi = -kw * out_grid.range[i]; + const auto phasor = std::complex(std::cos(phi), std::sin(phi)); + for (int j = 0; j < out_grid.sin_squint.size(); ++j) { + out[j * out_grid.width() + i] *= phasor; + } + } + + auto status = + all_converged ? ErrorCode::Success : ErrorCode::FailedToConverge; + return std::make_tuple(status, std::move(out), std::move(height)); +} + +PolarGrid +mergePolarGrids(const std::vector& grids, + const DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const std::optional& dq_min, + const std::optional& tq) +{ + if (grids.size() <= 0) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "can't find common grid among empty list"); + } else if (grids.size() == 1) { + return grids[0]; + } + + // reference ellipsoid + Ellipsoid ellipsoid = makeProjection(dem.epsgCode())->ellipsoid(); + + // Compute a bunch of stats with a first pass over the data. + // Average origin and axis, weighted by aperture duration. + Vec3 origin{0, 0, 0}, axis{0, 0, 0}; + // Inferred dimensionless Doppler spacing time constant + double tq_inferred = 0.0; + // Min range spacing (in case different among grids) + auto dr = grids[0].range.spacing(); + // Need total aperture size and sum of subaperture sizes. + // These are not equal if there are gaps or overlap between subapertures. + auto t_min = grids[0].aztime_start; // assume start > end + auto t_max = grids[0].aztime_end; // assume start > end + double sum_durations = 0; + const auto look_side = grids[0].look_side; + + for (const auto& grid : grids) { + const auto duration = grid.aztime_end - grid.aztime_start; + sum_durations += duration; + t_min = std::min(t_min, grid.aztime_start); // assume start > end + t_max = std::max(t_max, grid.aztime_end); // assume start > end + dr = std::min(dr, grid.range.spacing()); + origin += duration * grid.origin; + axis += duration * grid.axis; + tq_inferred += duration * (grid.sin_squint.spacing() * duration); + if (grid.look_side != look_side) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "inconsistent look_side among input polar grids"); + } + } + origin *= 1.0 / sum_durations; + axis *= 1.0 / axis.norm(); + tq_inferred /= sum_durations; + + // In general, figuring out the required Doppler spacing is pretty complex. + // You'd want to figure out the Doppler bandwidth observed by all targets + // across all grids, maxing out around the azimuth resolution. + // For now let's just just be conservative and increase it linearly. + auto dq = tq.value_or(tq_inferred) / (t_max - t_min); + + // But the user can override this. + if (dq_min) { + dq = std::max(dq_min.value(), dq); + } + + // Compute range & Doppler bounds of new grid using corners of each input. + // Use lambda to avoid copy/paste. + using isce3::geometry::detail::polar2polar_bracket; + auto polar2polar = [&](const PolarGrid& grid, double r, double ssq) { + auto csq = std::sqrt(1.0 - ssq * ssq); + double r_out, ssq_out; + auto ec = polar2polar_bracket(&ssq_out, &r_out, ssq, csq, r, + grid.origin, grid.axis, origin, axis, dem, ellipsoid, look_side, + r2g_params); + if (ec != ErrorCode::Success) { + throw isce3::except::DomainError(ISCE_SRCINFO(), + "polar2polar failed with ErrorCode (" + + isce3::error::getErrorString(ec) + ") for point at r=" + + std::to_string(r) + " sin_squint=" + std::to_string(ssq)); + } + return std::make_tuple(r_out, ssq_out); + }; + + // NOTE Use grid _edges_ for determining extent. Okay to initialize with + // center, though. + auto [r_min, q_min] = polar2polar(grids[0], grids[0].range[0], + grids[0].sin_squint[0]); + auto r_max = r_min, q_max = q_min; + for (const auto& grid : grids) { + for (const auto& ri : grid.range.bounds()) { + for (const auto& qi : grid.sin_squint.bounds()) { + const auto [ro, qo] = polar2polar(grid, ri, qi); + r_min = std::min(r_min, ro); + r_max = std::max(r_max, ro); + q_min = std::min(q_min, qo); + q_max = std::max(q_max, qo); + } + } + } + + const int nr = static_cast(std::ceil((r_max - r_min) / dr)); + const int nq = static_cast(std::ceil((q_max - q_min) / dq)); + + // The ceil() means potentially extra coverage. We'll center it so there's + // equal padding on both sides of the interval. Note also that min/max are + // bin edges while we're specifying bin centers, hence (N-1) instead of N + // in the formulas. + const auto r0 = r_min - ((nr - 1) * dr - (r_max - r_min)) / 2; + const auto q0 = q_min - ((nq - 1) * dq - (q_max - q_min)) / 2; + return PolarGrid{t_min, t_max, origin, axis, + Linspace(r0, dr, nr), + Linspace(q0, dq, nq), + look_side}; +} + + +void mergePolarImages( + const std::vector& grids, + const std::vector*>& image_interpolators, + const PolarGrid& output_grid, + Eigen::Ref>> output_image, + const double fc, + const DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + int az_block_size) +{ + // check that output grid dimensions match buffer size + const auto m = output_grid.length(), n = output_grid.width(); + if ((m != output_image.rows()) or (n != output_image.cols())) { + std::string msg = "Dimensions of image grid (" + std::to_string(m) + + ", " + std::to_string(n) + ") do not match dimensions of image " + "buffer (" + std::to_string(output_image.rows()) + ", " + + std::to_string(output_image.cols()) + ")"; + throw isce3::except::LengthError(ISCE_SRCINFO(), msg); + } + + // check that we have a grid for each input image + const auto num_images = image_interpolators.size(); + if (grids.size() != num_images) { + std::string msg = "Size mismatch: got " + std::to_string(num_images) + + " sub images but " + std::to_string(grids.size()) + " grids"; + throw isce3::except::LengthError(ISCE_SRCINFO(), msg); + } + + // check look directions for consistency + const auto look_side = output_grid.look_side; + for (const auto& grid : grids) { + if (grid.look_side != look_side) { + std::string msg = "Output grid look direction does not match " + "input grid look direction"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), msg); + } + } + + // Check block size and allocate scratch space. + if (az_block_size <= 0) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "azimuth block size must be positive"); + } + az_block_size = std::min(az_block_size, output_grid.sin_squint.size()); + + auto block_positions = isce3::core::EArray2D(); + block_positions.resize(az_block_size, output_grid.width()); + + // reference ellipsoid + Ellipsoid ellipsoid = makeProjection(dem.epsgCode())->ellipsoid(); + + // wavenumber + const double kw = 4 * M_PI * fc / isce3::core::speed_of_light; + + using isce3::geometry::detail::polar2polar_bracket; + using isce3::geometry::geo2polar; + + // loop over output blocks + auto n_blocks = (m + az_block_size - 1) / az_block_size; + for (auto i_block = decltype(n_blocks){0}; i_block < n_blocks; ++i_block) { + auto i_row0 = i_block * az_block_size; + auto i_row1 = std::min(i_row0 + az_block_size, m); + + // Compute output pixel 3D locations + using isce3::geometry::detail::polar2geo_bracket; + #pragma omp parallel for collapse(2) + for (auto i_row = i_row0; i_row < i_row1; ++i_row) { + for (auto j = decltype(n){0}; j < n; ++j) { + auto i = i_row - i_row0; + double look_angle; + const auto ssq = output_grid.sin_squint[i_row]; + const auto csq = std::sqrt(1.0 - ssq * ssq); + auto ec = polar2geo_bracket(&block_positions(i, j), &look_angle, + output_grid.origin, output_grid.axis, output_grid.range[j], + ssq, csq, dem, ellipsoid, output_grid.look_side, r2g_params); + if (ec != ErrorCode::Success) { + throw isce3::except::DomainError(ISCE_SRCINFO(), + "polar2geo failed with ErrorCode (" + + isce3::error::getErrorString(ec) + ") for point at r=" + + std::to_string(output_grid.range[j]) + " sin_squint=" + + std::to_string(ssq)); + } // err + } // columns + } // rows + + // loop over input images + for (auto i_img = decltype(num_images){0}; i_img < num_images; ++i_img) { + const auto& input_grid = grids[i_img]; + const auto* nfft = image_interpolators[i_img]; + const auto npix = static_cast(n) * (i_row1 - i_row0); + auto ec = accumulatePolarImageToGeoPoints( + output_image.row(i_row0).data(), + block_positions.data(), npix, input_grid, *nfft, kw); + if (ec != ErrorCode::Success) { + throw isce3::except::RuntimeError(ISCE_SRCINFO(), + "projectPolarToGeo failed with ErrorCode (" + + isce3::error::getErrorString(ec) + ")"); + } // error + } // images + } // blocks + + // Baseband. Note that we could do this at the same time as the + // reprojection but it'd require a fair bit of copy/paste. + Eigen::VectorXcf phasors(n); + #pragma omp parallel for + for (auto j = decltype(n){0}; j < n; ++j) { + const double arg = -kw * output_grid.range[j]; + phasors(j) = std::complex(std::cos(arg), std::sin(arg)); + } + #pragma omp parallel for collapse(2) + for (auto i = decltype(m){0}; i < m; ++i) { + for (auto j = decltype(n){0}; j < n; ++j) { + output_image(i, j) *= phasors(j); + } // columns + } // rows +} + + +// For now structure like backproject() with inner loop on target. +// Might make more sense to project one image at a time instead. +ErrorCode +accumulatePolarImagesToRadarGrid(std::complex* out, + const RadarGeometry& out_geometry, + const isce3::core::Orbit& in_orbit, + const isce3::core::LUT2d& in_doppler, + const std::vector& grids, + const std::vector*>& image_interpolators, + const DEMInterpolator& dem, double fc, double ds, + DryTroposphereModel dry_tropo_model, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params, + float* height) +{ + static constexpr double c = isce3::core::speed_of_light; + static constexpr auto nan = std::numeric_limits::quiet_NaN(); + + // check that dry_tropo_model is supported internally + if (not(dry_tropo_model == DryTroposphereModel::NoDelay or + dry_tropo_model == DryTroposphereModel::TSX)) { + + std::string errmsg = "unexpected dry troposphere model"; + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + // will search sorted intervals to figure out active sub images per target + auto starts = std::vector(grids.size()); + std::transform(grids.begin(), grids.end(), starts.begin(), + [](const PolarGrid& grid) { return grid.aztime_start; }); + auto ends = std::vector(grids.size()); + std::transform(grids.begin(), grids.end(), ends.begin(), + [](const PolarGrid& grid) { return grid.aztime_end; }); + + // get input & output radar grid azimuth time & slant range + Linspace out_azimuth_time = out_geometry.sensingTime(); + Linspace out_slant_range = out_geometry.slantRange(); + + // reference ellipsoid + int epsg = dem.epsgCode(); + Ellipsoid ellipsoid = makeProjection(epsg)->ellipsoid(); + + // carrier wavelength + const double wvl = c / fc; + const double kw = 4 * M_PI / wvl; + + const size_t nout = out_geometry.gridLength() * out_geometry.gridWidth(); + std::vector x(nout); + std::vector tstart(nout), tend(nout), dr_atm(nout); + + // loop over targets in output grid + bool all_converged = true; + #pragma omp parallel for + for (size_t iflat = 0; iflat < nout; ++iflat) { + const size_t j = iflat / out_slant_range.size(); + const size_t i = iflat % out_slant_range.size(); + + // Run rdr2geo using orbit and Doppler associated with output grid + // to get target position. Only need LLH if dumping height or + // using TSX atmosphere model, but just compute it unconditionally. + Vec3 llh; + { + double t = out_azimuth_time[j]; + double r = out_slant_range[i]; + double fD = out_geometry.doppler().eval(t, r); + + const int converged = rdr2geo_bracket(t, r, fD, + out_geometry.orbit(), dem, x[iflat], wvl, + out_geometry.lookSide(), r2g_params.tol_height, + r2g_params.look_min, r2g_params.look_max); + + llh = ellipsoid.xyzToLonLat(x[iflat]); + + if (height != nullptr) { + height[iflat] = llh[2]; + } + if (not converged) { + all_converged = false; + out[iflat] = {nan, nan}; + if (height != nullptr) { + height[iflat] = nan; + } + continue; + } + } + + // run geo2rdr to estimate the center of the coherent processing + // window for the target + double t, r; + { + auto converged = + geo2rdr_bracket(x[iflat], in_orbit, + in_doppler, t, r, wvl, + out_geometry.lookSide(), // assumed same side + g2r_params.tol_aztime, + g2r_params.time_start, g2r_params.time_end); + + if (not converged) { + all_converged = false; + out[iflat] = {nan, nan}; + continue; + } + } + + // get platform position and velocity at center of CPI + Vec3 p, v; + in_orbit.interpolate(&p, &v, t); + + // estimate synthetic aperture length required to achieve the + // desired azimuth resolution + double l = wvl * r * (p.norm() / x[iflat].norm()) / (2. * ds); + + // approximate CPI duration (assuming constant platform velocity) + double cpi = l / v.norm(); + + // get coherent integration bounds (pulse indices) + tstart[iflat] = t - cpi / 2; + tend[iflat] = tstart[iflat] + cpi; + + // Calculate dry troposphere delay (in units of one-way range). + if (dry_tropo_model == DryTroposphereModel::TSX) { + dr_atm[iflat] = dryTropoDelayTSX(p, llh, ellipsoid) * c / 2.; + } + // else zero-initialized by vector ctor + } + + // std::vector unsuitable due to bit packing optimizations + Eigen::Array mask(nout); + + // TODO reduce tstart & tend + // TODO check this O(log(n)) algorithm + //const auto kstart = std::distance(ends.begin(), + // std::lower_bound(ends.begin(), ends.end(), tstart)); + //const auto kstop = std::distance(starts.begin(), + // std::upper_bound(starts.start(), starts.end(), tstart + cpi)); + const auto num_images = image_interpolators.size(); + const decltype(num_images) kstart = 0, kstop = num_images; + + for (auto k = kstart; k < kstop; ++k) { + // check if we need to replan FFTs + const auto& grid = grids[k]; + const auto* nfft = image_interpolators[k]; + makeSubApertureMask(grid.aztime_start, grid.aztime_end, + nout, tstart.data(), tend.data(), mask.data()); + accumulatePolarImageToGeoPoints(out, x.data(), nout, grid, *nfft, kw, + mask.data(), dr_atm.data()); + } + + if (not all_converged) { + return ErrorCode::FailedToConverge; + } + return ErrorCode::Success; +} + +void +makeSubApertureMask( + const double subaperture_start, const double subaperture_end, + const size_t n, + const double* pixel_start, + const double* pixel_end, + bool* mask) +{ + #pragma omp parallel for + for (auto i = decltype(n){0}; i < n; ++i) { + mask[i] = (subaperture_end > pixel_start[i]) + and (subaperture_start < pixel_end[i]); + } +} + +ErrorCode +accumulatePolarImageToGeoPoints( + std::complex* image, + const Vec3* xyz, + const size_t n, + const PolarGrid& grid, + const NFFT2dResult& nfft, + const double kw, + const std::optional& mask, + const std::optional& dr_atm) +{ + #pragma omp parallel for + for (size_t i= 0; i < n; ++i) { + if (mask.has_value() and not mask.value()[i]) { + continue; + } + // compute target location in polar grid + double sin_squint, range; + geo2polar(&sin_squint, &range, xyz[i], grid.origin, grid.axis); + if (dr_atm.has_value()) { + range += dr_atm.value()[i]; + } + // convert to image index + const double ix = (range - grid.range.first()) / grid.range.spacing(), + iy = (sin_squint - grid.sin_squint.first()) / grid.sin_squint.spacing(); + // interpolate baseband data + const auto z = nfft.interp({iy, ix}, /* periodic */ false); + // compensate phase and sum contribution + const double phase = kw * range; + image[i] += + z * std::complex(std::cos(phase), std::sin(phase)); + } + return ErrorCode::Success; +} + +std::tuple +findPolarGridBoundingBoxInRadarCoord( + const PolarGrid& polar_grid, + const Orbit& orbit, + const LUT2d& doppler, + const double wavelength, + const LookSide lookside, + const DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params, + const int nextra) +{ + using isce3::geometry::detail::polar2geo_bracket; + if (nextra < 0) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "specified negative number of extra points"); + } + + // get (angle, range) points along perimeter of polar grid + const int n = 4 * (1 + nextra); + int nwritten = 0; + std::vector> points(n); + for (int i = 0; i <= nextra; ++i) { + const auto q = polar_grid.sin_squint.first(); + const auto dr = (polar_grid.range.last() - polar_grid.range.first()) / + (1 + nextra); + const auto r = polar_grid.range.first() + i * dr; + points[nwritten++] = {q, r}; + } + for (int i = 0; i <= nextra; ++i) { + const auto r = polar_grid.range.last(); + const auto dq = (polar_grid.sin_squint.last() - polar_grid.sin_squint.first()) / + (1 + nextra); + const auto q = polar_grid.sin_squint.first() + i * dq; + points[nwritten++] = {q, r}; + } + for (int i = 0; i <= nextra; ++i) { + const auto q = polar_grid.sin_squint.last(); + const auto dr = (polar_grid.range.last() - polar_grid.range.first()) / + (1 + nextra); + const auto r = polar_grid.range.last() - i * dr; + points[nwritten++] = {q, r}; + } + for (int i = 0; i <= nextra; ++i) { + const auto r = polar_grid.range.first(); + const auto dq = (polar_grid.sin_squint.last() - polar_grid.sin_squint.first()) / + (1 + nextra); + const auto q = polar_grid.sin_squint.last() - i * dq; + points[nwritten++] = {q, r}; + } + assert(nwritten == n); + + int epsg = dem.epsgCode(); + Ellipsoid ellipsoid = makeProjection(epsg)->ellipsoid(); + auto status = ErrorCode::Success; + + #pragma omp parallel for + for (int i = 0; i < n; ++i) { + // read polar coordinate + const double ssq = points[i][0]; + const double rin = points[i][1]; + // compute cos from sin assuming abs(squint) < 90 deg + const double csq = std::sqrt(1.0 - ssq * ssq); + // convert to xyz + Vec3 xyz; + double lookangle; + auto err = polar2geo_bracket(&xyz, &lookangle, polar_grid.origin, + polar_grid.axis, rin, ssq, csq, dem, ellipsoid, + lookside, r2g_params); + if (err != ErrorCode::Success) { + status = err; + } + // convert to stripmap radar coordinates + double tout, rout; + int success = geo2rdr_bracket(xyz, orbit, + doppler, tout, rout, wavelength, + lookside, g2r_params.tol_aztime, g2r_params.time_start, + g2r_params.time_end); + if (!success) { + status = ErrorCode::FailedToConverge; + } + // write back + points[i] = {tout, rout}; + } + + // find extrema + double tmin, tmax, rmin, rmax; + tmin = tmax = points[0][0]; + rmin = rmax = points[0][1]; + for (int i = 1; i < n; ++i) { + const double t = points[i][0], r = points[i][1]; + if (t > tmax) tmax = t; + if (t < tmin) tmin = t; + if (r > rmax) rmax = r; + if (r < rmin) rmin = r; + } + + return std::make_tuple(tmin, tmax, rmin, rmax, status); +} + +std::tuple +findPolarGridBoundingBoxInRadarGrid( + const PolarGrid& polar_grid, + const RadarGeometry& radar_geom, + const DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params, + const int nextra) +{ + auto [tmin, tmax, rmin, rmax, status] = + findPolarGridBoundingBoxInRadarCoord(polar_grid, radar_geom.orbit(), + radar_geom.doppler(), radar_geom.wavelength(), + radar_geom.lookSide(), dem, r2g_params, g2r_params, nextra); + + // too much typing + const auto t0 = radar_geom.sensingTime().first(); + const auto dt = radar_geom.sensingTime().spacing(); + const auto r0 = radar_geom.slantRange().first(); + const auto dr = radar_geom.slantRange().spacing(); + const int m = static_cast(radar_geom.gridLength()); + const int n = static_cast(radar_geom.gridWidth()); + + // convert extrema to indices in radar grid + int i0, j0, i1, j1; + i0 = static_cast(std::floor((tmin - t0) / dt)); + i1 = static_cast(std::ceil((tmax - t0) / dt)); + j0 = static_cast(std::floor((rmin - r0) / dr)); + j1 = static_cast(std::ceil((rmax - r0) / dr)); + + // return empty grid if non-overlapping + if ((i1 < 0) or (i0 >= m) or (j1 < 0) or (j0 >= n)) { + return std::make_tuple(0, 0, 0, 0, status); + } + + // otherwise clamp to grid bounds + i0 = std::max(0, std::min(i0, m - 1)); + i1 = std::max(0, std::min(i1, m)); + j0 = std::max(0, std::min(j0, n - 1)); + j1 = std::max(0, std::min(j1, n)); + + return std::make_tuple(i0, i1, j0, j1, status); +} + +std::tuple, ErrorCode> +computeRadarGridGeoPoints( + const RadarGeometry& geom, + const DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params) +{ + const size_t n = geom.gridLength() * geom.gridWidth(); + std::vector points(n); + auto status = computeRadarGridGeoPoints(points.data(), geom, dem, r2g_params); + return std::make_tuple(points, status); +} + +ErrorCode +computeRadarGridGeoPoints( + Vec3* points, + const RadarGeometry& geom, + const DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params) +{ + const size_t n = geom.gridLength() * geom.gridWidth(); + ErrorCode status = ErrorCode::Success; + #pragma omp parallel for + for (size_t k = 0; k < n; ++k) { + const int i = static_cast(k / geom.gridWidth()); + const int j = static_cast(k % geom.gridWidth()); + const double t = geom.sensingTime()[i]; + const double r = geom.slantRange()[j]; + const double fd = geom.doppler().eval(t, r); + const int success = isce3::geometry::rdr2geo_bracket(t, r, fd, + geom.orbit(), dem, points[k], geom.wavelength(), geom.lookSide(), + r2g_params.tol_height, r2g_params.look_min, r2g_params.look_max); + if (!success) { + // race condition okay since always pushing the same value + status = ErrorCode::FailedToConverge; + } + } + return status; +} + } // namespace focus } // namespace isce3 diff --git a/cxx/isce3/focus/Backproject.h b/cxx/isce3/focus/Backproject.h index d41d65fe6..2aec2f3fa 100644 --- a/cxx/isce3/focus/Backproject.h +++ b/cxx/isce3/focus/Backproject.h @@ -3,12 +3,16 @@ #include #include #include +#include #include +#include +#include #include #include #include +#include #include "DryTroposphereModel.h" @@ -46,5 +50,439 @@ backproject(std::complex* out, const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params = {}, float* height = nullptr); + +/** Structure describing the coordinate system of a polar image. */ +struct PolarGrid { + /** Start azimuth time of the synthetic aperture (s) */ + double aztime_start; + + /** End azimuth time of the synthetic aperture (s); one PRI past the last pulse */ + double aztime_end; + + /** ECEF origin of the polar coordinate system (m). + * + * This is the reference point from which polar coordinates are defined, + * typically computed as the mean platform position over the aperture. + */ + isce3::core::Vec3 origin; + + /** Unit vector along the azimuth (along-track) axis of the polar grid. + * + * This is the normalized mean platform velocity, used to define the + * azimuth direction of the polar coordinate system. + */ + isce3::core::Vec3 axis; + + /** Slant range grid (m) */ + isce3::core::Linspace range; + + /** Sine of the squint angle (dimensionless Doppler) grid */ + isce3::core::Linspace sin_squint; + + /** Side looking direction (left or right of flight track) */ + isce3::core::LookSide look_side; + + PolarGrid() = delete; + + CUDA_HOSTDEV auto width() const { return range.size(); } + CUDA_HOSTDEV auto length() const { return sin_squint.size(); } + + PolarGrid offsetAndResize(int q_off, int r_off, int nq, int nr) const + { + using LS = isce3::core::Linspace; + const auto q = LS(sin_squint[q_off], sin_squint.spacing(), nq); + const auto r = LS(range[r_off], range.spacing(), nr); + return PolarGrid { + aztime_start, aztime_end, origin, axis, r, q, look_side}; + } +}; + +/** + * @brief Backproject range-compressed signal into a polar grid. + * + * Performs time-domain backprojection of range-compressed SAR signal data + * onto a polar coordinate grid. For each output pixel in the polar grid, + * the signal is resampled from the input data by accumulating contributions + * from all pulses according to the instantaneous slant range, with phase + * compensation for motion and the troposphere. + * + * @param[in] in Input range-compressed signal data + * @param[in] in_slant_range Slant range grid of the input data (m) + * @param[in] pos Platform position vectors at each pulse + * (ECEF, m) + * @param[in] vel Platform velocity vectors at each pulse + * (ECEF, m/s) + * @param[in] out_grid Target polar grid to backproject onto + * @param[in] dem DEM + * @param[in] fc Center frequency (Hz) + * @param[in] kernel 1-D interpolation kernel + * @param[in] dry_tropo_model Dry troposphere path delay model + * @param[in] r2g_params rdr2geo configuration parameters + * + * @returns A tuple containing: + * - error code (non-zero if geometry fails to converge for + * any pixel, in which case values for those pixels are NaN) + * - focused signal data on the polar grid (size = + * out_grid.width() * out_grid.length()) + * - per-pixel height above the ellipsoid (m) + * + * @see PolarGrid for a description of the polar grid coordinate system. + */ +std::tuple[]>, + std::unique_ptr> +backprojectToPolarGrid(const std::complex* in, + const isce3::core::Linspace& in_slant_range, + const std::vector& pos, + const std::vector& vel, + const PolarGrid& out_grid, + const isce3::geometry::DEMInterpolator& dem, + double fc, + const isce3::core::Kernel& kernel, + DryTroposphereModel dry_tropo_model, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params = {}); + +/** + * @brief Accumulate polar grid images onto an output radar geometry grid. + * + * Combines multiple subaperture polar grid images together onto a stripmap + * radar geometry grid. For each pixel in the output grid, the target + * position is computed via rdr2geo, the corresponding coherent processing + * interval is determined via geo2rdr, and the polar image data is + * accumulated using NFFT-based interpolation. + * + * The caller is responsible for allocating the output arrays to the + * appropriate size (out_geometry.gridLength() * out_geometry.gridWidth()). + * + * @param[out] out Accumulated focused signal data + * @param[in] out_geometry Target output grid, orbit, and Doppler + * @param[in] in_orbit Input data orbit + * @param[in] in_doppler Input data Doppler centroid LUT + * @param[in] grids List of subaperture polar grids + * @param[in] image_interpolators NFFT interpolators for each polar grid + * @param[in] dem Digital elevation model (DEM) + * @param[in] fc Center frequency (Hz) + * @param[in] ds Desired azimuth resolution (m) + * @param[in] dry_tropo_model Dry troposphere path delay model + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * @param[in] g2r_params geo2rdr_bracket configuration parameters + * @param[out] height Height of each pixel (m) above the + * ellipsoid (optional, may be nullptr) + * + * @returns Non-zero error code if rdr2geo or geo2rdr fails to converge + * for any pixel, and the values for these pixels are set to NaN. + */ +isce3::error::ErrorCode +accumulatePolarImagesToRadarGrid(std::complex* out, + const isce3::container::RadarGeometry& out_geometry, + const isce3::core::Orbit& in_orbit, + const isce3::core::LUT2d& in_doppler, + const std::vector& grids, + const std::vector*>& image_interpolators, + const isce3::geometry::DEMInterpolator& dem, double fc, double ds, + const DryTroposphereModel dry_tropo_model = DryTroposphereModel::TSX, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params = {}, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params = {}, + float* height = nullptr); + + + +/** + * @brief Get the time constant associated with polar angle spacing + * + * @param fc Radar center frequency, Hz + * @param vs Satellite velocity (along azimuth axis), m/s + * @param bandwidth Radar bandwidth, Hz (defaults to zero, e.g., narrow band) + * @param c Speed of light, m/s (defaults to vacuum sol) + * @return Time constant $T_q$, s + * + * This time constant is used to determine the sampling requirement for the + * sine of the squint angle (dimensionless Doppler) + * $$ q = \frac{\vec{v}}{v} \cdot \hat{l} $$ + * where $\vec{v}$ is the velocity and $\hat{l}$ is the line-of-sight direction. + * Specifically, the Nyquist criterion is + * $$ \Delta q \leq \frac{T_q}{T_{sa}} $$ + * where $T_{sa}$ is the time duration of the synthetic aperture. + * + * Helps implement equation (11) in @cite yegulalp2013 + */ +double +getPolarAngleTimeConstant(const double fc, const double vs, + const double bandwidth = 0.0, const double c = isce3::core::speed_of_light); + +/** + * @brief Set up a polar grid for a group of pulses. + * + * Constructs a PolarGrid that covers the synthetic aperture defined by + * the given azimuth times, along with the platform position and velocity + * vectors interpolated at each pulse. The grid spacing in the sine of + * the squint angle and range is determined by the Doppler bandwidth, + * range bandwidth, and desired azimuth resolution. + * + * @param[in] in_geometry Input data grid, orbit, & doppler + * @param[in] azimuth_time Azimuth times at which to set up the grid (s) + * @param[in] range_bandwidth Radar range bandwidth (Hz) + * @param[in] azimuth_resolution Desired azimuth resolution (m) + * @param[in] oversample_range Range oversampling factor + * (defaults to 1.2) + * @param[in] oversample_azimuth Azimuth oversampling factor + * (defaults to 1.2) + * @param[in] num_doppler_eval Number of range locations to evaluate + * Doppler centroid for bandwidth estimation + * (defaults to 2) + * @param[in] pri Pulse repetition interval (s); if not + * provided, it is inferred from azimuth_time + * + * @returns A tuple containing: + * - the constructed PolarGrid + * - platform position vectors at each pulse (ECEF, m) + * - platform velocity vectors at each pulse (ECEF, m/s) + * + * @see PolarGrid for a description of the polar grid coordinate system. + */ +std::tuple, std::vector> +setupPolarGridForPulses( + const isce3::container::RadarGeometry& in_geometry, + const Eigen::Ref& azimuth_time, + double range_bandwidth, + double azimuth_resolution, + double oversample_range = 1.2, double oversample_azimuth = 1.2, + int num_doppler_eval = 2, std::optional pri = std::nullopt); + +/** + * @brief Create polar grid capable of sampling data from all input grids. + * + * @param grids List of subaperture grids. + * @param dem Digital elevation model reporting height (m) above the + * ellispoid associated with its CRS. + * @param r2g_params Root finding parameters for radar2geo + * @param dq_min Minimum allowed dimensionless Doppler spacing. + * Necessary for stripmap processing large subapertures. + * @param tq Time constant for dimensionless Doppler spacing. + * If not provided it will be inferred from input grids. + */ +PolarGrid +mergePolarGrids(const std::vector& grids, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params = {}, + const std::optional& dq_min = {}, + const std::optional& tq = {}); + +/** + * @brief Merge subaperture polar grid images into a single output grid. + * + * Combines multiple subaperture polar grid images onto a merged output + * polar grid. For each pixel in the output grid, the 3D target position + * is computed via polar2geo, and the input image data is accumulated + * via NFFT-based interpolation. The output image is expected to be + * zero-initialized by the caller. + * + * @param[in] grids List of subaperture input polar grids + * @param[in] image_interpolators NFFT interpolators for each input grid + * @param[in] output_grid Merged output polar grid + * @param[out] output_image Accumulated output image (must be + * zero-initialized); dimensions must + * match output_grid + * @param[in] fc Center frequency (Hz) + * @param[in] dem Digital elevation model (DEM) + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * @param[in] az_block_size Number of azimuth rows to process at + * a time (defaults to 1024) + * + * @throws isce3::except::LengthError if output image dimensions or + * grid/interpolator counts are inconsistent + * @throws isce3::except::InvalidArgument if look directions are + * inconsistent or az_block_size is negative + * @throws isce3::except::DomainError if polar2geo fails to converge + * @throws isce3::except::RuntimeError if NFFT interpolation fails + */ +void mergePolarImages( + const std::vector& grids, + const std::vector*>& image_interpolators, + const PolarGrid& output_grid, + Eigen::Ref>> output_image, + const double fc, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params = {}, + int az_block_size = 1024); + +/** + * @brief Compute a mask of pixels overlapping a subaperture. + * + * For each pixel, determines whether the pixel's coherent processing + * interval [pixel_start, pixel_end] overlaps with the subaperture + * [subaperture_start, subaperture_end]. + * + * @param[in] subaperture_start Start time of the subaperture (s) + * @param[in] subaperture_end End time of the subaperture (s) + * @param[in] n Number of pixels + * @param[in] pixel_start Start time of each pixel's CPI (s) + * @param[in] pixel_end End time of each pixel's CPI (s) + * @param[out] mask Output mask; true if the pixel + * overlaps the subaperture + */ +void +makeSubApertureMask( + const double subaperture_start, const double subaperture_end, + const size_t n, + const double* pixel_start, + const double* pixel_end, + bool* mask); + +/** + * @brief Interpolate a polar grid image to given XYZ positions. + * + * Accumulates (adds) contributions from a polar grid image into an + * output complex signal array at specified 3D positions. For each + * position, the target location in the polar grid is computed via + * geo2polar, and the image is interpolated using NFFT. The phase is + * compensated by the wavenumber-range product kw * range. + * + * The output image is accumulated, so it must be zero-initialized + * by the caller before calling this function. + * + * @param[out] image Output complex signal data (accumulates, + * so caller must init to zero) + * @param[in] xyz Target 3D positions (ECEF, m) + * @param[in] n Number of target positions + * @param[in] grid Polar grid containing the image data + * @param[in] nfft NFFT interpolator for the polar grid + * @param[in] kw Wavenumber (rad/m); equals 4*pi*fc/c + * @param[in] mask Optional mask; pixels with false are skipped + * @param[in] dr_atm Optional atmospheric path delay correction + * (m), added to the range before interpolation + * + * @returns Error code indicating success or failure of the interpolation + */ +isce3::error::ErrorCode +accumulatePolarImageToGeoPoints( + std::complex* image, + const isce3::core::Vec3* xyz, + const size_t n, + const PolarGrid& grid, + const isce3::signal::NFFT2dResult& nfft, + const double kw, + const std::optional& mask = std::nullopt, + const std::optional& dr_atm = std::nullopt); + +/** + * @brief Find the bounding box of a polar grid in radar coordinates. + * + * Converts the perimeter of the polar grid to stripmap radar coordinates + * (azimuth time, slant range) and finds the minimum bounding box. + * The perimeter is sampled with nextra+1 points along each edge + * (4 edges), yielding 4*(nextra+1) total perimeter points. + * + * @param[in] polar_grid Input polar grid + * @param[in] orbit Orbit used to convert XYZ to radar coords + * @param[in] doppler Doppler centroid LUT + * @param[in] wavelength Radar wavelength (m) + * @param[in] lookside Look side (left or right) + * @param[in] dem Digital elevation model (DEM) + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * @param[in] g2r_params geo2rdr_bracket configuration parameters + * @param[in] nextra Number of extra perimeter points per edge + * (defaults to 0) + * + * @returns A tuple of: + * - tmin: minimum azimuth time (s) + * - tmax: maximum azimuth time (s) + * - rmin: minimum slant range (m) + * - rmax: maximum slant range (m) + * - error code (non-zero if any perimeter point fails to + * converge in polar2geo or geo2rdr) + */ +std::tuple +findPolarGridBoundingBoxInRadarCoord( + const PolarGrid& polar_grid, + const isce3::core::Orbit& orbit, + const isce3::core::LUT2d& doppler, + const double wavelength, + const isce3::core::LookSide lookside, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params, + const int nextra = 0); + +/** + * @brief Find the subset of a radar grid covered by a polar grid. + * + * Computes the bounding box of the polar grid in stripmap radar + * coordinates using findPolarGridBoundingBoxInRadarCoord, then + * converts this bounding box to integer radar grid indices + * (azimuth line, range sample). If the polar grid does not + * overlap the radar grid at all, a zero-sized subset (0, 0, 0, 0) + * is returned. + * + * @param[in] polar_grid Input polar grid + * @param[in] radar_geom Target radar geometry grid + * @param[in] dem Digital elevation model (DEM) + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * @param[in] g2r_params geo2rdr_bracket configuration parameters + * @param[in] nextra Number of extra perimeter points per + * edge (defaults to 0) + * + * @returns A tuple of: + * - i0: starting azimuth line index + * - i1: ending azimuth line index (exclusive) + * - j0: starting range sample index + * - j1: ending range sample index (exclusive) + * - error code (from findPolarGridBoundingBoxInRadarCoord) + */ +std::tuple +findPolarGridBoundingBoxInRadarGrid( + const PolarGrid& polar_grid, + const isce3::container::RadarGeometry& radar_geom, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params, + const isce3::geometry::detail::Geo2RdrBracketParams& g2r_params, + const int nextra = 0); + +/** + * @brief Compute 3D geo coordinates for a radar grid. + * + * Computes the 3D XYZ position for every pixel in a radar geometry + * grid using rdr2geo_bracket. Returns the points in a vector along with an + * error code indicating success or failure. + * + * @param[in] geom Radar geometry grid, orbit, & doppler + * @param[in] dem Digital elevation model (DEM) + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * + * @returns A tuple containing: + * - vector of 3D XYZ positions (ECEF, m), one per pixel + * - error code (non-zero if rdr2geo fails to converge + * for any pixel) + */ +std::tuple, isce3::error::ErrorCode> +computeRadarGridGeoPoints( + const isce3::container::RadarGeometry& geom, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params); + +/** + * @brief Compute 3D geo coordinates for a radar grid (pre-allocated). + * + * Computes the 3D XYZ position for every pixel in a radar geometry + * grid using rdr2geo_bracket, writing results into a pre-allocated buffer. + * The caller must allocate the buffer to at least + * geom.gridLength() * geom.gridWidth() elements. + * + * @param[out] points Output 3D XYZ positions (ECEF, m) + * @param[in] geom Radar geometry grid, orbit, and doppler + * @param[in] dem Digital elevation model (DEM) + * @param[in] r2g_params rdr2geo_bracket configuration parameters + * + * @returns Error code indicating success or failure of the + * computation (non-zero if rdr2geo fails to converge + * for any pixel) + */ +isce3::error::ErrorCode +computeRadarGridGeoPoints( + isce3::core::Vec3* points, + const isce3::container::RadarGeometry& geom, + const isce3::geometry::DEMInterpolator& dem, + const isce3::geometry::detail::Rdr2GeoBracketParams& r2g_params); + } // namespace focus } // namespace isce3 diff --git a/cxx/isce3/geometry/detail/Rdr2Geo.h b/cxx/isce3/geometry/detail/Rdr2Geo.h index d14b2ba72..61923edcc 100644 --- a/cxx/isce3/geometry/detail/Rdr2Geo.h +++ b/cxx/isce3/geometry/detail/Rdr2Geo.h @@ -130,6 +130,62 @@ rdr2geo_bracket(isce3::core::Vec3* xyz, double wavelength, isce3::core::LookSide side, const Rdr2GeoBracketParams& params = {}); + +/** + * \internal + * Lower level version of rdr2geo_bracket that works directly in polar + * coordinates. Also avoids repeated Orbit interpolations and trig calls. + * + * @param[out] xyz Output target ECEF XYZ position (m) + * @param[out] lookAngle Output pseudo-look angle of target (rad) + * @param[in] origin Origin of the polar grid, ECEF XYZ (m) + * @param[in] axis Along-track axis of polar grid, unit ECEF XYZ + * @param[in] slantRange Distance from origin to target (m) + * @param[in] sinSquint Sine of squint angle + * @param[in] cosSquint Cosine of squint angle + * @param[in] dem Digital elevation model (m above ellipsoid) + * @param[in] ellipsoid Ellipsoid associated with DEM + * @param[in] side Look direction (Left or Right) + * @param[in] params Root finding algorithm parameters + */ +template +CUDA_HOSTDEV isce3::error::ErrorCode +polar2geo_bracket(isce3::core::Vec3* xyz, double* lookAngle, + const isce3::core::Vec3& origin, const isce3::core::Vec3& axis, + const double slantRange, const double sinSquint, const double cosSquint, + const DEMInterpolator& dem, const isce3::core::Ellipsoid& ellipsoid, + isce3::core::LookSide side, const Rdr2GeoBracketParams& params); + + +/** + * \internal + * Low level version of polar2polar_bracket. + * + * @param[in] outSinSquint Sine of squint angle in output grid. + * @param[in] outRange Distance from output origin to target (m) + * @param[in] inSinSquint Sine of squint angle in input grid. + * @param[in] inCosSquint Cosine of squint angle in input grid + * Equal to sqrt(1 - inSinSquint^2) + * @param[in] inRange Distance from input origin to target (m) + * @param[in] inOrigin Origin of the input polar grid, ECEF XYZ (m) + * @param[in] inAxis Along-track axis of input polar grid, unit ECEF XYZ + * @param[in] outOrigin Origin of the output polar grid, ECEF XYZ (m) + * @param[in] outAxis Along-track axis of output polar grid, unit ECEF XYZ + * @param[in] dem Digital elevation model (m above ellipsoid) + * @param[in] ellipsoid Ellipsoid associated with DEM + * @param[in] side Look direction (Left or Right) + * @param[in] params Root finding algorithm parameters + */ +template +CUDA_HOSTDEV isce3::error::ErrorCode +polar2polar_bracket(double* outSinSquint, double* outRange, + const double inSinSquint, const double inCosSquint, + const double inRange, + const isce3::core::Vec3& inOrigin, const isce3::core::Vec3& inAxis, + const isce3::core::Vec3& outOrigin, const isce3::core::Vec3& outAxis, + const DEMInterpolator& dem, const isce3::core::Ellipsoid& ellipsoid, + const isce3::core::LookSide& side, const Rdr2GeoBracketParams& params); + }}} // namespace isce3::geometry::detail #include "Rdr2Geo.icc" diff --git a/cxx/isce3/geometry/detail/Rdr2Geo.icc b/cxx/isce3/geometry/detail/Rdr2Geo.icc index f5fdac880..93fb4e8c4 100644 --- a/cxx/isce3/geometry/detail/Rdr2Geo.icc +++ b/cxx/isce3/geometry/detail/Rdr2Geo.icc @@ -170,41 +170,27 @@ rdr2geo(isce3::core::Vec3* llh, const isce3::core::Pixel& pixel, NVCC_HD_WARNING_DISABLE -template +template CUDA_HOSTDEV isce3::error::ErrorCode -rdr2geo_bracket(isce3::core::Vec3* xyz, - double aztime, double slantRange, double doppler, const Orbit& orbit, +polar2geo_bracket(isce3::core::Vec3* xyz, double* lookAngle, + const isce3::core::Vec3& origin, const isce3::core::Vec3& axis, + const double slantRange, const double sinSquint, const double cosSquint, const DEMInterpolator& dem, const isce3::core::Ellipsoid& ellipsoid, - double wavelength, isce3::core::LookSide side, - const Rdr2GeoBracketParams& params) + isce3::core::LookSide side, const Rdr2GeoBracketParams& params) { using namespace isce3::core; using isce3::error::ErrorCode; - // Interpolate orbit to get radar position and velocity. - Vec3 radarXYZ, velocity; - const auto errCode = orbit.interpolate(&radarXYZ, &velocity, aztime); - if (errCode != ErrorCode::Success) { - return errCode; - } - const double speed = velocity.norm(); - - // Construct some useful basis vectors. + // Construct some useful basis vectors. Naming assumes axis is along-track. // For simplicity we'll use geocentric nadir (not geodetic), // which determines our definition of look angle. - const Vec3 alongTrack = velocity / speed; - const Vec3 right = alongTrack.cross(radarXYZ).normalized(); - const Vec3 down = alongTrack.cross(right); + const Vec3 right = axis.cross(origin).normalized(); + const Vec3 down = axis.cross(right); const Vec3 horizontal = (side == LookSide::Right) ? right : -right; - // Convert Doppler into equivalent squint angle. - const double sinSquint = doppler * wavelength / (2 * speed); - // NOTE squint is in [-90, 90] by definition. - const double cosSquint = std::sqrt(1.0 - sinSquint * sinSquint); - - // Doppler cone and range sphere intersect to make a circle with - // the following parameters. - const Vec3 center = radarXYZ + sinSquint * slantRange * alongTrack; + // Axis & angle define a cone, and range defines a sphere. + // These intersect to make a circle with the following parameters. + const Vec3 center = origin + sinSquint * slantRange * axis; const double radius = cosSquint * slantRange; // Parameterize points on this circle with an angle similar to look angle. @@ -221,9 +207,9 @@ rdr2geo_bracket(isce3::core::Vec3* xyz, // Now we just need to find the angle on this circle where it intersects // the DEM. Create an error function we'll feed to a root finder. auto dh = [&](double look) { - const Vec3 xyz = getXYZ(look); + const Vec3 xyz_ = getXYZ(look); Vec3 lonLatH; - ellipsoid.xyzToLonLat(xyz, lonLatH); + ellipsoid.xyzToLonLat(xyz_, lonLatH); return lonLatH[2] - dem.interpolateLonLat(lonLatH[0], lonLatH[1]); }; @@ -231,14 +217,72 @@ rdr2geo_bracket(isce3::core::Vec3* xyz, const double tolLook = params.tol_height / radius; // Solve. - double lookSolution = 0.0; + *lookAngle = 0.0; const auto err = isce3::math::find_zero_brent( - params.look_min, params.look_max, dh, tolLook, &lookSolution); + params.look_min, params.look_max, dh, tolLook, lookAngle); if (err != ErrorCode::Success) { return err; } - *xyz = getXYZ(lookSolution); + *xyz = getXYZ(*lookAngle); + return ErrorCode::Success; + +} + + +NVCC_HD_WARNING_DISABLE +template +CUDA_HOSTDEV isce3::error::ErrorCode +polar2polar_bracket(double* outSinSquint, double* outRange, + const double inSinSquint, const double inCosSquint, + const double inRange, + const isce3::core::Vec3& inOrigin, const isce3::core::Vec3& inAxis, + const isce3::core::Vec3& outOrigin, const isce3::core::Vec3& outAxis, + const DEMInterpolator& dem, const isce3::core::Ellipsoid& ellipsoid, + const isce3::core::LookSide& side, const Rdr2GeoBracketParams& params) +{ + using isce3::error::ErrorCode; + isce3::core::Vec3 xyz; + double lookAngle; + auto errCode = polar2geo_bracket(&xyz, &lookAngle, inOrigin, inAxis, + inRange, inSinSquint, inCosSquint, dem, ellipsoid, side, params); + if (errCode != ErrorCode::Success) { + return errCode; + } + // geo2polar + const isce3::core::Vec3 lookvec = xyz - outOrigin; + *outRange = lookvec.norm(); + *outSinSquint = lookvec.dot(outAxis) / (*outRange); return ErrorCode::Success; } + +NVCC_HD_WARNING_DISABLE +template +CUDA_HOSTDEV isce3::error::ErrorCode +rdr2geo_bracket(isce3::core::Vec3* xyz, + double aztime, double slantRange, double doppler, const Orbit& orbit, + const DEMInterpolator& dem, const isce3::core::Ellipsoid& ellipsoid, + double wavelength, isce3::core::LookSide side, + const Rdr2GeoBracketParams& params) +{ + using isce3::error::ErrorCode; + + // Interpolate orbit to get radar position and velocity. + isce3::core::Vec3 radarXYZ, velocity; + const auto errCode = orbit.interpolate(&radarXYZ, &velocity, aztime); + if (errCode != ErrorCode::Success) { + return errCode; + } + const auto speed = velocity.norm(); + + // Convert Doppler into equivalent squint angle. + const double sinSquint = doppler * wavelength / (2 * speed); + // NOTE squint is in [-90, 90] by definition. + const double cosSquint = std::sqrt(1.0 - sinSquint * sinSquint); + + double lookAngle = 0; + return polar2geo_bracket(xyz, &lookAngle, radarXYZ, velocity / speed, + slantRange, sinSquint, cosSquint, dem, ellipsoid, side, params); +} + }}} // namespace isce3::geometry::detail diff --git a/cxx/isce3/geometry/geometry.cpp b/cxx/isce3/geometry/geometry.cpp index a48881b5f..dbe307f66 100644 --- a/cxx/isce3/geometry/geometry.cpp +++ b/cxx/isce3/geometry/geometry.cpp @@ -642,4 +642,16 @@ double isce3::geometry::compute_mean_dem(const DEMInterpolator& dem) } return dem.refHeight(); } -// end of file + + +CUDA_HOSTDEV +isce3::error::ErrorCode +isce3::geometry::geo2polar(double* sinSquint, double* range, + const Vec3& xyz, const Vec3& origin, + const Vec3& axis) +{ + const Vec3 lookvec = xyz - origin; + *range = lookvec.norm(); + *sinSquint = lookvec.dot(axis) / (*range); + return isce3::error::ErrorCode::Success; +} diff --git a/cxx/isce3/geometry/geometry.h b/cxx/isce3/geometry/geometry.h index 1171443aa..0c7c8f897 100644 --- a/cxx/isce3/geometry/geometry.h +++ b/cxx/isce3/geometry/geometry.h @@ -16,6 +16,7 @@ #include "forward.h" #include +#include #include #include @@ -405,5 +406,20 @@ std::tuple lookIncAngFromSlantRange( */ double compute_mean_dem(const DEMInterpolator& dem); +/** Convert 3D position to polar coordinates. + * + * @param[out] sinSquint Sine of the squint angle (complement of the angle + * between the azimuth axis and the radar-to-target + * line of sight). + * @param[out] range Distance to the target (m) + * @param[in] xyz Target ECEF XYZ position (m) + * @param[in] origin Origin of polar coordinate system, ECEF XYZ (m) + * @param[in] axis Azimuth axis, unit ECEF XYZ + */ +CUDA_HOSTDEV +isce3::error::ErrorCode +geo2polar(double* sinSquint, double* range, const isce3::core::Vec3& xyz, + const isce3::core::Vec3& origin, const isce3::core::Vec3& axis); + } // namespace geometry } // namespace isce3 diff --git a/cxx/isce3/signal/NFFT.cpp b/cxx/isce3/signal/NFFT.cpp index 978183754..4e65d8c70 100644 --- a/cxx/isce3/signal/NFFT.cpp +++ b/cxx/isce3/signal/NFFT.cpp @@ -98,12 +98,11 @@ set_spectrum(const std::valarray> &x) template std::complex isce3::signal::NFFT:: -interp(double t) const +interp(double t, bool periodic) const { // scale time index to account for zero-padding of spectrum. t *= (double)_fft_size / (double)_n; - return isce3::core::interp1d>(_kernel, _xt, t, - /*periodic*/true); + return isce3::core::interp1d>(_kernel, _xt, t, periodic); } template diff --git a/cxx/isce3/signal/NFFT.h b/cxx/isce3/signal/NFFT.h index 63c4d3046..910db3499 100644 --- a/cxx/isce3/signal/NFFT.h +++ b/cxx/isce3/signal/NFFT.h @@ -153,12 +153,14 @@ class isce3::signal::NFFT { /** Interpolate the transformed signal. * - * @param[in] t Location in [0,n) to sample the time-domain signal. + * @param[in] t Location in [0,n) to sample the time-domain + * signal. + * @param[in] periodic Whether to treat the domain as periodic. * * @see execute is an alternative strategy. * @see set_spectrum must be called first. */ - std::complex interp(double t) const; + std::complex interp(double t, bool periodic = true) const; size_t size_kernel() const {return 2*_m+1;} size_t size_spectrum() const {return _n;} diff --git a/cxx/isce3/signal/NFFT2d.cpp b/cxx/isce3/signal/NFFT2d.cpp new file mode 100644 index 000000000..9fc16e4ce --- /dev/null +++ b/cxx/isce3/signal/NFFT2d.cpp @@ -0,0 +1,222 @@ +#include "NFFT2d.h" + +#include +#include + +using isce3::core::NFFTKernel; + +namespace isce3::signal { + +// constructor +template +NFFT2d::NFFT2d( + const dims_t& m, const dims_t& sizes, const dims_t& fft_sizes) + : m_(m), sizes_(sizes), fft_sizes_(fft_sizes), kernels_( + {NFFTKernel(m[0], sizes[0], fft_sizes[0]), + NFFTKernel(m[1], sizes[1], fft_sizes[1])}) + +{ + size_t nout = static_cast(fft_sizes[0]) * fft_sizes[1]; + xf_.resize(nout); + + // Pre-compute spectral weights (1/phi_hat in NFFT papers). + // Also include factor of n since FFTW does not normalize DFT. + for (int idim = 0; idim < ndims; ++idim) { + weights_[idim].resize(sizes[idim]); + T b = M_PI * (2.0 - 1.0 * sizes[idim] / fft_sizes[idim]); + T norm = isce3::math::bessel_i0(b * m[idim]) / sizes[idim]; + size_t n2 = (sizes[idim] - 1) / 2 + 1; + for (size_t i = 0; i < n2; ++i) { + double f = 2 * M_PI * i / fft_sizes_[idim]; + weights_[idim][i] = norm / + isce3::math::bessel_i0(m[idim] * std::sqrt(b * b - f * f)); + } + for (size_t i = n2; i < sizes[idim]; ++i) { + double f = 2 * M_PI * ((double)i - sizes[idim]) / fft_sizes[idim]; + weights_[idim][i] = norm / + isce3::math::bessel_i0(m[idim] * std::sqrt(b * b - f * f)); + } + } +} + +// Digest some data. +template +NFFT2dResult +NFFT2d::transform(const dims_t& sizes, + const dims_t& strides, const std::complex *x) +{ + for (int idim = 0; idim < ndims; ++idim) { + if (sizes[idim] != sizes_[idim]) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "Spectrum size != NFFT size."); + } + } + // Clear any old data. + size_t nout = static_cast(fft_sizes_[0]) * fft_sizes_[1]; + xf_.assign(nout, std::complex(0, 0)); + + const size_t m2 = sizes_[0] / 2; + const size_t n2 = sizes_[1] / 2; + + // Zero-pad and scale spectrum. + // For even n the non-zero intervals are [0, n // 2) and [-n // 2, n - 1] + // (Assuming we don't bother splitting the Nyquist bin). For odd n it's + // the same except it's symmetric [0, n // 2] and [-n // 2, n - 1] with + // both intervals *closed*. + const auto row_end = (sizes_[0] % 2 == 0) ? m2 : m2 + 1; + const auto col_end = (sizes_[1] % 2 == 0) ? n2 : n2 + 1; + #pragma omp parallel for + for (size_t i = 0; i < row_end; ++i) { + const auto wi = weights_[0][i]; + // pointer to row i of input data + const auto pxi = x + (strides[0] * i); + // pointer to row i of ifft buffer + const auto pxfi = xf_.data() + (fft_sizes_[1] * i); + // columns [0, n2) or [0, n2] + for (size_t j = 0; j < col_end; ++j) { + const auto wj = weights_[1][j]; + pxfi[j] = wi * wj * pxi[strides[1] * j]; + } + // columns [-n2, 0) + for (size_t j = n2; j > 0; --j) { + const auto wj = weights_[1][sizes_[1] - j]; + pxfi[fft_sizes_[1] - j] = wi * wj * pxi[strides[1] * (sizes_[1] - j)]; + } + } + #pragma omp parallel for + for (size_t i = m2; i > 0; --i) { + const auto wi = weights_[0][sizes_[0] - i]; + // pointer to row (ny - i) + const auto pxi = x + (strides[0] * (sizes_[0] - i)); + const auto pxfi = xf_.data() + (fft_sizes_[1] * (fft_sizes_[0] - i)); + for (size_t j = 0; j < col_end; ++j) { + const auto wj = weights_[1][j]; + pxfi[j] = wi * wj * pxi[strides[1] * j]; + } + for (size_t j = n2; j > 0; --j) { + const auto wj = weights_[1][sizes_[1] - j]; + pxfi[fft_sizes_[1] - j] = wi * wj * pxi[strides[1] * (sizes_[1] - j)]; + } + } + + // Allocate result object, with friend access to storage. + auto result = NFFT2dResult(m_, sizes_, fft_sizes_, kernels_); + + // NOTE For even lengths we're not splitting Nyquist bin. + // Transform to (expanded) time-domain. + const int dims[] = {fft_sizes_[0], fft_sizes_[1]}; + isce3::fft::ifft2d(result.data(), xf_.data(), dims); + + return result; +} + +template +std::complex NFFT2dResult::interp(const std::array& t, bool periodic) const +{ + constexpr int xdim = 1, ydim = 0; + + // scale time index to account for zero-padding of spectrum. + double x = t[xdim] * fft_sizes_[xdim] / sizes_[xdim]; + double y = t[ydim] * fft_sizes_[ydim] / sizes_[ydim]; + + return isce3::core::interp2d>(kernels_[xdim], + kernels_[ydim], xt_.data(), fft_sizes_[xdim], /* stridex */ 1, + fft_sizes_[ydim], /* stridey */ fft_sizes_[xdim], x, y, periodic); +} + + +template +NFFT2dResult makeImageNFFT2d( + const Eigen::Ref>>& image, + const NFFT2dParams& params, + bool pad_input) +{ + using isce3::fft::nextFastPower; + + // NFFTKernel width = 2*m+1. The CPU interp uses heap allocation so large m + // is okay, but the CUDA interp2d uses fixed-size stack arrays + // (MAX_WIDTH=16), so reject m > 7 here to avoid confusion when the same + // params work on CPU but overflow on GPU. + constexpr int MAX_NFFT_M = 7; + if (params.rows.m > MAX_NFFT_M or params.cols.m > MAX_NFFT_M) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "NFFT kernel m must be <= 7 (kernel width = 2*m+1 must fit in " + "CUDA device stack arrays of size 16)"); + } + + auto rows_in = image.rows(); + auto cols_in = image.cols(); + using image_t = isce3::core::EArray2D>; + auto image_copy = image_t(0, 0); + + // Pointer to input image or padded/copied version so we can have fewer + // conditionals later. + // FIXME figure out how to do this with an Eigen type... + auto image_ptr = image.data(); + + // Need to copy if image is not contiguous row-major since we don't have + // high-level interface for strided FFTs. + bool need_copy = (image.innerStride() != 1) or (image.outerStride() != cols_in); + if (need_copy) { + image_copy.resize(rows_in, cols_in); + // assign later + } + + if (pad_input) { + auto padded_rows_in = nextFastPower(rows_in); + auto padded_cols_in = nextFastPower(cols_in); + if ((rows_in == padded_rows_in) and (cols_in == padded_cols_in)) { + // User asked for padding but we don't actually need it. + pad_input = false; + } else { + image_copy.resize(padded_rows_in, padded_cols_in); + image_copy.setZero(); + rows_in = padded_rows_in; + cols_in = padded_cols_in; + // assign later + } + } + + if (need_copy or pad_input) { + // This way NFFT2d::interp() coordinates are preserved, though user + // will be able to get some extra data. + image_copy.topLeftCorner(image.rows(), image.cols()) = image; + image_ptr = image_copy.data(); + } + + // Use fft2 b/c planfft2d could modify inputs and we won't reuse it anyway. + using dims_t = typename NFFT2d::dims_t; + dims_t dims = { + static_cast(rows_in), + static_cast(cols_in)}; + auto spectrum = image_t(dims[0], dims[1]); + isce3::fft::fft2d(spectrum.data(), image_ptr, {dims[0], dims[1]}); + + // Calculate sizes for padded inverse transform. + dims_t dims_out = { + nextFastPower(static_cast(std::round(params.rows.s * dims[0]))), + nextFastPower(static_cast(std::round(params.cols.s * dims[1])))}; + + const dims_t m = {params.rows.m, params.cols.m}; + auto plan = NFFT2d(m, dims, dims_out); + return plan.transform(dims, {dims[1], 1}, spectrum.data()); +} + +} + +template class isce3::signal::NFFT2d; +template class isce3::signal::NFFT2d; +template class isce3::signal::NFFT2dResult; +template class isce3::signal::NFFT2dResult; + +template isce3::signal::NFFT2dResult +isce3::signal::makeImageNFFT2d( + const Eigen::Ref>>& image, + const isce3::signal::NFFT2dParams& params, + bool pad_input); + +template isce3::signal::NFFT2dResult +isce3::signal::makeImageNFFT2d( + const Eigen::Ref>>& image, + const isce3::signal::NFFT2dParams& params, + bool pad_input); diff --git a/cxx/isce3/signal/NFFT2d.h b/cxx/isce3/signal/NFFT2d.h new file mode 100644 index 000000000..4701419c1 --- /dev/null +++ b/cxx/isce3/signal/NFFT2d.h @@ -0,0 +1,156 @@ +#pragma once + +#include "forward.h" + +#include +#include +#include + +#include +#include + + +namespace isce3::signal { + +template +class NFFT2d { + public: + static constexpr int ndims = 2; + using dims_t = std::array; + + NFFT2d() = delete; + + /** + * @brief Construct a new NFFT2d object + * + * @param m Interpolator half-length along {rows, columns} + * @param sizes Image spectrum dimensions {rows, columns} + * @param fft_sizes Transform sizes along {rows, columns}. + * Usually larger than image size. + */ + NFFT2d(const dims_t& m, const dims_t& sizes, + const dims_t& fft_sizes); + + /** + * @brief Ingest the image spectrum. + * + * @param sizes Image spectrum dimensions {rows, columns}. + * Must match dimensions provided in ctor. + * @param strides Strides (in pixels) along each dimension {rows, columns}. + * @param x Image spectrum. + * + * The spectrum will be zero-padded, pre-filtered, and transformed + * to the time-domain. + */ + NFFT2dResult transform(const dims_t& sizes, const dims_t& strides, + const std::complex *x); + + /** Image spectrum dimensions */ + const dims_t& sizes() const { return sizes_; } + + /** Transform sizes */ + const dims_t& fft_sizes() const { return fft_sizes_; } + + /** Pointer to most recent spectral data (filtered and padded) */ + const std::complex* spectrum() const { return xf_.data(); } + + private: + dims_t m_, sizes_, fft_sizes_; + std::vector> xf_; + std::array, 2> weights_; + std::array, 2> kernels_; +}; + + +template +class NFFT2dResult { + + public: + using dims_t = typename NFFT2d::dims_t; + + NFFT2dResult() = delete; + + NFFT2dResult( + const dims_t& m, + const dims_t& sizes, + const dims_t& fft_sizes, + const std::array, 2>& kernels, + const std::complex* xt = nullptr) : + m_{m}, sizes_{sizes}, fft_sizes_{fft_sizes}, kernels_{kernels} + { + const auto n = static_cast(fft_sizes[0]) * fft_sizes[1]; + if (xt == nullptr) { + xt_.resize(n); + } else { + xt_.assign(xt, xt + n); + } + }; + + /** + * @brief Interpolate the image + * + * @param t Desired pixel location {row, column} + * @param periodic Whether to use a periodic boundary condition. + * @return Interpolated value. + */ + std::complex interp(const std::array& t, + bool periodic = true) const; + + const std::complex* data() const { return xt_.data(); } + std::complex* data() { return xt_.data(); } + + const dims_t& kernel_radii() const { return m_; } + const dims_t& sizes() const { return sizes_; } + const dims_t& fft_sizes() const { return fft_sizes_; } + const auto& kernels() const { return kernels_; } + + private: + dims_t m_, sizes_, fft_sizes_; + std::array, 2> kernels_; + std::vector> xt_; +}; + + +struct NFFTParams { + int m; /// half width of interpolator + double s; /// oversampling factor + + NFFTParams(int m_ = 2, double s_ = 2.0) : m{m_}, s{s_} { + if (m_ < 1) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Need interpolator size m >= 1 for NFFT"); + } + if (s_ <= 1.0) { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "Need oversampling ratio s > 1.0 for NFFT"); + } + } +}; + +struct NFFT2dParams { + NFFTParams rows, cols; +}; + + +/** + * @brief Create an NFFT2d object for interpolating an image. + * + * @tparam T Format of real/imag pixel data, typically float or double + * @param image Input time-domain image. A temporary copy will be made if + * it is not row-major with a column stride of one. + * @param m Half-length of interpolator along {rows, columns} + * @param s Minimum factors (> 1) for frequency-domain zero-padding + * along {rows, columns}. Actual padding may be larger to + * achieve efficient inverse transform size. + * @param pad_input Whether to also zero-pad input data to an efficient + * forward transform size. Requires extra memory. + * + * @return NFFT2d object for interpolating the image. + */ +template +NFFT2dResult makeImageNFFT2d( + const Eigen::Ref>>& image, + const NFFT2dParams& params = {}, + bool pad_input = false); + +} // namespace isce3::signal \ No newline at end of file diff --git a/cxx/isce3/signal/forward.h b/cxx/isce3/signal/forward.h index 883e54292..359121057 100644 --- a/cxx/isce3/signal/forward.h +++ b/cxx/isce3/signal/forward.h @@ -3,10 +3,14 @@ namespace isce3 { namespace signal { class Crossmul; + class NFFTParams; + class NFFT2dParams; template class Covariance; template class Filter; template class Looks; template class NFFT; + template class NFFT2d; + template class NFFT2dResult; template class Signal; template class FilterData; }} diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 3e253c187..f2a68f277 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -210,3 +210,15 @@ @inbook{brent year={1973}, publisher={Prentice-Hall}, } + +@inproceedings{yegulalp1999, + author={Yegulalp, A.F.}, + booktitle={Proceedings of the 1999 IEEE Radar Conference. Radar into the Next Millennium (Cat. No.99CH36249)}, + title={Fast backprojection algorithm for synthetic aperture radar}, + year={1999}, + volume={}, + number={}, + pages={60-65}, + doi={10.1109/NRC.1999.767270}, +} + diff --git a/extern/CMakeLists.txt b/extern/CMakeLists.txt index 3cafc9928..517035d85 100644 --- a/extern/CMakeLists.txt +++ b/extern/CMakeLists.txt @@ -95,12 +95,17 @@ endmacro() macro(getpackage_openmp_optional) # Check for OpenMP (optional dependency). # If not found, default to an empty placeholder target. - find_package(OpenMP OPTIONAL_COMPONENTS CXX) + find_package(OpenMP OPTIONAL_COMPONENTS CXX CUDA) add_library(OpenMP::OpenMP_CXX_Optional INTERFACE IMPORTED) + add_library(OpenMP::OpenMP_CUDA_Optional INTERFACE IMPORTED) if(TARGET OpenMP::OpenMP_CXX) target_link_libraries(OpenMP::OpenMP_CXX_Optional INTERFACE OpenMP::OpenMP_CXX) endif() + if(TARGET OpenMP::OpenMP_CUDA) + target_link_libraries(OpenMP::OpenMP_CUDA_Optional + INTERFACE OpenMP::OpenMP_CUDA) + endif() endmacro() macro(getpackage_pybind11) diff --git a/python/extensions/pybind_isce3/Sources.cmake b/python/extensions/pybind_isce3/Sources.cmake index 6732d3dcc..10bf21220 100644 --- a/python/extensions/pybind_isce3/Sources.cmake +++ b/python/extensions/pybind_isce3/Sources.cmake @@ -54,6 +54,7 @@ geometry/RTC.cpp geometry/metadataCubes.cpp geometry/ltpcoordinates.cpp geometry/pntintersect.cpp +geogrid/geogrid_ecef_coords.cpp geogrid/getRadarGrid.cpp geogrid/relocateRaster.cpp geogrid/geogrid.cpp @@ -83,6 +84,7 @@ signal/CrossMultiply.cpp signal/flatten.cpp signal/filter2D.cpp signal/multilook.cpp +signal/NFFT2d.cpp product/GeoGridParameters.cpp product/product.cpp product/RadarGridParameters.cpp @@ -107,7 +109,7 @@ if(WITH_CUDA) cuda/geometry/geometry.cpp cuda/geometry/geo2rdr.cpp cuda/geometry/rdr2geo.cpp - cuda/focus/Backproject.cpp + cuda/focus/Backproject.cu cuda/focus/focus.cpp cuda/image/image.cpp cuda/image/Resample.cpp @@ -116,5 +118,6 @@ if(WITH_CUDA) cuda/matchtemplate/pycuampcor.cpp cuda/signal/signal.cpp cuda/signal/Crossmul.cpp + cuda/signal/NFFT2d.cu ) endif() diff --git a/python/extensions/pybind_isce3/core/Interp1d.cpp b/python/extensions/pybind_isce3/core/Interp1d.cpp index 3ec026d38..005867f52 100644 --- a/python/extensions/pybind_isce3/core/Interp1d.cpp +++ b/python/extensions/pybind_isce3/core/Interp1d.cpp @@ -15,7 +15,8 @@ using isce3::except::RuntimeError; template static py::object -interp_duckt(const Kernel & kernel, py::buffer_info & info, py::object t) +interp_duckt(const Kernel & kernel, py::buffer_info & info, + py::object t, bool periodic) { DataType* data = static_cast(info.ptr); int stride = info.strides[0] / sizeof(DataType); @@ -31,12 +32,12 @@ interp_duckt(const Kernel & kernel, py::buffer_info & info, py::obje py::gil_scoped_release release; #pragma omp parallel for for (size_t i=0; i < ta.size(); ++i) { - outbuf[i] = interp1d(kernel, data, n, stride, ta(i)); + outbuf[i] = interp1d(kernel, data, n, stride, ta(i), periodic); } } else { // can't release GIL since kernel is a Python object for (size_t i=0; i < ta.size(); ++i) { - outbuf[i] = interp1d(kernel, data, n, stride, ta(i)); + outbuf[i] = interp1d(kernel, data, n, stride, ta(i), periodic); } } return out; @@ -47,7 +48,7 @@ interp_duckt(const Kernel & kernel, py::buffer_info & info, py::obje template static py::object -interp_duckbuf(Kernel & kernel, py::buffer buf, py::object t) +interp_duckbuf(Kernel & kernel, py::buffer buf, py::object t, bool periodic) { py::buffer_info info = buf.request(); using C8 = std::complex; @@ -56,30 +57,30 @@ interp_duckbuf(Kernel & kernel, py::buffer buf, py::object t) throw RuntimeError(ISCE_SRCINFO(), "data buffer must be 1-D"); } if (info.format == py::format_descriptor::format()) { - return interp_duckt(kernel, info, t); + return interp_duckt(kernel, info, t, periodic); } else if (info.format == py::format_descriptor::format()) { - return interp_duckt(kernel, info, t); + return interp_duckt(kernel, info, t, periodic); } else if (info.format == py::format_descriptor::format()) { - return interp_duckt(kernel, info, t); + return interp_duckt(kernel, info, t, periodic); } else if (info.format == py::format_descriptor::format()) { - return interp_duckt(kernel, info, t); + return interp_duckt(kernel, info, t, periodic); } throw RuntimeError(ISCE_SRCINFO(), "Unsupported types for interp1d"); } void addbinding_interp1d(py::module & m) { - m.def("interp1d", [](py::object pyKernel, py::buffer buf, py::object t) { + m.def("interp1d", [](py::object pyKernel, py::buffer buf, py::object t, bool periodic) { if (py::isinstance>(pyKernel)) { auto kernel = pyKernel.cast *>(); - return interp_duckbuf(*kernel, buf, t); + return interp_duckbuf(*kernel, buf, t, periodic); } else if (py::isinstance>(pyKernel)) { auto kernel = pyKernel.cast *>(); - return interp_duckbuf(*kernel, buf, t); + return interp_duckbuf(*kernel, buf, t, periodic); } throw RuntimeError(ISCE_SRCINFO(), "Expected Kernel or KernelF32"); }, @@ -88,5 +89,5 @@ void addbinding_interp1d(py::module & m) units are sample numbers (starting at zero), and `time` may be a scalar or an array. )", - py::arg("kernel"), py::arg("data"), py::arg("time")); + py::arg("kernel"), py::arg("data"), py::arg("time"), py::arg("periodic") = false); } diff --git a/python/extensions/pybind_isce3/core/Linspace.cpp b/python/extensions/pybind_isce3/core/Linspace.cpp index 8f168b2f6..d5353546e 100644 --- a/python/extensions/pybind_isce3/core/Linspace.cpp +++ b/python/extensions/pybind_isce3/core/Linspace.cpp @@ -108,6 +108,10 @@ void addbinding(py::class_>& pyLinspace) return std::make_tuple(self.size()); }) + .def_property_readonly("bounds", [](const Linspace& self) { + + return self.bounds(); + }) // methods .def("resize", [](Linspace& self, int size) { @@ -121,6 +125,12 @@ void addbinding(py::class_>& pyLinspace) .def("search", [](const Linspace& self, T x) { return self.search(x); }, "Return the position where the specified value would be inserted " "in the sequence in order to maintain sorted order.") + + .def("__repr__", [](const Linspace& self) { + return std::string("Linspace(" + std::to_string(self.first()) + + ", " + std::to_string(self.spacing()) + ", " + + std::to_string(self.size()) + ")"); + }) ; } diff --git a/python/extensions/pybind_isce3/cuda/focus/Backproject.cpp b/python/extensions/pybind_isce3/cuda/focus/Backproject.cpp deleted file mode 100644 index 0655cf6e2..000000000 --- a/python/extensions/pybind_isce3/cuda/focus/Backproject.cpp +++ /dev/null @@ -1,117 +0,0 @@ -#include "Backproject.h" -#include "pybind_isce3/focus/Backproject.h" // parse parameter dicts - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace py = pybind11; - -using namespace isce3::cuda::focus; -using namespace isce3::except; - -using isce3::container::RadarGeometry; -using isce3::core::Kernel; -using isce3::error::ErrorCode; -using isce3::focus::parseDryTropoModel; -using isce3::geometry::DEMInterpolator; - -void addbinding_cuda_backproject(py::module& m) -{ - m.def("backproject", []( - py::array_t, py::array::c_style> out, - const RadarGeometry& out_geometry, - py::array_t, py::array::c_style> in, - const RadarGeometry& in_geometry, - const DEMInterpolator& dem, - double fc, - double ds, - const Kernel& kernel, - const std::string& dry_tropo_model, - py::dict rdr2geo_params, - py::dict geo2rdr_params, - int batch, - std::optional> height) { - - if (out.ndim() != 2) { - throw InvalidArgument(ISCE_SRCINFO(), "output array must be 2-D"); - } - - if (out.shape()[0] != out_geometry.gridLength() or - out.shape()[1] != out_geometry.gridWidth()) { - - std::string errmsg = "output array shape must match output " - "radar grid shape"; - throw InvalidArgument(ISCE_SRCINFO(), errmsg); - } - - if (in.ndim() != 2) { - throw InvalidArgument(ISCE_SRCINFO(), "input signal data must be 2-D"); - } - - if (in.shape()[0] != in_geometry.gridLength() or - in.shape()[1] != in_geometry.gridWidth()) { - - std::string errmsg = "input signal data shape must match " - "input radar grid shape"; - throw InvalidArgument(ISCE_SRCINFO(), errmsg); - } - - std::complex* out_data = out.mutable_data(); - const std::complex* in_data = in.data(); - float* height_data = nullptr; - - if (height.has_value()) { - auto h = height.value(); - if (h.shape()[0] != out_geometry.gridLength() or - h.shape()[1] != out_geometry.gridWidth()) { - - std::string errmsg = "height array shape must match output " - "radar grid shape"; - throw InvalidArgument(ISCE_SRCINFO(), errmsg); - } - height_data = h.mutable_data(); - } - - DryTroposphereModel atm = parseDryTropoModel(dry_tropo_model); - - const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); - const auto g2rparams = parse_geo2rdr_params(geo2rdr_params); - - if (batch < 1) { - throw DomainError(ISCE_SRCINFO(), "batch size must be > 0"); - } - - ErrorCode err; - { - py::gil_scoped_release release; - err = backproject(out_data, out_geometry, in_data, in_geometry, - dem, fc, ds, kernel, atm, r2gparams, g2rparams, batch, - height_data); - } - // TODO bind ErrorCode class. For now return nonzero on failure. - return err != ErrorCode::Success; - }, - R"( - Focus in azimuth via time-domain backprojection. - )", - py::arg("out"), - py::arg("out_geometry"), - py::arg("in"), - py::arg("in_geometry"), - py::arg("dem"), - py::arg("fc"), - py::arg("ds"), - py::arg("kernel"), - py::arg("dry_tropo_model") = "tsx", - py::arg("rdr2geo_params") = py::dict(), - py::arg("geo2rdr_params") = py::dict(), - py::arg("batch") = 1024, - py::arg("height") = py::none()); -} diff --git a/python/extensions/pybind_isce3/cuda/focus/Backproject.cu b/python/extensions/pybind_isce3/cuda/focus/Backproject.cu new file mode 100644 index 000000000..6eec94b52 --- /dev/null +++ b/python/extensions/pybind_isce3/cuda/focus/Backproject.cu @@ -0,0 +1,368 @@ +#include "Backproject.h" +#include "pybind_isce3/focus/Backproject.h" // parse parameter dicts +#include "pybind_isce3/signal/NFFT2d.h" // parse NFFT2d parameters + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +using namespace isce3::cuda::focus; +using namespace isce3::except; + +using isce3::container::RadarGeometry; +using isce3::core::Kernel; +using isce3::error::ErrorCode; +using isce3::focus::parseDryTropoModel; +using isce3::geometry::DEMInterpolator; + + +void addbinding_cuda_backproject(py::module& m) +{ + m.def("backproject", []( + py::array_t, py::array::c_style> out, + const RadarGeometry& out_geometry, + py::array_t, py::array::c_style> in, + const RadarGeometry& in_geometry, + const DEMInterpolator& dem, + double fc, + double ds, + const Kernel& kernel, + const std::string& dry_tropo_model, + py::dict rdr2geo_params, + py::dict geo2rdr_params, + int batch, + std::optional> height) { + + if (out.ndim() != 2) { + throw InvalidArgument(ISCE_SRCINFO(), "output array must be 2-D"); + } + + if (out.shape()[0] != out_geometry.gridLength() or + out.shape()[1] != out_geometry.gridWidth()) { + + std::string errmsg = "output array shape must match output " + "radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + if (in.ndim() != 2) { + throw InvalidArgument(ISCE_SRCINFO(), "input signal data must be 2-D"); + } + + if (in.shape()[0] != in_geometry.gridLength() or + in.shape()[1] != in_geometry.gridWidth()) { + + std::string errmsg = "input signal data shape must match " + "input radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + std::complex* out_data = out.mutable_data(); + const std::complex* in_data = in.data(); + float* height_data = nullptr; + + if (height.has_value()) { + auto h = height.value(); + if (h.shape()[0] != out_geometry.gridLength() or + h.shape()[1] != out_geometry.gridWidth()) { + + std::string errmsg = "height array shape must match output " + "radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + height_data = h.mutable_data(); + } + + DryTroposphereModel atm = parseDryTropoModel(dry_tropo_model); + + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + const auto g2rparams = parse_geo2rdr_params(geo2rdr_params); + + if (batch < 1) { + throw DomainError(ISCE_SRCINFO(), "batch size must be > 0"); + } + + ErrorCode err; + { + py::gil_scoped_release release; + err = backproject(out_data, out_geometry, in_data, in_geometry, + dem, fc, ds, kernel, atm, r2gparams, g2rparams, batch, + height_data); + } + // TODO bind ErrorCode class. For now return nonzero on failure. + return err != ErrorCode::Success; + }, + R"( + Focus in azimuth via time-domain backprojection. + )", + py::arg("out"), + py::arg("out_geometry"), + py::arg("in"), + py::arg("in_geometry"), + py::arg("dem"), + py::arg("fc"), + py::arg("ds"), + py::arg("kernel"), + py::arg("dry_tropo_model") = "tsx", + py::arg("rdr2geo_params") = py::dict(), + py::arg("geo2rdr_params") = py::dict(), + py::arg("batch") = 1024, + py::arg("height") = py::none()); + + m.def("backproject_to_polar_grid", []( + const py::array_t, py::array::c_style> in, + const isce3::core::Linspace& in_slant_range, + const std::vector& pos, + const std::vector& vel, + const isce3::focus::PolarGrid& grid, + const DEMInterpolator& dem, + double fc, + const Kernel& kernel, + const std::string& dry_tropo_model, + py::dict rdr2geo_params) { + + if (in.ndim() != 2) { + throw InvalidArgument(ISCE_SRCINFO(), "input signal data must be 2-D"); + } + + if (in.shape()[0] != pos.size() or + in.shape()[1] != in_slant_range.size()) { + + std::string errmsg = "input signal data shape must match " + "input radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + DryTroposphereModel atm = parseDryTropoModel(dry_tropo_model); + + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + + const std::complex* in_data = in.data(); + + auto [err, outp, heightp] = [&]() { + py::gil_scoped_release release; + return isce3::cuda::focus::backprojectToPolarGrid(in_data, + in_slant_range, pos, vel, grid, + dem, fc, kernel, atm, r2gparams); + }(); + + // TODO bind ErrorCode class. For now return nonzero on failure. + bool status = err == ErrorCode::Success; + + auto out = move_to_numpy(std::move(outp), + {static_cast(grid.length()), + static_cast(grid.width())}); + auto height = move_to_numpy(std::move(heightp), + {static_cast(grid.length()), + static_cast(grid.width())}); + + return std::make_tuple(status, out, height); + }, + R"( + Focus in azimuth via time-domain backprojection. + )", + py::arg("in"), + py::arg("in_slant_range"), + py::arg("position"), + py::arg("velocity"), + py::arg("out_grid"), + py::arg("dem"), + py::arg("fc"), + py::arg("kernel"), + py::arg("dry_tropo_model") = "nodelay", // off here, on later + py::arg("rdr2geo_params") = py::dict()); + + m.def("project_polar_to_geo", []( + py::array_t>& geo_image, + const py::array_t& geo_points, + const isce3::focus::PolarGrid& grid, + const py::array_t>& polar_image, + const double wavelength, + py::dict nfft2_params) { + + // get root finding parameters + const auto params = parse_nfft2d_params(nfft2_params); + if (geo_points.size() != 3 * geo_image.size()) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "shape mismatch between geo image and position arrays"); + } + auto n = static_cast(geo_image.size()); + if ((polar_image.shape(0) != grid.length()) + or (polar_image.shape(1) != grid.width())) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "shape mismatch between polar image array and grid"); + } + + // XXX type cast after checking sizes, assume alignment is okay + // TODO redo with Eigen::Map or change interface from Vec3 to double[3]? + using isce3::core::Vec3; + static_assert(sizeof(Vec3) == (sizeof(double[3]))); + const auto ptr = reinterpret_cast(geo_points.data()); + + auto status = isce3::cuda::focus::projectPolarToGeo(geo_image.mutable_data(), ptr, + n, grid, polar_image.data(), wavelength, params); + + if (status != ErrorCode::Success) { + throw isce3::except::RuntimeError(ISCE_SRCINFO(), + "Could not compute map projection of polar grid coords."); + } + }, + py::arg("geo_image"), + py::arg("geo_points"), + py::arg("grid"), + py::arg("polar_image"), + py::arg("wavelength"), + py::arg("nfft2_params") = py::dict()); + + m.def("accumulate_polar_images_to_radar_grid", []( + py::array_t, py::array::c_style> out, + const RadarGeometry& out_geometry, + const isce3::core::Orbit& in_orbit, + const isce3::core::LUT2d& in_doppler, + const std::vector& grids, + const py::sequence& py_image_interpolators, + const DEMInterpolator& dem, + double fc, + double ds, + const std::string& dry_tropo_model, + py::dict rdr2geo_params, + py::dict geo2rdr_params, + std::optional> height) { + + if (out.ndim() != 2) { + throw InvalidArgument(ISCE_SRCINFO(), "output array must be 2-D"); + } + + if (out.shape()[0] != out_geometry.gridLength() or + out.shape()[1] != out_geometry.gridWidth()) { + + std::string errmsg = "output array shape must match output " + "radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + const auto nimg = py_image_interpolators.size(); + if (grids.size() != nimg) { + throw InvalidArgument(ISCE_SRCINFO(), "must have grid for each sub-image"); + } + + std::complex* out_data = out.mutable_data(); + float* height_data = nullptr; + + if (height.has_value()) { + auto h = height.value(); + if (h.shape()[0] != out_geometry.gridLength() or + h.shape()[1] != out_geometry.gridWidth()) { + + std::string errmsg = "height array shape must match output " + "radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + height_data = h.mutable_data(); + } + + DryTroposphereModel atm = parseDryTropoModel(dry_tropo_model); + + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + const auto g2rparams = parse_geo2rdr_params(geo2rdr_params); + + // Convert Python sequence to std::vector. Element type is + // pointer so this shouldn't involve any copy. Keep an owning + // reference to each element alive for the duration of the call, + // since the GIL is released below and another thread could + // otherwise mutate py_image_interpolators and drop the last + // reference to one of its elements while we hold a raw pointer + // into it. + using T = isce3::cuda::signal::NFFT2dResult; + std::vector interpolator_owners( + py_image_interpolators.begin(), py_image_interpolators.end()); + auto interpolators = std::vector(nimg); + std::transform(interpolator_owners.begin(), + interpolator_owners.end(), interpolators.begin(), + [](const py::handle& py_itp) -> const T* { + return &(py_itp.cast()); + }); + + ErrorCode err; + { + py::gil_scoped_release release; + err = isce3::cuda::focus::accumulatePolarImagesToRadarGrid( + out_data, out_geometry, in_orbit, in_doppler, grids, + interpolators, dem, fc, ds, atm, r2gparams, g2rparams, + height_data); + } + // TODO bind ErrorCode class. For now return nonzero on failure. + return err != ErrorCode::Success; + }, + py::arg("out"), + py::arg("out_geometry"), + py::arg("in_orbit"), + py::arg("in_doppler"), + py::arg("grids"), + py::arg("image_interpolators"), + py::arg("dem"), + py::arg("fc"), + py::arg("ds"), + py::arg("dry_tropo_model") = "tsx", + py::arg("rdr2geo_params") = py::dict(), + py::arg("geo2rdr_params") = py::dict(), + py::arg("height") = py::none()); + + m.def("merge_polar_images", []( + const std::vector& grids, + py::sequence py_image_interpolators, + const isce3::focus::PolarGrid& output_grid, + Eigen::Ref>> output_image, + const double fc, + const isce3::geometry::DEMInterpolator& dem, + const py::dict rdr2geo_params, // only difference for python + int az_block_size) + { + // Convert Python sequence to std::vector. Element type is + // pointer so this shouldn't involve any copy. Keep an owning + // reference to each element alive for the duration of the call, + // since the GIL is released below and another thread could + // otherwise mutate py_image_interpolators and drop the last + // reference to one of its elements while we hold a raw pointer + // into it. + const auto nimg = py_image_interpolators.size(); + using T = isce3::cuda::signal::NFFT2dResult; + std::vector interpolator_owners( + py_image_interpolators.begin(), py_image_interpolators.end()); + auto interpolators = std::vector(nimg); + std::transform(interpolator_owners.begin(), + interpolator_owners.end(), interpolators.begin(), + [](const py::handle& py_itp) -> const T* { + return &(py_itp.cast()); + }); + + const auto r2g_params = parse_rdr2geo_params(rdr2geo_params); + + { + py::gil_scoped_release release; + mergePolarImages(grids, interpolators, output_grid, + output_image, fc, dem, r2g_params, az_block_size); + } + }, + py::arg("grids"), + py::arg("image_interpolators"), + py::arg("output_grid"), + py::arg("output_image"), + py::arg("fc"), + py::arg("dem") = DEMInterpolator(), + py::arg("rdr2geo_parameters") = py::dict(), + py::arg("az_block_size") = 1024 + ); +} diff --git a/python/extensions/pybind_isce3/cuda/signal/NFFT2d.cu b/python/extensions/pybind_isce3/cuda/signal/NFFT2d.cu new file mode 100644 index 000000000..e30e7ff50 --- /dev/null +++ b/python/extensions/pybind_isce3/cuda/signal/NFFT2d.cu @@ -0,0 +1,56 @@ +#include "NFFT2d.h" +#include +#include +#include +#include + +using namespace isce3::cuda::signal; +namespace py = pybind11; + +template +void addbinding(py::class_>& pyNFFT2d) +{ + using dims_t = typename NFFT2d::dims_t; + pyNFFT2d + .def(py::init(), + py::arg("m"), py::arg("sizes"), py::arg("fft_sizes")) + ; + // TODO more methods +} + +template +void addbinding(py::class_>& pyNFFT2dResult) +{ + using dims_t = typename NFFT2d::dims_t; + pyNFFT2dResult + .def(py::init&>()) + .def("copy_to_host", [](const NFFT2dResult& self) { + return isce3::signal::NFFT2dResult(self); + }) + ; +} + +// instantiate +template void addbinding(py::class_>&); +template void addbinding(py::class_>&); +template void addbinding(py::class_>&); +template void addbinding(py::class_>&); + +void addbinding_make_image_nfft2d_gpu(pybind11::module& m) +{ + // TODO generalize to CF32 and CF64 + using T = float; + using array_t = isce3::core::EArray2D>; + m.def("make_image_nfft2d", [](Eigen::Ref image, + py::dict params, + bool pad_input) + { + const auto params_ = parse_nfft2d_params(params); + return isce3::cuda::signal::makeImageNFFT2d(image, params_, + pad_input); + }, + py::arg("image"), + py::arg("params") = py::dict{}, + py::arg("pad_input") = false) + ; +} \ No newline at end of file diff --git a/python/extensions/pybind_isce3/cuda/signal/NFFT2d.h b/python/extensions/pybind_isce3/cuda/signal/NFFT2d.h new file mode 100644 index 000000000..82b4cf587 --- /dev/null +++ b/python/extensions/pybind_isce3/cuda/signal/NFFT2d.h @@ -0,0 +1,12 @@ +#pragma once + +#include +#include + +template +void addbinding(pybind11::class_>& pyNFFT2d); + +template +void addbinding(pybind11::class_>& pyNFFT2dResult); + +void addbinding_make_image_nfft2d_gpu(pybind11::module& m); \ No newline at end of file diff --git a/python/extensions/pybind_isce3/cuda/signal/signal.cpp b/python/extensions/pybind_isce3/cuda/signal/signal.cpp index fa7bdde29..4da394773 100644 --- a/python/extensions/pybind_isce3/cuda/signal/signal.cpp +++ b/python/extensions/pybind_isce3/cuda/signal/signal.cpp @@ -1,6 +1,7 @@ #include "signal.h" #include "Crossmul.h" +#include "NFFT2d.h" namespace py = pybind11; @@ -10,7 +11,16 @@ void addsubmodule_cuda_signal(py::module & m) // forward declare bound classes py::class_ pyCrossmul(m_signal, "Crossmul"); + py::class_> pyNFFT2dF32(m_signal, "NFFT2dF32"); + py::class_> pyNFFT2dF64(m_signal, "NFFT2dF64"); + py::class_> pyNFFT2dF32Result(m_signal, "NFFT2dF32Result"); + py::class_> pyNFFT2dF64Result(m_signal, "NFFT2dF64Result"); // add bindings addbinding(pyCrossmul); + addbinding(pyNFFT2dF32); + addbinding(pyNFFT2dF64); + addbinding(pyNFFT2dF32Result); + addbinding(pyNFFT2dF64Result); + addbinding_make_image_nfft2d_gpu(m_signal); } diff --git a/python/extensions/pybind_isce3/focus/Backproject.cpp b/python/extensions/pybind_isce3/focus/Backproject.cpp index 1f62c6f42..381e3a373 100644 --- a/python/extensions/pybind_isce3/focus/Backproject.cpp +++ b/python/extensions/pybind_isce3/focus/Backproject.cpp @@ -1,17 +1,22 @@ #include "Backproject.h" +#include "pybind_isce3/signal/NFFT2d.h" // parse NFFT2d parameters +#include #include +#include #include #include #include #include +#include #include #include #include #include #include #include +#include namespace py = pybind11; @@ -19,11 +24,13 @@ using namespace isce3::focus; using isce3::container::RadarGeometry; using isce3::core::Kernel; +using isce3::core::EArray2D; using isce3::error::ErrorCode; using isce3::except::InvalidArgument; using isce3::geometry::DEMInterpolator; using isce3::geometry::detail::Rdr2GeoBracketParams; using isce3::geometry::detail::Geo2RdrBracketParams; +using isce3::signal::NFFT2dResult; Rdr2GeoBracketParams parse_rdr2geo_params(const py::dict& params) @@ -78,6 +85,74 @@ Geo2RdrBracketParams parse_geo2rdr_params(const py::dict& params) } +void addbinding(py::class_& pyPolarGrid) +{ + using isce3::core::Vec3; + using isce3::core::Linspace; + using isce3::core::LookSide; + + pyPolarGrid + .def(py::init, Linspace, LookSide>(), + py::arg("aztime_start"), + py::arg("aztime_end"), + py::arg("origin"), + py::arg("axis"), + py::arg("range"), + py::arg("sin_squint"), + py::arg("look_side") + ) + .def_readonly("aztime_start", &PolarGrid::aztime_start) + .def_readonly("aztime_end", &PolarGrid::aztime_end) + .def_readonly("origin", &PolarGrid::origin) + .def_readonly("axis", &PolarGrid::axis) + .def_readonly("range", &PolarGrid::range) + .def_readonly("sin_squint", &PolarGrid::sin_squint) + .def_readonly("look_side", &PolarGrid::look_side) + .def_property_readonly("shape", [](const PolarGrid& self) { + return std::make_tuple(self.sin_squint.size(), self.range.size()); + }) + .def("__repr__", [](const py::object self) { + std::vector keys {"aztime_start", "aztime_end", + "origin", "axis", "range", "sin_squint", "look_side"}; + std::string out("PolarGrid("); + for (auto it = keys.begin(); it != keys.end(); ++it) { + auto key = *it; + auto ckey = key.c_str(); + out += key + "=" + std::string(py::str(self.attr(ckey))); + if (it != keys.end() - 1) + out += ", "; + } + return out + ")"; + }) + // all properties are read-only, so instances are hashable + .def_property_readonly("_members", [](const py::object self) { + const auto q = py::getattr(self, "sin_squint").cast>(); + const auto r = py::getattr(self, "range").cast>(); + return py::make_tuple( + py::getattr(self, "aztime_start"), + py::getattr(self, "aztime_end"), + py::tuple(py::getattr(self, "origin")), + py::tuple(py::getattr(self, "axis")), + py::make_tuple(q.first(), q.spacing(), q.size()), + py::make_tuple(r.first(), r.spacing(), r.size()), + py::getattr(self, "look_side") + ); + }) + .def("__hash__", [](const py::object self) { + return py::hash(py::getattr(self, "_members")); + }) + .def("__eq__", [](const py::object self, const PolarGrid& typed_other) { + // Strongly-typed function signature means we don't have to check + // type. The _members method is only defined in Python, though, so + // cast to py::object. + const py::object other = py::cast(typed_other); + const auto a = getattr(self, "_members"); + const auto b = getattr(other, "_members"); + return a.equal(b); + }) + ; +} + void addbinding_backproject(py::module& m) { m.def("backproject", []( @@ -151,6 +226,40 @@ void addbinding_backproject(py::module& m) }, R"( Focus in azimuth via time-domain backprojection. + + Parameters + ---------- + out : numpy.ndarray[complex64] + Output 2D array of focused signal data. + out_geometry : isce3.container.RadarGeometry + Target output grid, orbit, and Doppler. + in : numpy.ndarray[complex64] + Input 2D array of range-compressed signal data. + in_geometry : isce3.container.RadarGeometry + Input data grid, orbit, and Doppler. + dem : isce3.geometry.DEMInterpolator + Digital elevation model. + fc : float + Radar center frequency (Hz). + ds : float + Desired azimuth resolution (m). + kernel : isce3.core.Kernel + 1-D interpolation kernel. + dry_tropo_model : str, optional + Dry troposphere path delay model (defaults to "tsx"). + rdr2geo_params : dict, optional + rdr2geo_bracket configuration keyword arguments. + geo2rdr_params : dict, optional + geo2rdr_bracket configuration keyword arguments. + height : numpy.ndarray[float32], optional + Output array to store height of each pixel in meters above + the ellipsoid. + + Returns + ------- + bool + True if successful, False if geometry fails to converge for + any pixel (those pixels are set to NaN). )", py::arg("out"), py::arg("out_geometry"), @@ -164,4 +273,636 @@ void addbinding_backproject(py::module& m) py::arg("rdr2geo_params") = py::dict(), py::arg("geo2rdr_params") = py::dict(), py::arg("height") = py::none()); + + m.def("setup_polar_grid_for_pulses", &setupPolarGridForPulses, + R"( + Setup a polar (range-Doppler) grid corresponding to a set of pulses. + + Parameters + ---------- + in_geometry : isce3.container.RadarGeometry + azimuth_time : Sequence[float] + range_bandwidth : float + azimuth_resolution : float + oversample_range : float, optional + oversample_azimuth : float, optional + num_doppler_eval : int, optional + Number of points across swath to evaluate Doppler centroid to + bound the variation of the centroid. Default = 2 + pri : float, optional + Pulse repetition interval in s. If variable, provide the PRI + between the last pulse and the next one. If not provided the + average PRI will be used. + + Returns + ------- + polar_grid : isce3.focus.PolarGrid + Polar grid that efficiently samples the raw data. + position : list[numpy.ndarray] + Sensor position at each input pulse time. + velocity : list[numpy.ndarray] + Sensor velocity at each input pulse time. + )", + py::arg("in_geometry"), + py::arg("azimuth_time"), + py::arg("range_bandwidth"), + py::arg("azimuth_resolution"), + py::arg("oversample_range") = 1.2, + py::arg("oversample_azimuth") = 1.2, + py::arg("num_doppler_eval") = 2, + py::arg("pri") = py::none()); + + m.def("get_polar_angle_time_constant", &getPolarAngleTimeConstant, + R"( + Get the time constant associated with polar angle spacing + + Parameters + ---------- + fc : float + Radar center frequency, Hz + vs : float + Satellite velocity (along azimuth axis), m/s + bandwidth : float, optional + Radar bandwidth, Hz (defaults to zero, e.g., narrow band) + c : float, optional + Speed of light, m/s (defaults to vacuum sol) + + Returns + ------- + tq : float + Time constant $T_q$, s + + This time constant is used to determine the sampling requirement for the + sine of the squint angle (dimensionless Doppler) + $$ q = \frac{\vec{v}}{v} \cdot \hat{l} $$ + where $\vec{v}$ is the velocity and $\hat{l}$ is the line-of-sight direction. + Specifically, the Nyquist criterion is + $$ \Delta q \leq \frac{T_q}{T_{sa}} $$ + where $T_{sa}$ is the time duration of the synthetic aperture. + + Helps implement equation (11) in @cite yegulalp2013 + )", + py::arg("fc"), + py::arg("vs"), + py::arg("bandwidth") = 0.0, + py::arg("c") = isce3::core::speed_of_light); + + m.def("backproject_to_polar_grid", []( + const py::array_t, py::array::c_style> in, + const isce3::core::Linspace& in_slant_range, + const std::vector& pos, + const std::vector& vel, + const PolarGrid& grid, + const DEMInterpolator& dem, + double fc, + const Kernel& kernel, + const std::string& dry_tropo_model, + py::dict rdr2geo_params) { + + if (in.ndim() != 2) { + throw InvalidArgument(ISCE_SRCINFO(), "input signal data must be 2-D"); + } + + if (in.shape()[0] != pos.size() or + in.shape()[1] != in_slant_range.size()) { + + std::string errmsg = "input signal data shape must match " + "input radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + if (pos.size() != vel.size()) { + throw InvalidArgument(ISCE_SRCINFO(), "must provide same " + "number of position and velocity vectors"); + } + + DryTroposphereModel atm = parseDryTropoModel(dry_tropo_model); + + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + + const std::complex* in_data = in.data(); + + auto [err, outp, heightp] = [&]() { + py::gil_scoped_release release; + return isce3::focus::backprojectToPolarGrid(in_data, + in_slant_range, pos, vel, grid, + dem, fc, kernel, atm, r2gparams); + }(); + + // TODO bind ErrorCode class. For now return nonzero on failure. + bool status = err == ErrorCode::Success; + + auto out = move_to_numpy(std::move(outp), + {static_cast(grid.length()), + static_cast(grid.width())}); + auto height = move_to_numpy(std::move(heightp), + {static_cast(grid.length()), + static_cast(grid.width())}); + + return std::make_tuple(status, out, height); + }, + R"( + Focus in azimuth via time-domain backprojection onto a + polar grid. + + Parameters + ---------- + in : numpy.ndarray[complex64] + Input 2D array of range-compressed signal data. + in_slant_range : isce3.core.Linspace + Slant range grid of the input data (m). + position : list[numpy.ndarray] + Platform position vectors at each pulse (ECEF, m). + velocity : list[numpy.ndarray] + Platform velocity vectors at each pulse (ECEF, m/s). + out_grid : isce3.focus.PolarGrid + Target polar grid to backproject onto. + dem : isce3.geometry.DEMInterpolator + Digital elevation model. + fc : float + Radar center frequency (Hz). + kernel : isce3.core.Kernel + 1D interpolation kernel. + dry_tropo_model : str, optional + Dry troposphere path delay model (defaults to "nodelay"). + rdr2geo_params : dict, optional + rdr2geo_bracket configuration keyword arguments. + + Returns + ------- + tuple + - success : bool + True if successful, False if geometry fails + to converge for any pixel. + - out : numpy.ndarray[complex64] + Focused signal data on the polar grid. + - height : numpy.ndarray[float32] + Per-pixel height above the ellipsoid (m). + )", + py::arg("in"), + py::arg("in_slant_range"), + py::arg("position"), + py::arg("velocity"), + py::arg("out_grid"), + py::arg("dem"), + py::arg("fc"), + py::arg("kernel"), + py::arg("dry_tropo_model") = "nodelay", // off here, on later + py::arg("rdr2geo_params") = py::dict()); + + m.def("merge_polar_grids", [](const std::vector& grids, + const DEMInterpolator& dem, + py::dict rdr2geo_params, + const std::optional& dq_min, + const std::optional& tq) { + const auto r2g_params = parse_rdr2geo_params(rdr2geo_params); + return mergePolarGrids(grids, dem, r2g_params, dq_min, tq); + }, + R"( + Create polar grid capable of sampling data from all input grids. + + Parameters + ---------- + grids : list[isce3.focus.PolarGrid] + List of subaperture grids. + dem : isce3.geometry.DEMInterpolator, optional + Digital elevation model reporting height (m) above the + ellipsoid associated with its CRS. + rdr2geo_params : dict, optional + rdr2geo_bracket configuration keyword arguments. + dq_min : float, optional + Minimum allowed dimensionless Doppler spacing. + Necessary for stripmap processing large subapertures. + tq : float, optional + Time constant for dimensionless Doppler spacing. + If not provided it will be inferred from input grids. + + Returns + ------- + isce3.focus.PolarGrid + Merged output polar grid. + )", + py::arg("grids"), + py::arg("dem") = DEMInterpolator(), + py::arg("rdr2geo_params") = py::dict(), + py::arg("dq_min") = py::none(), + py::arg("tq") = py::none() + ); + + m.def("merge_polar_images", []( + const std::vector& grids, + const py::sequence& py_image_interpolators, + const PolarGrid& output_grid, + Eigen::Ref>> output_image, + const double fc, + const isce3::geometry::DEMInterpolator& dem, + const py::dict rdr2geo_params, // only difference for python + int az_block_size) + { + const auto r2g_params = parse_rdr2geo_params(rdr2geo_params); + // Convert Python sequence to std::vector of pointers + const auto nimg = py_image_interpolators.size(); + using T = NFFT2dResult; + auto interpolators = std::vector(nimg); + std::transform(py_image_interpolators.begin(), + py_image_interpolators.end(), interpolators.begin(), + [](const py::handle& py_itp) -> const T* { + return &(py_itp.cast()); + }); + return mergePolarImages(grids, interpolators, output_grid, + output_image, fc, dem, r2g_params, az_block_size); + }, + R"( + Merge subaperture polar grid images into a single output grid. + + Combines multiple subaperture polar grid images onto a merged + output polar grid. For each pixel in the output grid, the 3D + target position is computed via polar2geo, and the input image + data is accumulated via NFFT-based interpolation. The output + image is expected to be zero-initialized by the caller. + + Parameters + ---------- + grids : list[isce3.focus.PolarGrid] + List of subaperture input polar grids. + image_interpolators : list[numpy.ndarray] + NFFT interpolators for each input grid. + output_grid : isce3.focus.PolarGrid + Merged output polar grid. + output_image : numpy.ndarray[complex64] + Accumulated output image (must be zero-initialized); + dimensions must match output_grid. + fc : float + Center frequency (Hz). + dem : isce3.geometry.DEMInterpolator, optional + Digital elevation model. + rdr2geo_parameters : dict, optional + rdr2geo_bracket configuration keyword arguments. + az_block_size : int, optional + Number of azimuth rows to process at a time + (defaults to 1024). + + Raises + ------ + isce3.except.LengthError + If output image dimensions or grid/interpolator + counts are inconsistent. + isce3.except.InvalidArgument + If look directions are inconsistent or + az_block_size is negative. + isce3.except.DomainError + If polar2geo fails to converge. + isce3.except.RuntimeError + If NFFT interpolation fails. + )", + py::arg("grids"), + py::arg("image_interpolators"), + py::arg("output_grid"), + py::arg("output_image"), + py::arg("fc"), + py::arg("dem") = DEMInterpolator(), + py::arg("rdr2geo_parameters") = py::dict(), + py::arg("az_block_size") = 1024 + ); + + m.def("accumulate_polar_images_to_radar_grid", []( + py::array_t, py::array::c_style> out, + const RadarGeometry& out_geometry, + const isce3::core::Orbit& in_orbit, + const isce3::core::LUT2d& in_doppler, + const std::vector& grids, + const py::sequence& py_image_interpolators, + const DEMInterpolator& dem, + double fc, + double ds, + const std::string& dry_tropo_model, + py::dict rdr2geo_params, + py::dict geo2rdr_params, + std::optional> height) { + + if (out.ndim() != 2) { + throw InvalidArgument(ISCE_SRCINFO(), "output array must be 2-D"); + } + + if (out.shape()[0] != out_geometry.gridLength() or + out.shape()[1] != out_geometry.gridWidth()) { + + std::string errmsg = "output array shape must match output " + "radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + + const auto nimg = py_image_interpolators.size(); + if (grids.size() != nimg) { + throw InvalidArgument(ISCE_SRCINFO(), "must have grid for each sub-image"); + } + + std::complex* out_data = out.mutable_data(); + float* height_data = nullptr; + + if (height.has_value()) { + auto h = height.value(); + if (h.shape()[0] != out_geometry.gridLength() or + h.shape()[1] != out_geometry.gridWidth()) { + + std::string errmsg = "height array shape must match output " + "radar grid shape"; + throw InvalidArgument(ISCE_SRCINFO(), errmsg); + } + height_data = h.mutable_data(); + } + + DryTroposphereModel atm = parseDryTropoModel(dry_tropo_model); + + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + const auto g2rparams = parse_geo2rdr_params(geo2rdr_params); + + // Convert Python sequence to std::vector of pointers + using T = NFFT2dResult; + auto interpolators = std::vector(nimg); + std::transform(py_image_interpolators.begin(), + py_image_interpolators.end(), interpolators.begin(), + [](const py::handle& py_itp) -> const T* { + return &(py_itp.cast()); + }); + + ErrorCode err; + { + py::gil_scoped_release release; + err = accumulatePolarImagesToRadarGrid(out_data, out_geometry, + in_orbit, in_doppler, grids, interpolators, dem, fc, + ds, atm, r2gparams, g2rparams, height_data); + } + // TODO bind ErrorCode class. For now return nonzero on failure. + return err != ErrorCode::Success; + }, + R"( + Accumulate polar grid images onto an output stripmap radar grid. + + Combines multiple subaperture polar grid images together onto a + stripmap radar geometry grid. For each pixel in the output grid, + the target position is computed via rdr2geo, the corresponding + coherent processing interval is determined via geo2rdr, and the + polar image data is accumulated using NFFT-based interpolation. + + Parameters + ---------- + out : numpy.ndarray[complex64] + Output 2D array of focused signal data. + out_geometry : isce3.container.RadarGeometry + Target output grid, orbit, and Doppler. + in_orbit : isce3.core.Orbit + Input data orbit. + in_doppler : isce3.core.LUT2d + Input data Doppler centroid LUT. + grids : list[isce3.focus.PolarGrid] + List of subaperture polar grids. + image_interpolators : list[numpy.ndarray] + NFFT interpolators for each polar grid. + dem : isce3.geometry.DEMInterpolator, optional + Digital elevation model. + fc : float + Center frequency (Hz). + ds : float + Desired azimuth resolution (m). + dry_tropo_model : str, optional + Dry troposphere path delay model (defaults to "tsx"). + rdr2geo_params : dict, optional + rdr2geo_bracket configuration keyword arguments. + geo2rdr_params : dict, optional + geo2rdr_bracket configuration keyword arguments. + height : numpy.ndarray[float32], optional + Output array to store height of each pixel in meters + above the ellipsoid. + + Returns + ------- + bool + True if successful, False if rdr2geo or geo2rdr fails + to converge for any pixel (those pixels are set to + NaN). + )", + py::arg("out"), + py::arg("out_geometry"), + py::arg("in_orbit"), + py::arg("in_doppler"), + py::arg("grids"), + py::arg("image_interpolators"), + py::arg("dem"), + py::arg("fc"), + py::arg("ds"), + py::arg("dry_tropo_model") = "tsx", + py::arg("rdr2geo_params") = py::dict(), + py::arg("geo2rdr_params") = py::dict(), + py::arg("height") = py::none()); + + m.def("find_polar_grid_bbox_in_radar_grid", []( + const PolarGrid& polar_grid, + const RadarGeometry& radar_geom, + const DEMInterpolator& dem, + py::dict rdr2geo_params, + py::dict geo2rdr_params, + int nextra) { + + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + const auto g2rparams = parse_geo2rdr_params(geo2rdr_params); + + auto [i0, i1, j0, j1, status] = findPolarGridBoundingBoxInRadarGrid( + polar_grid, radar_geom, dem, r2gparams, g2rparams, nextra); + + if (status != ErrorCode::Success) { + throw isce3::except::RuntimeError(ISCE_SRCINFO(), + "Could not determine polar grid bounds within radar grid."); + } + auto rows = py::slice( + static_cast(i0), + static_cast(i1), + std::nullopt); + auto cols = py::slice( + static_cast(j0), + static_cast(j1), + std::nullopt); + return std::make_tuple(rows, cols); + }, + R"( + Find the subset of a radar grid covered by a polar grid. + + Computes the bounding box of the polar grid in stripmap radar + coordinates, then converts this bounding box to integer radar + grid indices (azimuth line, range sample). If the polar grid + does not overlap the radar grid at all, a zero-sized subset + (0, 0, 0, 0) is returned. + + Parameters + ---------- + polar_grid : isce3.focus.PolarGrid + Input polar grid. + radar_geom : isce3.container.RadarGeometry + Target radar geometry grid. + dem : isce3.geometry.DEMInterpolator + Digital elevation model. + rdr2geo_params : dict, optional + rdr2geo_bracket configuration keyword arguments. + geo2rdr_params : dict, optional + geo2rdr_bracket configuration keyword arguments. + nextra : int, optional + Number of extra perimeter points per edge (defaults to 0). + + Returns + ------- + tuple[slice, slice] + - rows : slice + Azimuth line range (start, end). + - cols : slice + Range sample range (start, end). + + Raises + ------ + isce3.except.RuntimeError + If the polar grid bounds cannot be determined. + )", + py::arg("polar_grid"), + py::arg("radar_geom"), + py::arg("dem"), + py::arg("rdr2geo_params") = py::dict(), + py::arg("geo2rdr_params") = py::dict(), + py::arg("nextra") = 0); + + m.def("computeRadarGridGeoPoints", []( + const RadarGeometry& geom, + const DEMInterpolator& dem, + py::dict rdr2geo_params) { + + // get root finding parameters + const auto r2gparams = parse_rdr2geo_params(rdr2geo_params); + + // allocate memory + const py::ssize_t m = geom.gridLength(), n = geom.gridWidth(); + auto points = py::array_t({m, n, 3L}); + + // XXX type cast after checking sizes, assume alignment is okay + // TODO redo with Eigen::Map or change interface from Vec3 to double[3]? + using isce3::core::Vec3; + static_assert(sizeof(Vec3) == (sizeof(double[3]))); + auto ptr = reinterpret_cast(points.mutable_data()); + + // run the thing + auto status = computeRadarGridGeoPoints(ptr, geom, dem, r2gparams); + + if (status != ErrorCode::Success) { + throw isce3::except::RuntimeError(ISCE_SRCINFO(), + "Could not compute map projection of polar grid coords."); + } + return points; + }, + R"( + Compute 3D geo coordinates for a radar grid. + + Computes the 3D XYZ position for every pixel in a radar + geometry grid using rdr2geo_bracket. + + Parameters + ---------- + geom : isce3.container.RadarGeometry + Radar geometry grid, orbit, and Doppler. + dem : isce3.geometry.DEMInterpolator + Digital elevation model. + rdr2geo_params : dict, optional + rdr2geo_bracket configuration keyword arguments. + + Returns + ------- + numpy.ndarray[float64] + 3D XYZ positions (ECEF, m) with shape (m, n, 3), where m is the + grid length and n is the grid width. + + Raises + ------ + isce3.except.RuntimeError + If rdr2geo fails to converge for any pixel. + )", + py::arg("geom"), + py::arg("dem"), + py::arg("rdr2geo_params") = py::dict()); + + m.def("accumulate_polar_image_to_geo_points", []( + py::array_t, py::array::c_style>& image, + const py::array_t& xyz, + const PolarGrid& grid, + const NFFT2dResult& nfft, + const double wavelength, + const std::optional>& mask) { + + const auto n = image.size(); + if (xyz.size() != 3 * n) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "shape mismatch between geo image and position arrays"); + } + if (xyz.shape(xyz.ndim() - 1) != 3) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "expected trailing dimension size == 3 for XYZ points"); + } + if (mask.has_value() and (mask.value().size() != n)) { + throw isce3::except::LengthError(ISCE_SRCINFO(), + "pixel mask size does not equal image size"); + } + + // XXX type cast after checking sizes + using isce3::core::Vec3; + static_assert(sizeof(Vec3) == (sizeof(double[3]))); + const auto ptr = reinterpret_cast(xyz.data()); + + std::optional mask_ptr = std::nullopt; + if (mask.has_value()) { + mask_ptr = mask.value().data(); + } + + auto status = accumulatePolarImageToGeoPoints( + image.mutable_data(), ptr, n, grid, nfft, wavelength, mask_ptr); + + if (status != ErrorCode::Success) { + throw isce3::except::RuntimeError(ISCE_SRCINFO(), + "Could not compute map projection of polar grid coords."); + } + }, + R"( + Interpolate a polar grid image to given XYZ positions. + + Accumulates (adds) contributions from a polar grid image into + an output complex signal array at specified 3D positions. For + each position, the target location in the polar grid is + computed via geo2polar, and the image is interpolated using + NFFT. The phase is compensated by the wavenumber-range + product kw * range. + + Parameters + ---------- + image : numpy.ndarray[complex64] + Output complex signal data (accumulates, so caller + must init to zero). + xyz : numpy.ndarray[float64] + Target 3D positions (ECEF, m) with shape (n, 3). + grid : isce3.focus.PolarGrid + Polar grid containing the image data. + nfft : numpy.ndarray + NFFT interpolator for the polar grid. + wavelength : float + Radar wavelength (m). + mask : numpy.ndarray[bool], optional + Pixel mask; pixels with false are skipped. + + Raises + ------ + isce3.except.LengthError + If shape mismatch between geo image and position + arrays, or if mask size does not equal image size. + isce3.except.RuntimeError + If the computation fails. + )", + py::arg("image"), + py::arg("xyz"), + py::arg("grid"), + py::arg("nfft"), + py::arg("wavelength"), + py::arg("mask") = py::none()); } diff --git a/python/extensions/pybind_isce3/focus/Backproject.h b/python/extensions/pybind_isce3/focus/Backproject.h index 37ab65616..3630c8003 100644 --- a/python/extensions/pybind_isce3/focus/Backproject.h +++ b/python/extensions/pybind_isce3/focus/Backproject.h @@ -1,9 +1,14 @@ #pragma once +#include +#include #include +#include +#include #include #include +void addbinding(pybind11::class_& pyPolarGrid); void addbinding_backproject(pybind11::module& m); isce3::geometry::detail::Rdr2GeoBracketParams @@ -11,3 +16,34 @@ parse_rdr2geo_params(const pybind11::dict& params); isce3::geometry::detail::Geo2RdrBracketParams parse_geo2rdr_params(const pybind11::dict& params); + +/** + * Transfer ownership of a unique_ptr buffer to a numpy array without + * copying. The array takes ownership via a capsule whose destructor deletes[] + * the buffer when the last reference dies. + * + * \param ptr Unique pointer to transfer (will be released) + * \param shape Array shape (C-order / row-major is assumed) + * \return numpy array wrapping the transferred buffer + */ +template +pybind11::array_t move_to_numpy( + std::unique_ptr ptr, + std::vector shape) +{ + // Compute C-order strides. + std::vector strides(shape.size()); + pybind11::ssize_t stride = sizeof(T); + for (int i = static_cast(shape.size()) - 1; i >= 0; --i) { + strides[i] = stride; + stride *= shape[i]; + } + + T* data = ptr.release(); + + pybind11::capsule owner(data, [](void *p) { + delete[] static_cast(p); + }); + + return pybind11::array_t(shape, strides, data, owner); +} diff --git a/python/extensions/pybind_isce3/focus/focus.cpp b/python/extensions/pybind_isce3/focus/focus.cpp index 436d833e5..37bfb9ec2 100644 --- a/python/extensions/pybind_isce3/focus/focus.cpp +++ b/python/extensions/pybind_isce3/focus/focus.cpp @@ -15,6 +15,7 @@ void addsubmodule_focus(py::module & m) // forward declare bound enums py::enum_ pyDryTropoModel(m_focus, "DryTroposphereModel"); + py::class_ pyPolarGrid(m_focus, "PolarGrid"); py::class_ pyRangeComp(m_focus, "RangeComp"); py::enum_ pyMode(pyRangeComp, "Mode"); @@ -27,4 +28,5 @@ void addsubmodule_focus(py::module & m) addbinding_tsx_delay(m_focus); addbindings_presum(m_focus); addbinding(pyRangeComp); + addbinding(pyPolarGrid); } diff --git a/python/extensions/pybind_isce3/geogrid/geogrid.cpp b/python/extensions/pybind_isce3/geogrid/geogrid.cpp index 95a53923a..42cd2271a 100644 --- a/python/extensions/pybind_isce3/geogrid/geogrid.cpp +++ b/python/extensions/pybind_isce3/geogrid/geogrid.cpp @@ -1,6 +1,7 @@ #include "geogrid.h" #include "getRadarGrid.h" #include "relocateRaster.h" +#include "geogrid_ecef_coords.h" void addsubmodule_geogrid(py::module & m) { @@ -8,4 +9,5 @@ void addsubmodule_geogrid(py::module & m) addbinding_get_radar_grid(m_geogrid); addbinding_relocate_raster(m_geogrid); + addbinding_get_geogrid_ecef_coords(m_geogrid); } diff --git a/python/extensions/pybind_isce3/geogrid/geogrid_ecef_coords.cpp b/python/extensions/pybind_isce3/geogrid/geogrid_ecef_coords.cpp new file mode 100644 index 000000000..dca5e1182 --- /dev/null +++ b/python/extensions/pybind_isce3/geogrid/geogrid_ecef_coords.cpp @@ -0,0 +1,57 @@ +#include "geogrid_ecef_coords.h" + +#include +#include +#include +#include + +#include + +namespace py = pybind11; +using isce3::geometry::DEMInterpolator; +using isce3::product::GeoGridParameters; +using isce3::core::makeProjection; +using isce3::except::InvalidArgument; + +void addbinding_get_geogrid_ecef_coords(pybind11::module& m) +{ + m.def("get_geogrid_ecef_coords", + [](const GeoGridParameters& grid, const DEMInterpolator& dem) + { + const long m = grid.length(); + const long n = grid.width(); + + auto proj_in = makeProjection(grid.epsg()); + auto proj_out = makeProjection(4978); + + if (not dem.haveStats()) { + throw InvalidArgument(ISCE_SRCINFO(), + "Input DEM does not have stats."); + } + + const double h = dem.meanHeight(); + + auto out = py::array_t({m, n, 3L}); + auto r = out.mutable_unchecked<3>(); + + #pragma omp parallel for collapse(2) + for (long i = 0; i < m; ++i) { + for (long j = 0; j < n; ++j) { + isce3::core::Vec3 pos, xyz; + pos[0] = grid.startX() + grid.spacingX() * j; + pos[1] = grid.startY() + grid.spacingY() * i; + pos[2] = h; + isce3::core::projTransform(proj_in.get(), proj_out.get(), + pos, xyz); + #pragma unroll + for (int k = 0; k < 3; ++k) { + r(i, j, k) = xyz[k]; + } + } + } + return out; + }, + py::arg("grid"), + py::arg("dem") + ); +} diff --git a/python/extensions/pybind_isce3/geogrid/geogrid_ecef_coords.h b/python/extensions/pybind_isce3/geogrid/geogrid_ecef_coords.h new file mode 100644 index 000000000..adefb0d71 --- /dev/null +++ b/python/extensions/pybind_isce3/geogrid/geogrid_ecef_coords.h @@ -0,0 +1,5 @@ +#pragma once + +#include + +void addbinding_get_geogrid_ecef_coords(pybind11::module& m); diff --git a/python/extensions/pybind_isce3/geometry/rdr2geo_roots.cpp b/python/extensions/pybind_isce3/geometry/rdr2geo_roots.cpp index b2dd0404f..0a2d6d30d 100644 --- a/python/extensions/pybind_isce3/geometry/rdr2geo_roots.cpp +++ b/python/extensions/pybind_isce3/geometry/rdr2geo_roots.cpp @@ -6,14 +6,17 @@ #include #include #include +#include #include #include +#include #include namespace py = pybind11; using namespace isce3::core; using namespace isce3::geometry; +using isce3::geometry::detail::polar2geo_bracket; void addbinding_rdr2geo_roots(py::module& m) { @@ -91,4 +94,35 @@ void addbinding_rdr2geo_roots(py::module& m) of the nadir vector into a plane perpendicular to the velocity. For simplicity, we use the geocentric nadir definition. )"); + + + m.def("polar2geo_bracket", + [](const Vec3& origin, const Vec3& axis, const double slant_range, + const double sin_squint, + py::object py_side, + const DEMInterpolator& dem, + double tol_height, + double look_min, double look_max) { + Vec3 target_xyz; + double look_angle; + const auto side = duck_look_side(py_side); + const auto csq = std::sqrt(1.0 - sin_squint * sin_squint); + const auto ellipsoid = makeProjection(dem.epsgCode())->ellipsoid(); + auto ec = polar2geo_bracket(&target_xyz, &look_angle, origin, axis, + slant_range, sin_squint, csq, dem, ellipsoid, side, + {tol_height, look_min, look_max}); + if (ec != isce3::error::ErrorCode::Success) { + throw std::runtime_error("failed to converge"); + } + return std::make_tuple(target_xyz, look_angle); + }, + py::arg("origin"), + py::arg("axis"), + py::arg("slant_range"), + py::arg("sin_squint"), + py::arg("side"), + py::arg("dem") = DEMInterpolator(), + py::arg("tol_height") = isce3::geometry::detail::DEFAULT_TOL_HEIGHT, + py::arg("look_min") = 0.0, + py::arg("look_max") = M_PI / 2); } diff --git a/python/extensions/pybind_isce3/signal/NFFT2d.cpp b/python/extensions/pybind_isce3/signal/NFFT2d.cpp new file mode 100644 index 000000000..7d7210db7 --- /dev/null +++ b/python/extensions/pybind_isce3/signal/NFFT2d.cpp @@ -0,0 +1,106 @@ +#include "NFFT2d.h" +#include +#include +#include + +using namespace isce3::signal; +namespace py = pybind11; + +NFFT2dParams parse_nfft2d_params(const py::dict& params) +{ + auto parse_ms = [](const py::dict& d) { + NFFTParams out; + for (auto item : d) { + auto key = item.first.cast(); + if (key == "m") { + out.m = item.second.cast(); + } + else if (key == "s") { + out.s = item.second.cast(); + } + else { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "unexpected NFFT keyword: " + key); + } + } + return out; + }; + NFFT2dParams out; + for (auto item : params) { + auto key = item.first.cast(); + if (key == "rows") { + out.rows = parse_ms(item.second.cast()); + } + else if (key == "cols") { + out.cols = parse_ms(item.second.cast()); + } + else { + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), + "unexpected NFFT2dParms keyword: " + key); + } + } + return out; +} + +template +void addbinding(py::class_>& pyNFFT2d) +{ + using dims_t = typename NFFT2d::dims_t; + pyNFFT2d + .def(py::init(), + py::arg("m"), py::arg("sizes"), py::arg("fft_sizes")) + .def_property_readonly("spectrum", [](const NFFT2d& self) { + const auto ptr = self.spectrum(); + const auto dims = self.fft_sizes(); + // property implies reference_internal return value policy + return py::array_t>(dims, ptr); + }) + .def("transform", [](NFFT2d& self, const py::array_t>& z) { + const dims_t shape { + static_cast(z.shape(0)), + static_cast(z.shape(1))}; + const auto itemsize = sizeof(std::complex); + const dims_t strides { + static_cast(z.strides(0) / itemsize), + static_cast(z.strides(1) / itemsize)}; + return self.transform(shape, strides, z.data()); + }) + .def_property_readonly("sizes", &NFFT2d::sizes) + .def_property_readonly("fft_sizes", &NFFT2d::fft_sizes) + ; +} + +template +void addbinding(py::class_>& pyNFFT2dResult) +{ + using dims_t = typename NFFT2d::dims_t; + pyNFFT2dResult + .def("interp", &NFFT2dResult::interp, + py::arg("t"), py::arg("periodic") = true) + ; + // TODO more methods +} + +// instantiate +template void addbinding(py::class_>&); +template void addbinding(py::class_>&); +template void addbinding(py::class_>&); +template void addbinding(py::class_>&); + +void addbinding_make_image_nfft2d(pybind11::module& m) +{ + // TODO generalize to CF32 and CF64 + using T = float; + using array_t = isce3::core::EArray2D>; + m.def("make_image_nfft2d", [](Eigen::Ref image, + py::dict params, + bool pad_input) + { + const auto params_ = parse_nfft2d_params(params); + return makeImageNFFT2d(image, params_, pad_input); + }, + py::arg("image"), + py::arg("params") = py::dict{}, + py::arg("pad_input") = false) + ; +} \ No newline at end of file diff --git a/python/extensions/pybind_isce3/signal/NFFT2d.h b/python/extensions/pybind_isce3/signal/NFFT2d.h new file mode 100644 index 000000000..f3d5aeb52 --- /dev/null +++ b/python/extensions/pybind_isce3/signal/NFFT2d.h @@ -0,0 +1,14 @@ +#pragma once + +#include +#include + +isce3::signal::NFFT2dParams parse_nfft2d_params(const pybind11::dict& params); + +template +void addbinding(pybind11::class_>& pyNFFT2d); + +template +void addbinding(pybind11::class_>& pyNFFT2dResult); + +void addbinding_make_image_nfft2d(pybind11::module& m); \ No newline at end of file diff --git a/python/extensions/pybind_isce3/signal/signal.cpp b/python/extensions/pybind_isce3/signal/signal.cpp index ee62ff26e..d8407253c 100644 --- a/python/extensions/pybind_isce3/signal/signal.cpp +++ b/python/extensions/pybind_isce3/signal/signal.cpp @@ -6,6 +6,7 @@ #include "flatten.h" #include "filter2D.h" #include "multilook.h" +#include "NFFT2d.h" namespace py = pybind11; @@ -19,10 +20,18 @@ void addsubmodule_signal(py::module & m) py::class_ pyCrossmul(m_signal, "Crossmul"); py::class_ pyCrossMultiply(m_signal, "CrossMultiply"); + py::class_> pyNFFT2dF32(m_signal, "NFFT2dF32"); + py::class_> pyNFFT2dF64(m_signal, "NFFT2dF64"); + py::class_> pyNFFT2dF32Result(m_signal, "NFFT2dF32Result"); + py::class_> pyNFFT2dF64Result(m_signal, "NFFT2dF64Result"); // add bindings addbinding(pyCrossmul); addbinding(pyCrossMultiply); + addbinding(pyNFFT2dF32); + addbinding(pyNFFT2dF64); + addbinding(pyNFFT2dF32Result); + addbinding(pyNFFT2dF64Result); addbinding_flatten(m_signal); addbinding_filter2D(m_signal); addbinding_convolve2D(m_signal); @@ -33,4 +42,5 @@ void addsubmodule_signal(py::module & m) addbinding_multilook>>(m_signal); addbinding_multilook>(m_signal); addbinding_multilook>>(m_signal); + addbinding_make_image_nfft2d(m_signal); } diff --git a/python/packages/isce3/focus/__init__.py b/python/packages/isce3/focus/__init__.py index d63c0c2f2..1b893b634 100644 --- a/python/packages/isce3/focus/__init__.py +++ b/python/packages/isce3/focus/__init__.py @@ -1,4 +1,6 @@ from isce3.ext.isce3.focus import * +from . import azcomp_bp +from .serialization import save_polar_grid_to_h5, save_polar_image_to_h5 from .caltone import ToneRemover from .sar_duration import (get_sar_duration, get_radar_velocities, predict_azimuth_envelope) diff --git a/python/packages/isce3/focus/azcomp_bp.py b/python/packages/isce3/focus/azcomp_bp.py new file mode 100644 index 000000000..7575f51a1 --- /dev/null +++ b/python/packages/isce3/focus/azcomp_bp.py @@ -0,0 +1,447 @@ +#!/usr/bin/env python3 +""" +Azimuth compression algorithms for SAR focusing. + +Provides standard and factorized backprojection implementations. +""" +from __future__ import annotations +from collections import defaultdict +from functools import lru_cache +import logging +import numpy as np +import h5py +import isce3 +from isce3.core import LUT2d +from isce3.focus.serialization import BackprojectionStageParameters +from isce3.geometry import DEMInterpolator +from isce3.product import RadarGridParameters +from typing import Optional + +log = logging.getLogger("isce3.focus.azcomp_bp") + +# Type aliases for processing block plan +Selection2d = tuple[slice, slice] +TimeBounds = tuple[float, float] +BlockPlan = list[tuple[Selection2d, TimeBounds]] + + +def is_overlapping(a, b, c, d): + """ + Check if two intervals overlap. + + Parameters + ---------- + a : float + Start of first interval. + b : float + End of first interval. + c : float + Start of second interval. + d : float + End of second interval. + + Returns + ------- + bool + True if intervals [a, b] and [c, d] overlap, False otherwise. + """ + assert (b >= a) and (d >= c) + return (d >= a) and (c <= b) + + +class Task: + """ + Deferred task wrapper for lazy evaluation. + + Stores a function and its arguments for later execution, allowing + tasks to be defined without immediate evaluation. Used in factorized + backprojection to reduce memory pressure by only computing intermediate + results when needed. + + Parameters + ---------- + function : callable + Function to execute when result() is called. + *args + Positional arguments to pass to function. + **kwargs + Keyword arguments to pass to function. + """ + def __init__(self, function, *args, **kwargs): + self.function = function + self.args = args + self.kwargs = kwargs + + def result(self): + """ + Execute the stored function with its arguments. + + Returns + ------- + Any + Result of calling function(*args, **kwargs). + """ + return self.function(*self.args, **self.kwargs) + + +def find_min_cache_size(key_lists): + """ + Determine the minimum LRU cache size needed to avoid cache misses. + + Parameters + ---------- + key_lists: Iterable[Iterable[Hashable]] + List of jobs, where each job is a list of task keys, and + keys may be shared between jobs. + + Returns + ------- + n : int + Minimum cache size required to hold shared keys in memory, + assuming jobs are executed in the given order. + """ + # Flatten the key sequence into a single ordered list of key accesses + access_sequence = [key for keylist in key_lists for key in keylist] + + # For each key, record the indices where it's accessed + access_indices = defaultdict(list) + for i, key in enumerate(access_sequence): + access_indices[key].append(i) + + min_size = 1 + + for key, indices in access_indices.items(): + # Only care about keys accessed more than once (re-use case) + for j in range(1, len(indices)): + prev_idx = indices[j - 1] + curr_idx = indices[j] + + # Count distinct keys in the window [prev_idx, curr_idx] inclusive. + # If this many distinct keys were accessed, the LRU cache must hold + # at least this many entries to avoid evicting `key` before reuse. + window = access_sequence[prev_idx:curr_idx + 1] + distinct_in_window = len(set(window)) + min_size = max(min_size, distinct_in_window) + + return min_size + + +def nfft_params_dict(p: isce3.focus.serialization.NonUniformFFT2DParams): + return dict( + rows = dict( + m = p.azimuth.kernel_halfwidth, + s = p.azimuth.zero_padding_factor), + cols = dict( + m = p.range.kernel_halfwidth, + s = p.range.zero_padding_factor)) + + +def azcomp_bp(azres, kernel, blocks_bounds, igeom, rcdata, ogrid, writer, + height=None, dem=isce3.geometry.DEMInterpolator(), + rdr2geo_params=dict(), geo2rdr_params=dict(), atmos="nodelay", + use_gpu=False): + """ + Perform azimuth compression using standard backprojection algorithm. + + Parameters + ---------- + azres : float + Desired azimuth resolution, in meters. + kernel : isce3.core.Kernel + Interpolation kernel for backprojection. + blocks_bounds : BlockPlan + List of tuples containing ((row_slice, col_slice), (t0, t1)) for each + processing block, where slices define the output grid region and (t0, t1) + are the required raw data time bounds in seconds. + igeom : isce3.container.RadarGeometry + Input radar geometry for range-compressed data. + rcdata : array-like + Range-compressed data, shape (azimuth, range). + ogrid : RadarGridParameters + Output zero-Doppler radar grid parameters. + writer : BackgroundWriter + Writer object. Must have a method `queue_write(z, block)` for writing + out image subset `z` into selection `block` of output image where `z` + is a 2D numpy array and `block` is a tuple[slice, slice]. + height : array-like, optional + Optional storage for height above ellipsoid (in meters) for each output + pixel, shape matching ogrid. + dem : isce3.geometry.DEMInterpolator, optional + Digital elevation model. Default is ellipsoid (height=0). + rdr2geo_params : dict, optional + Parameters for rdr2geo_bracket solver. + geo2rdr_params : dict, optional + Parameters for geo2rdr_bracket solver. + atmos : str, optional + Atmospheric delay model. Default is "nodelay". + use_gpu : bool, optional + Use GPU acceleration if available. Default is False. + """ + if use_gpu: + backproject = isce3.cuda.focus.backproject + else: + backproject = isce3.focus.backproject + fc = isce3.core.speed_of_light / ogrid.wavelength + zerodop = isce3.core.LUT2d() + for block, (t0, t1) in blocks_bounds: + description = f"(i, j) = ({block[0].start}, {block[1].start})" + if not is_overlapping(t0, t1, igeom.radar_grid.sensing_start, + igeom.radar_grid.sensing_stop): + log.info(f"Skipping inactive azcomp block at {description}") + continue + log.info(f"Azcomp block at {description}") + bgrid = ogrid[block] + ogeom = isce3.container.RadarGeometry(bgrid, igeom.orbit, zerodop) + z = np.zeros(bgrid.shape, 'c8') + hgt = height[block] if height is not None else None + err = backproject(z, ogeom, rcdata, igeom, dem, fc, azres, kernel, + atmos, rdr2geo_params, geo2rdr_params, height=hgt) + if err: + log.warning("azcomp block contains some invalid pixels") + writer.queue_write(z, block) + + +def azcomp_fbp(factors: BackprojectionStageParameters, + azres, kernel, blocks_bounds, igeom, + rcdata, ogrid, writer, height=None, dem=isce3.geometry.DEMInterpolator(), + rdr2geo_params=dict(), geo2rdr_params=dict(), atmos="nodelay", + use_gpu=False, bandwidth=0.0, debugfile=None): + """ + Perform azimuth compression using factorized backprojection algorithm. + + The factorized backprojection algorithm processes data in multiple stages, + first focusing blocks of pulses to intermediate polar grids, then merging + those grids hierarchically, and finally accumulating the results to the + output zero-Doppler radar grid. + + Parameters + ---------- + factors : BackprojectionStageParameters + List of factorization stage parameters defining the processing hierarchy. + Each stage specifies the number of pulses or polar images to combine, + oversample factors, and NFFT interpolation parameters. + azres : float + Desired azimuth resolution, in meters. + kernel : isce3.core.Kernel + Interpolation kernel for backprojection. + blocks_bounds : BlockPlan + List of tuples containing ((row_slice, col_slice), (t0, t1)) for each + processing block, where slices define the output grid region and (t0, t1) + are the required raw data time bounds in seconds. + igeom : isce3.container.RadarGeometry + Input radar geometry for range-compressed data. + rcdata : array-like + Range-compressed data, shape (azimuth, range). + ogrid : RadarGridParameters + Output zero-Doppler radar grid parameters. + writer : BackgroundWriter + Writer object. Must have a method `queue_write(z, block)` for writing + out image subset `z` into selection `block` of output image where `z` + is a 2D numpy array and `block` is a tuple[slice, slice]. + height : array-like, optional + Height above ellipsoid (in meters) for each output pixel, shape matching + ogrid. If None, uses DEM. + dem : isce3.geometry.DEMInterpolator, optional + Digital elevation model. Default is ellipsoid (height=0). + rdr2geo_params : dict, optional + Parameters for rdr2geo_bracket solver. + geo2rdr_params : dict, optional + Parameters for geo2rdr_bracket solver. + atmos : str, optional + Atmospheric delay model. Applied at final stage only to avoid phase + modulation artifacts from DEM sampling across subimages. Default is "nodelay". + use_gpu : bool, optional + Use GPU acceleration if available. Default is False. + bandwidth : float, optional + Signal bandwidth in Hz, used to determine polar grid angular extent. + Default is 0.0. + debugfile : file-like, optional + HDF5 file handle for writing intermediate polar grids for debugging. + If None, no debug output is written. + """ + fc = isce3.core.speed_of_light / ogrid.wavelength + zerodop = isce3.core.LUT2d() + + if use_gpu: + bp_to_polar_grid = isce3.cuda.focus.backproject_to_polar_grid + merge_polar_images = isce3.cuda.focus.merge_polar_images + add_to_radar_grid = isce3.cuda.focus.accumulate_polar_images_to_radar_grid + # TODO could make this a separate option to conserve GPU memory. + # That'd require a little finess in the bindings and CUDA side, though. + make_image_nfft2d = isce3.cuda.signal.make_image_nfft2d + else: + bp_to_polar_grid = isce3.focus.backproject_to_polar_grid + merge_polar_images = isce3.focus.merge_polar_images + add_to_radar_grid = isce3.focus.accumulate_polar_images_to_radar_grid + make_image_nfft2d = isce3.signal.make_image_nfft2d + + _, v = igeom.orbit.interpolate(igeom.orbit.mid_time) + vs = np.linalg.norm(v) + tq_max = isce3.focus.get_polar_angle_time_constant(fc, vs, bandwidth) + + if debugfile is not None: + log.debug(f"Writing FBP metadata to file {debugfile.name}") + with h5py.File(debugfile, "w") as h5: + epoch = igeom.reference_epoch + igeom.orbit.save_to_h5(h5.require_group("orbit")) + igeom.doppler.save_to_h5(h5.require_group("doppler"), "doppler", + epoch, "Hz") + h5.create_dataset("epoch", data=np.bytes_(epoch)) + h5.create_dataset("wavelength", data=igeom.radar_grid.wavelength) + + # Focus to intermediate grids. + # NOTE We'll actually just define the tasks and only evaluate them as needed + # in order to reduce memory pressure. We could process them in + # parallel using concurrent.futures or dask, for for now just store them in + # a dict keyed by the PolarGrid associated with the imagelets. + tasks = dict() + aztimes = np.array(igeom.radar_grid.sensing_times) + pris = np.hstack((np.diff(aztimes), aztimes[-1] - aztimes[-2])) + stage = factors[0] + nfft2d_params = nfft_params_dict(stage.interpolation) + pulse_starts = range(0, igeom.radar_grid.length, stage.size) + log.info(f"Beginning initial factorizations of {stage.size} pulses") + nblocks = len(pulse_starts) + + def process_pulses(iblock, nblocks, debugfile, fdata, sr, x, v, + polar_grid, dem, fc, kernel, atmos, rdr2geo_params): + log.info(f"Focusing {len(x)} pulses to polar image {iblock + 1} of {nblocks}") + # NOTE Atmosphere will get applied (if requested) at final stage to + # avoid phase modulation from DEM sampling issues across subimages. + # Always "nodelay" in this stage. + _, img, _ = bp_to_polar_grid(fdata, sr, x, v, + polar_grid, dem, fc, kernel, "nodelay", rdr2geo_params) + if debugfile is not None: + import h5py + log.debug(f"Dumping FBP factor with shape = {img.shape} to file.") + with h5py.File(debugfile, "w") as h5: # okay to reopen stream + g = h5.require_group(f"stage_00/block_{iblock:06d}") + isce3.focus.save_polar_image_to_h5(img, polar_grid, g) + log.info("NFFT upsampling and filtering") + return make_image_nfft2d(img, nfft2d_params, pad_input=True) + + for i in pulse_starts: + pulses = slice(i, i + stage.size) + ti = aztimes[pulses] + if len(ti) < 2: + log.info("Skipping FBP block containing only a single pulse.") + continue + fgrid = igeom.radar_grid[pulses, :] + fgeom = isce3.container.RadarGeometry(fgrid, igeom.orbit, igeom.doppler) + fdata = rcdata[pulses, :] + iblock = i // stage.size + polar_grid, x, v = isce3.focus.setup_polar_grid_for_pulses(fgeom, ti, + bandwidth, azres, stage.oversample_range, stage.oversample_azimuth, + pri=pris[pulses][-1]) + tasks[polar_grid] = Task(process_pulses, iblock, nblocks, + debugfile, fdata, fgrid.slant_ranges, x, v, + polar_grid, dem, fc, kernel, atmos, rdr2geo_params) + + grids = sorted(tasks.keys(), key = lambda grid: grid.aztime_start) + + # Merge polar grids to make bigger polar grids. + # With Python 3.12 we could use itertools.batched + def process_merge(i_stage, i_block, nblocks, in_grids, out_grid, nfft2d_params): + # Process the input data we need. Middle stages don't overlap, so no + # harm in popping the task off the stack. + in_images = [tasks.pop(grid).result() for grid in in_grids] + log.info(f"Merging {len(in_grids)} polar images stage {i_stage} block " + f"{i_block + 1} / {nblocks}") + out_image = np.zeros(out_grid.shape, np.complex64) + merge_polar_images(in_grids, in_images, out_grid, out_image, + fc, dem, rdr2geo_params) + if debugfile is not None: + import h5py + name = f"stage_{i_stage + 1:02d}/block_{i_block:06d}" + with h5py.File(debugfile, "w") as h5: # okay to reopen stream + g = h5.require_group(name) + isce3.focus.save_polar_image_to_h5(out_image, out_grid, g) + log.info("NFFT upsampling and filtering") + return make_image_nfft2d(out_image, nfft2d_params, pad_input=True) + + num_middle_stages = len(factors[1:]) + for i_stage, stage in enumerate(factors[1:]): + log.info("Planning intermediate factorization stage " + + f"{i_stage + 1} / {num_middle_stages}") + + # Don't let azimuth resolution grow finer than user requested one. + dq_min = azres / (ogrid.slant_ranges[-1] * stage.oversample_azimuth) + tq = tq_max / stage.oversample_azimuth + + nfft2d_params = nfft_params_dict(stage.interpolation) + input_block_starts = range(0, len(grids), stage.size) + nblocks = len(input_block_starts) + stage_grids = [] + + for i in input_block_starts: + i_block = i // stage.size + mask = slice(i, i + stage.size) + input_grids = grids[mask] + if len(input_grids) == 1: + # No need to merge. Grid is already tasked, though be sure to + # carry it forward to next stage. + stage_grids.append(input_grids[0]) + continue + my_grid = isce3.focus.merge_polar_grids(input_grids, dem, + rdr2geo_params, dq_min, tq) + tasks[my_grid] = Task(process_merge, i_stage, i_block, nblocks, + input_grids, my_grid, nfft2d_params) + stage_grids.append(my_grid) + + # Use this stage's grids as input for next stage. + grids = stage_grids + + # Plan final stage to get bound on LRU cache size, assuming FIFO access + # pattern. + blocks_grids = list() + for block, (t0, t1) in blocks_bounds: + description = f"(i, j) = ({block[0].start}, {block[1].start})" + if not is_overlapping(t0, t1, igeom.radar_grid.sensing_start, + igeom.radar_grid.sensing_stop): + log.info(f"Will skip inactive azcomp block at {description}") + continue + active_grids = [grid for grid in grids + if is_overlapping(t0, t1, grid.aztime_start, grid.aztime_end)] + blocks_grids.append((block, active_grids)) + + max_images = max(len(grids) for (_, grids) in blocks_grids) + log.info(f"Proceeding to final stage with max {max_images} sub-images per block") + # Required cache size may be smaller than max_images when not all + # sub-images in one block are used in the next block. However, it can + # also be more when we subdivide in range, since a far-range block may + # need all the sub-images of a near-range block plus a few more. + cache_size = find_min_cache_size([grid for (_, grid) in blocks_grids]) + log.info(f"Calculated min cache size = {cache_size}") + + @lru_cache(maxsize=cache_size) + def get_image_iterpolator(polar_grid): + # Using pop() to remove from stack requires that cache size is adequate + # to avoid redundant computations, which we prioritize over generality. + try: + return tasks.pop(polar_grid).result() + except KeyError as err: + msg = ("Failed to retrieve sub-image spanning time interval " + f"[{polar_grid.aztime_start}, {polar_grid.aztime_end}). This " + "could mean that the stripmap assumption was violated.") + log.error(msg) + raise + + # sum factors into final image + for block, active_grids in blocks_grids: + description = f"(i, j) = ({block[0].start}, {block[1].start})" + active_images = [get_image_iterpolator(grid) for grid in active_grids] + bgrid = ogrid[block] + ogeom = isce3.container.RadarGeometry(bgrid, igeom.orbit, zerodop) + z = np.zeros(bgrid.shape, 'c8') + hgt = height[block] if height is not None else None + log.info(f"Azcomp final sums for block at {description} using " + f"{len(active_images)} sub-apertures") + err = add_to_radar_grid( + z, ogeom, igeom.orbit, igeom.doppler, active_grids, active_images, + dem, fc, azres, atmos, rdr2geo_params, geo2rdr_params, hgt) + if err: + log.warning("azcomp block contains some invalid pixels") + writer.queue_write(z, block) + + if len(tasks) != 0: + log.warning(f"Queued {len(tasks)} tasks that were never needed.") diff --git a/python/packages/isce3/focus/serialization.py b/python/packages/isce3/focus/serialization.py new file mode 100644 index 000000000..b397efbe9 --- /dev/null +++ b/python/packages/isce3/focus/serialization.py @@ -0,0 +1,252 @@ +from dataclasses import dataclass +import h5py +from isce3.core import Linspace, LookSide +from isce3.focus import PolarGrid +import numpy as np + +# FIXME Not sure where to put this stuff. Maybe monkey patch the classes? + +def overwrite(group: h5py.Group, key: str, value): + """Overwrite or create a dataset in an HDF5 group. + + Parameters + ---------- + group : h5py.Group + HDF5 group to modify. + key : str + Name of the dataset. + value + Data to store in the dataset. + """ + if key in group: + del group[key] + group.create_dataset(key, data=value) + + +def save_linspace_to_h5(x: Linspace, group: h5py.Group): + """Save a Linspace object to an HDF5 group. + + Parameters + ---------- + x : Linspace + Linspace object to save. + group : h5py.Group + HDF5 group where the linspace will be stored. + """ + for key in ("first", "spacing", "size"): + val = getattr(x, key) + overwrite(group, key, val) + + +def load_linspace_from_h5(group: h5py.Group) -> Linspace: + """Load a Linspace object from an HDF5 group. + + Parameters + ---------- + group : h5py.Group + HDF5 group containing the linspace data. + + Returns + ------- + Linspace + Reconstructed Linspace object. + """ + args = [group[key][()] for key in ("first", "spacing", "size")] + return Linspace(*args) + + +def save_lookside_to_h5(side: LookSide, group: h5py.Group): + """Save a LookSide enum to an HDF5 group. + + Parameters + ---------- + side : LookSide + LookSide enum value to save. + group : h5py.Group + HDF5 group where the look side will be stored. + """ + side_str = str(side).split(".")[1] + overwrite(group, "look_side", np.bytes_(side_str)) + + +def load_lookside_from_h5(group: h5py.Group) -> LookSide: + """Load a LookSide enum from an HDF5 group. + + Parameters + ---------- + group : h5py.Group + HDF5 group containing the look side data. + + Returns + ------- + LookSide + LookSide enum value. + """ + side_str_lower = group["look_side"][()].decode("utf-8").lower() + valid_sides = {"left": LookSide.Left, "right": LookSide.Right} + return valid_sides[side_str_lower] + + +def save_polar_grid_to_h5(grid: PolarGrid, group: h5py.Group): + """Save a PolarGrid object to an HDF5 group. + + Parameters + ---------- + grid : PolarGrid + PolarGrid object to save. + group : h5py.Group + HDF5 group where the polar grid will be stored. + """ + for key in ("aztime_start", "aztime_end", "origin", "axis"): + val = getattr(grid, key) + overwrite(group, key, val) + for key in ("range", "sin_squint"): + g = group.require_group(key) + val = getattr(grid, key) + save_linspace_to_h5(val, g) + save_lookside_to_h5(grid.look_side, group) + + +def load_polar_grid_from_h5(group: h5py.Group) -> PolarGrid: + """Load a PolarGrid object from an HDF5 group. + + Parameters + ---------- + group : h5py.Group + HDF5 group containing the polar grid data. + + Returns + ------- + PolarGrid + Reconstructed PolarGrid object. + """ + keys = ("aztime_start", "aztime_end", "origin", "axis") + args = [group[key][()] for key in keys] + args.append(load_linspace_from_h5(group["range"])) + args.append(load_linspace_from_h5(group["sin_squint"])) + args.append(load_lookside_from_h5(group)) + return PolarGrid(*args) + + +def save_polar_image_to_h5(z: np.ndarray, grid: PolarGrid, group: h5py.Group): + """Save a polar image and its grid to an HDF5 group. + + Parameters + ---------- + z : np.ndarray + Image data array. + grid : PolarGrid + PolarGrid defining the image geometry. + group : h5py.Group + HDF5 group where the image and grid will be stored. + """ + overwrite(group, "image", z) + g = group.require_group("polar_grid") + save_polar_grid_to_h5(grid, g) + + +def load_polar_image_from_h5(group: h5py.Group) -> tuple[np.ndarray, PolarGrid]: + """Load a polar image and its grid from an HDF5 group. + + Parameters + ---------- + group : h5py.Group + HDF5 group containing the image and polar grid data. + + Returns + ------- + z : np.ndarray + Image array. + grid : PolarGrid + PolarGrid object. + """ + z = group["image"][:] + grid = load_polar_grid_from_h5(group["polar_grid"]) + return (z, grid) + +@dataclass(frozen=True) +class NonUniformFFTParameters: + """Parameters for non-uniform FFT interpolation. + + Attributes + ---------- + zero_padding_factor : float + Oversampling factor for FFT grid (must be > 1.0). + kernel_halfwidth : int + Half-width of interpolation kernel (must be >= 1). + """ + zero_padding_factor: float = 2.0 + kernel_halfwidth: int = 2 + + def __post_init__(self): + if self.zero_padding_factor <= 1.0: + raise ValueError("require NFFT zero_padding_factor > 1.0") + if self.kernel_halfwidth < 1: + raise ValueError("require NFFT kernel_halfwidth >= 1") + + @classmethod + def from_dict(cls, d: dict): + default = cls() + return cls( + float(d.get("zero_padding_factor", default.zero_padding_factor)), + int(d.get("kernel_halfwidth", default.kernel_halfwidth))) + +@dataclass(frozen=True) +class NonUniformFFT2DParameters: + """Parameters for 2D non-uniform FFT interpolation. + + Attributes + ---------- + range : NonUniformFFTParameters + NFFT parameters for range dimension. + azimuth : NonUniformFFTParameters + NFFT parameters for azimuth dimension. + """ + range: NonUniformFFTParameters = NonUniformFFTParameters() + azimuth: NonUniformFFTParameters = NonUniformFFTParameters() + + @classmethod + def from_dict(cls, d: dict): + default = cls() + T = NonUniformFFTParameters + rg = T.from_dict(d["range"]) if "range" in d else default.range + az = T.from_dict(d["azimuth"]) if "azimuth" in d else default.azimuth + return cls(rg, az) + +@dataclass(frozen=True) +class BackprojectionStageParameters: + """Parameters for a factorized backprojection processing stage. + + Attributes + ---------- + size : int + Number of pulses/subapertures per stage (must be >= 1). + oversample_range : float + Range oversampling factor (must be >= 1.0). + oversample_azimuth : float + Azimuth oversampling factor (must be >= 1.0). + interpolation : NonUniformFFT2DParameters + NFFT parameters for 2D interpolation. + """ + size: int = 1 + oversample_range: float = 1.2 + oversample_azimuth: float = 1.2 + interpolation: NonUniformFFT2DParameters = NonUniformFFT2DParameters() + + def __post_init__(self): + if self.size < 1: + raise ValueError("require at least 1 pulse/subaperture per stage") + if self.oversample_range < 1.0: + raise ValueError("must sample range at or above Nyquist limit") + if self.oversample_azimuth < 1.0: + raise ValueError("must sample azimuth at or above Nyquist limit") + + @classmethod + def from_dict(cls, d: dict): + default = cls() + key, T = "interpolation", NonUniformFFT2DParameters + return cls( + d.get("size", default.size), + d.get("oversample_range", default.oversample_range), + d.get("oversample_azimuth", default.oversample_azimuth), + T.from_dict(d[key]) if key in d else default.interpolation) \ No newline at end of file diff --git a/python/packages/nisar/workflows/focus.py b/python/packages/nisar/workflows/focus.py index 596b4dbb0..8702a698a 100644 --- a/python/packages/nisar/workflows/focus.py +++ b/python/packages/nisar/workflows/focus.py @@ -35,7 +35,9 @@ from isce3.core import DateTime, TimeDelta, LUT2d, Attitude, Orbit from isce3.focus import (make_los_luts, fill_gaps, make_cal_luts, Notch, find_bad_rangline_slices) -from isce3.geometry import los2doppler +from isce3.focus.azcomp_bp import (azcomp_bp, azcomp_fbp, BlockPlan, + TimeBounds) +from isce3.focus.serialization import BackprojectionStageParameters from isce3.io.gdal import Raster, GDT_CFloat32 from isce3.product import (RadarGridParameters, get_radar_grid_nominal_ground_spacing) @@ -739,10 +741,6 @@ def get_geo2rdr_params(cfg: Struct, orbit: Optional[Orbit] = None) -> dict: return geo2rdr_params -Selection2d = tuple[slice, slice] -TimeBounds = tuple[float, float] -BlockPlan = list[tuple[Selection2d, TimeBounds]] - def plan_processing_blocks(cfg: Struct, grid: RadarGridParameters, doppler: LUT2d, dem: isce3.geometry.DEMInterpolator, orbit: Orbit, pad: float = 0.1) -> BlockPlan: @@ -809,10 +807,6 @@ def total_bounds(blocks_bounds: BlockPlan) -> TimeBounds: return (begin, end) -def is_overlapping(a, b, c, d): - assert (b >= a) and (d >= c) - return (d >= a) and (c <= b) - def get_kernel(cfg: Struct): # TODO opt = cfg.processing.azcomp.kernel @@ -1706,6 +1700,15 @@ def get_focused_sub_swaths(rawlist, out_chan, grid, orbit, doppler, dem, azres, return swaths +def get_azcomp_stage_config(cfg: Struct): + factors = cfg.processing.azcomp.factorization + if not isinstance(factors, Iterable) or len(factors) < 1: + raise ValueError("Must specify at least one factorization stage " + "in config file.") + T = isce3.focus.serialization.BackprojectionStageParameters + return [T.from_dict(struct2dict(factor)) for factor in factors] + + def get_caltone_algorithm(cfg, fc, fs, n, is_dithered): """Helper for configuring caltone removal. @@ -1788,10 +1791,6 @@ def focus(runconfig, runconfig_path=""): log.info(f"Processing using CUDA device {device.id} ({device.name})") - backproject = isce3.cuda.focus.backproject - else: - backproject = isce3.focus.backproject - # Generate output grids. grid_epoch, t0, t1, r0, r1 = get_total_grid_bounds(rawnames) log.info(f"Raw data time spans [{t0}, {t1}] seconds since {grid_epoch}.") @@ -2323,27 +2322,24 @@ def temp(suffix): # Do azimuth compression. igeom = isce3.container.RadarGeometry(rc_grid, orbit, dop[frequency]) - - for block, (t0, t1) in blocks_bounds[frequency]: - description = f"(i, j) = ({block[0].start}, {block[1].start})" - if not cfg.processing.is_enabled.azcomp: - continue - if not is_overlapping(t0, t1, - rc_grid.sensing_start, rc_grid.sensing_stop): - log.info(f"Skipping inactive azcomp block at {description}") - continue - log.info(f"Azcomp block at {description}") - bgrid = ogrid[frequency][block] - ogeom = isce3.container.RadarGeometry(bgrid, orbit, zerodop) - z = np.zeros(bgrid.shape, 'c8') - hgt = hgt_mm[block] if dump_height else None - err = backproject(z, ogeom, rcfile.data, igeom, dem, - channel_out.band.center, azres, - kernel, atmos, get_rdr2geo_params(cfg), - get_geo2rdr_params(cfg, orbit), height=hgt) - if err: - log.warning("azcomp block contains some invalid pixels") - writer.queue_write(z, block) + if cfg.processing.is_enabled.azcomp: + factors = get_azcomp_stage_config(cfg) + if factors[0].size > 1: + debugfile = (temp(f"_{frequency}{pol}_fbp_factors.h5") + if not cfg.processing.delete_tempfiles else None) + azcomp_fbp(factors, azres, kernel, + blocks_bounds[frequency], igeom, + rcfile.data, ogrid[frequency], writer, + hgt_mm if dump_height else None, dem, + get_rdr2geo_params(cfg), get_geo2rdr_params(cfg, orbit), + atmos, use_gpu, channel_out.band.width, + debugfile) + else: + azcomp_bp(azres, kernel, blocks_bounds[frequency], igeom, + rcfile.data, ogrid[frequency], writer, + hgt_mm if dump_height else None, dem, + get_rdr2geo_params(cfg), get_geo2rdr_params(cfg, orbit), + atmos, use_gpu) # Raster/GDAL creates a .hdr file we have to clean up manually. hdr = fd.name.replace(".c8", ".hdr") @@ -2385,7 +2381,7 @@ def configure_logging(): sh.setFormatter(fmt) log.addHandler(sh) for friend in ("Raw", "SLCWriter", "nisar.antenna.pattern", "rslc_cal", - "isce3.focus.notch"): + "isce3.focus.notch", "isce3.focus.azcomp_bp"): l = logging.getLogger(friend) l.setLevel(log_level) l.addHandler(sh) diff --git a/share/isce3/factorized_bp_gif.py b/share/isce3/factorized_bp_gif.py new file mode 100755 index 000000000..f99b4e91d --- /dev/null +++ b/share/isce3/factorized_bp_gif.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +import argparse +import h5py +import numpy as np +from PIL import Image +from tqdm import tqdm + +parser = argparse.ArgumentParser() +parser.add_argument("factors", help="HDF5 file containing FBP sub-images") +parser.add_argument("-r", "--looks-range", type=int, default=0) +parser.add_argument("-a", "--looks-azimuth", type=int, default=0) +parser.add_argument("-o", help="output animation", default="factors.gif") +parser.add_argument("--stage", type=int, default=0) +parser.add_argument("--cw", type=float, default=0.5) +parser.add_argument("--exp", type=float, default=1.0) +parser.add_argument("--duration", type=int, default=10) +args = parser.parse_args() + +def multilook(z, ny=1, nx=1, f=lambda z: z): + m, n = z.shape + mout, nout = m // ny, n // nx + x = f(z[:(mout * ny), :(nout * nx)]) + x.shape = mout, ny, nout, nx + return x.mean(axis=(1, 3)) + +def powlooks(z, ny=1, nx=1): + return multilook(z, ny, nx, f = lambda z: z.real**2 + z.imag**2) + +h5 = h5py.File(args.factors, mode="r") +group_name = f"stage_{args.stage:02d}" +group = h5[group_name] +blocks = sorted([key for key in group if key.startswith("block_")]) + +nr, na = args.looks_range, args.looks_azimuth +if args.looks_range == 0 or args.looks_azimuth == 0: + grid = group[blocks[0]]["polar_grid"] + origin = grid["origin"][:] + # law of cosines + o = np.linalg.norm(origin) + r = grid["range/first"][()] + a = 6378137. + look = np.arccos((o**2 + r**2 - a**2) / (2 * o * r)) + print("look angle (deg) =", np.rad2deg(look)) + dg = grid["range/spacing"][()] / np.sin(look) + ds = grid["sin_squint/spacing"] * r + if ds > dg: + na = 1 + nr = round(ds / dg) + else: + nr = 1 + na = round(dg / ds) + print("looks range =", nr) + print("looks azimuth =", na) + +# Figure out max image size to use for all frames. +image_shapes = [] +for block_name in blocks: + image_shapes.append(group[block_name]["image"].shape) +rows, cols = np.max(np.array(image_shapes), axis=0) +buf = np.zeros((rows, cols), dtype="c8") + +power_images = [] +for block_name in tqdm(blocks, "reading"): + buf[:] = np.nan + slc = group[block_name]["image"][:] + # Assume we should center each image in the frame. + i0 = (rows - slc.shape[0]) // 2 + i1 = i0 + slc.shape[0] + j0 = (cols - slc.shape[1]) // 2 + j1 = j0 + slc.shape[1] + buf[i0:i1, j0:j1] = slc + power_images.append(powlooks(buf, na, nr)) + +mean = np.nanmean(power_images) +cw_scale = 1.0 +if mean > 0: + cw_scale = 0.7 * args.cw / mean + +bitmaps = [] +for zpp in tqdm(power_images, "scaling"): + x = cw_scale * zpp + x **= args.exp + np.clip(x, 0, 1, x) + arr = (255 * np.nan_to_num(x)).astype(np.uint8) + img = Image.fromarray(arr, mode="L") + bitmaps.append(img) + +bitmaps[0].save(args.o, save_all=True, append_images=bitmaps[1:], + optimize=False, duration=args.duration, loop=0) diff --git a/share/nisar/defaults/focus.yaml b/share/nisar/defaults/focus.yaml index 8bbac52dc..cd4c909a6 100644 --- a/share/nisar/defaults/focus.yaml +++ b/share/nisar/defaults/focus.yaml @@ -429,7 +429,7 @@ runconfig: # Azimuth compression can be tiled arbitrarily, though # dimensions will affect runtime. block_size: - range: 32768 + range: 65536 azimuth: 1024 # Desired azimuth resolution in meters. @@ -446,6 +446,11 @@ runconfig: fit: Table # null or Cheby or Table fit_order: 2048 + # List of aperture factorization parameters. + # Default to two-stage with 128 pulses per subaperture. + factorization: + - size: 128 + dry_troposphere_model: tsx dem: @@ -563,7 +568,7 @@ runconfig: compression_type: gzip # Level of compression applied to raster - compression_level: 4 + compression_level: 2 # Chunk size of raster. Enter [-1, -1] to disable chunks. chunk_size: [512, 512] diff --git a/share/nisar/schemas/focus.yaml b/share/nisar/schemas/focus.yaml index 76d5a64e6..6cb27af41 100644 --- a/share/nisar/schemas/focus.yaml +++ b/share/nisar/schemas/focus.yaml @@ -219,6 +219,8 @@ runconfig: fit: enum('Cheby', 'Table', required=False) fit_order: int(min=1, required=False) + factorization: list(include('fbp_parameters'), min=1, required=False) + dry_troposphere_model: enum('nodelay', 'tsx', required=False) @@ -662,3 +664,34 @@ output_options: # Ideally at least large enough to hold `chunk_size` amount of uncompressed # data as well as its metadata. fs_page_size: int(min=512, max=1073741824, required=False) + + +nfft_parameters: + # Minimum spectral zero padding factor (ratio of transform size to data + # size). Actual padding may be larger to achive efficient transform size. + zero_padding_factor: num(min=1.0, required=False) + + # Size of interpolator is equal to (1 + 2 * kernel_halfwidth) + kernel_halfwidth: int(min=1, required=False) + +nfft2d_parameters: + range: include("nfft_parameters", required=False) + azimuth: include("nfft_parameters", required=False) + +fbp_parameters: + # For the initial factorization stage, this is the number of pulses to + # merge into each subaperture. For later stages, this is this number of + # subapertures to merge together. + size: int(min=1, required=False) + + # Minimum ratio of sample rate to Nyquist limit in range in polar grid. + # Actual ratio may be larger to achieve efficient transform sizes. + oversample_range: num(min=1.0, required=False) + + # Minimum ratio of sample rate to Nyquist limit in azimuth in polar grid. + # Actual ratio may be larger to achieve efficient transform sizes. + oversample_azimuth: num(min=1.0, required=False) + + # Settings for interpolating the resulting polar grid using 2D non-uniform + # fast Fourier transform (NFFT). + interpolation: include("nfft2d_parameters", required=False) \ No newline at end of file diff --git a/tests/cxx/isce3/core/linspace/linspace.cpp b/tests/cxx/isce3/core/linspace/linspace.cpp index e2f12699a..16c130f2c 100644 --- a/tests/cxx/isce3/core/linspace/linspace.cpp +++ b/tests/cxx/isce3/core/linspace/linspace.cpp @@ -198,6 +198,20 @@ TEST(LinspaceTest, Comparison) EXPECT_TRUE( x1 != x3 ); } +TEST(LinspaceTest, Bounds) +{ + double first = 1.0; + double spacing = 2.0; + int size = 5; + + isce3::core::Linspace x(first, spacing, size); + + auto [leading, trailing] = x.bounds(); + + EXPECT_DOUBLE_EQ( leading, 0.0 ); + EXPECT_DOUBLE_EQ( trailing, 10.0 ); +} + int main(int argc, char * argv[]) { testing::InitGoogleTest(&argc, argv); diff --git a/tests/cxx/isce3/cuda/Sources.cmake b/tests/cxx/isce3/cuda/Sources.cmake index dbb6f96a6..28e2cb1cb 100644 --- a/tests/cxx/isce3/cuda/Sources.cmake +++ b/tests/cxx/isce3/cuda/Sources.cmake @@ -27,4 +27,5 @@ signal/gpuCrossMul.cpp signal/gpuFilter.cpp signal/gpuLooks.cpp signal/gpuSignal.cpp +signal/nfft2d.cu ) diff --git a/tests/cxx/isce3/cuda/signal/nfft2d.cu b/tests/cxx/isce3/cuda/signal/nfft2d.cu new file mode 100644 index 000000000..07b9a7b20 --- /dev/null +++ b/tests/cxx/isce3/cuda/signal/nfft2d.cu @@ -0,0 +1,56 @@ +#include +#include +#include + +using isce3::cuda::signal::NFFT2d; +using isce3::cuda::signal::NFFT2dResultView; + +__global__ void +interp(const NFFT2dResultView result, thrust::complex* z0) +{ + const auto tid = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + + if (tid > 0) { + return; + } + + std::array t = {0.0, 0.0}; + const auto z = result.interp(t, true); + *z0 = z; +} + +TEST(nfft2d, ctor) +{ + using T = float; + using dims_t = NFFT2d::dims_t; + const dims_t m = {2, 2}; + const dims_t sizes = {201, 80}; + const dims_t fft_sizes = {512, 256}; + auto ft = NFFT2d(m, sizes, fft_sizes); + + // create a spectrum equal to one everywhere + auto npix = sizes[0] * sizes[1]; + auto spectrum = std::vector>(npix, 1.0f); + + // expect sinc in time domain, centered at [0, 0] since no phase above. + auto result = ft.transform_host(sizes, {sizes[1], 1}, spectrum.data()); + auto view = NFFT2dResultView(result); + auto results_d = thrust::device_vector>(1); + + interp<<<1, 1>>>(view, results_d.data().get()); + + checkCudaErrors(cudaPeekAtLastError()); + checkCudaErrors(cudaDeviceSynchronize()); + + thrust::host_vector> results_h = results_d; + auto z0 = results_h[0]; + + EXPECT_NEAR(z0.real(), 1.0f, 1e-4); + EXPECT_NEAR(z0.imag(), 0.0f, 1e-4); +} + +int main(int argc, char* argv[]) +{ + testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} \ No newline at end of file diff --git a/tests/cxx/isce3/geometry/geometry/geometry.cpp b/tests/cxx/isce3/geometry/geometry/geometry.cpp index 82db903f7..040575c92 100644 --- a/tests/cxx/isce3/geometry/geometry/geometry.cpp +++ b/tests/cxx/isce3/geometry/geometry/geometry.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -304,6 +305,19 @@ TEST(Geometry, SrLkvHeadDemNed) << "Wrong S/C Vel in ENU"; } +TEST(Geometry, GeoToPolar) +{ + double sinSquint = 100.0, range = 0.0; + isce3::core::Vec3 target_xyz{-1, 0, 0}, origin{0, 0, 1}, axis{0, 1, 0}; + auto status = isce3::geometry::geo2polar(&sinSquint, &range, target_xyz, + origin, axis); + + EXPECT_EQ(status, isce3::error::ErrorCode::Success); + const double atol = 1e-14; + EXPECT_NEAR(sinSquint, 0.0, atol); + EXPECT_NEAR(range, std::sqrt(2.0), atol); +} + int main(int argc, char* argv[]) { testing::InitGoogleTest(&argc, argv); diff --git a/tests/cxx/isce3/signal/nfft.cpp b/tests/cxx/isce3/signal/nfft.cpp index 9af148b2f..4ca9b9a58 100644 --- a/tests/cxx/isce3/signal/nfft.cpp +++ b/tests/cxx/isce3/signal/nfft.cpp @@ -6,6 +6,7 @@ #include #include #include "isce3/signal/NFFT.h" +#include "isce3/signal/NFFT2d.h" #include "isce3/signal/Filter.h" const int seed = 1234; @@ -155,6 +156,33 @@ TEST(Kernel, Singularity) EXPECT_GT(window(m), window(m+dx)); } +TEST(NFFT2d, IRF) +{ + using T = float; + using dims_t = isce3::signal::NFFT2d::dims_t; + dims_t dims = {32, 84}; + constexpr int s = 2, my = 4, mx = 4; + dims_t fft_dims = {dims[0] * s, dims[1] * s}; + auto nfft = isce3::signal::NFFT2d({my, mx}, dims, fft_dims); + + // Set spectrum to all ones. + size_t nimg = static_cast(dims[0]) * dims[1]; + std::vector> z(nimg); + z.assign(nimg, std::complex(1.0, 0.0)); + + // Transform is the impulse response, which for NFFT should be approximately + // a sinc. + const auto result = nfft.transform(dims, {dims[1], 1}, z.data()); + + // Zero phase on spectrum, so IRF should be centered at (0, 0). + const auto Z0 = result.interp({0.0, 0.0}, true); + + const auto err_real = std::abs(std::real(Z0) - 1.0); + const auto err_imag = std::abs(std::imag(Z0)); + EXPECT_LT(err_real, 1e-4); + EXPECT_LT(err_imag, 1e-4); +} + int main(int argc, char *argv[]) { diff --git a/tests/python/extensions/pybind/CMakeLists.txt b/tests/python/extensions/pybind/CMakeLists.txt index 1af56362f..ddcb91b9e 100644 --- a/tests/python/extensions/pybind/CMakeLists.txt +++ b/tests/python/extensions/pybind/CMakeLists.txt @@ -24,6 +24,7 @@ focus/presum.py focus/rangecomp.py geocode/geocodeCov.py geocode/radar_grid_cube.py +geogrid/geogrid_ecef_coords.py geometry/bbox.py geometry/dem.py geometry/geo2rdr.py @@ -42,6 +43,7 @@ signal/crossmul.py signal/crossmultiply.py signal/filter2D.py signal/multilook.py +signal/nfft2d.py product/generic_product.py product/geogridparameters.py product/radargridparameters.py diff --git a/tests/python/extensions/pybind/focus/backproject.py b/tests/python/extensions/pybind/focus/backproject.py index 5b671002f..c4e6aed94 100644 --- a/tests/python/extensions/pybind/focus/backproject.py +++ b/tests/python/extensions/pybind/focus/backproject.py @@ -4,7 +4,9 @@ import numpy as np import numpy.testing as npt import isce3.ext.isce3 as isce +import isce3 from isce3.core import load_orbit_from_h5_group +from isce3.focus.serialization import BackprojectionStageParameters from iscetest import data as test_data_dir from pathlib import Path import json @@ -70,6 +72,7 @@ def load_h5(filename): "target_azimuth": target_azimuth, "target_range": target_range} + def test_backproject(): # load point target simulation data filename = Path(test_data_dir) / "point-target-sim-rc.h5" @@ -168,3 +171,136 @@ def test_backproject(): # threshold is slightly higher - see # https://github.jpl.nasa.gov/bhawkins/nisar-notebooks/blob/master/Azimuth%20Resolution.ipynb assert(azimuth_width <= 6.62) + + +class DummyWriter: + def __init__(self, shape): + self.shape = shape + self.data = np.zeros(shape, dtype="c8") + + def queue_write(self, z, block): + self.data[block] = z + + +# Copy/paste of existing BP test for FBP. Uses the higher-level azcomp_fbp +# interface instead of the pybind11 bindings directly, which would involve +# even more copy/paste... +def test_azcomp_fbp(): + # load point target simulation data + filename = Path(test_data_dir) / "point-target-sim-rc.h5" + d = load_h5(filename) + + # eww gross + signal_data = d["signal_data"] + radar_grid = d["radar_grid"] + orbit = d["orbit"] + doppler = d["doppler"] + range_sampling_rate = d["range_sampling_rate"] + dem = d["dem"] + dry_tropo_model = d["dry_tropo_model"] + target_azimuth = d["target_azimuth"] + target_range = d["target_range"] + + # range bandwidth (Hz) + B = 20e6 + + # desired azimuth resolution (m) + azimuth_res = 6. + + # output chip size + nchip = 129 + + # how much to upsample the output for point target analysis + upsample_factor = 128 + + # create 9-point Knab kernel + # use tabulated kernel for performance + kernel = isce.core.KnabKernel(9., B / range_sampling_rate) + kernel = isce.core.TabulatedKernelF32(kernel, 2048) + + # create output radar grid centered on the target + dt = radar_grid.az_time_interval + dr = radar_grid.range_pixel_spacing + t0 = target_azimuth - 0.5 * (nchip - 1) * dt + r0 = target_range - 0.5 * (nchip - 1) * dr + out_grid = isce.product.RadarGridParameters( + t0, radar_grid.wavelength, radar_grid.prf, r0, dr, + radar_grid.lookside, nchip, nchip, orbit.reference_epoch) + + # init output buffer + out = np.empty((nchip, nchip), np.complex64) + # and debug height layer + height = np.empty(out.shape, np.float32) + + # collect input & output radar_grid, orbit, and Doppler + in_geometry = isce.container.RadarGeometry(radar_grid, orbit, doppler) + out_geometry = isce.container.RadarGeometry(out_grid, orbit, doppler) + + fbp_factors = [ + BackprojectionStageParameters(size=64), + BackprojectionStageParameters(size=2), + ] + blocks_bounds = [ + ( + (slice(None), slice(None)), + (radar_grid.sensing_start, radar_grid.sensing_stop), + ), + ] + writer = DummyWriter(out_grid.shape) + + isce3.focus.azcomp_bp.azcomp_fbp(fbp_factors, azimuth_res, kernel, + blocks_bounds, in_geometry, signal_data, out_grid, writer, + height=height, dem=dem, atmos=dry_tropo_model, bandwidth=B) + + out[...] = writer.data + + # We used a constant DEM height, so make sure the debug height layer + # contains that value everywhere. + npt.assert_allclose(height, dem.ref_height) + + # remove range carrier + kr = 4. * np.pi / out_grid.wavelength + r = np.array(out_geometry.slant_range) + out *= np.exp(-1j * kr * r) + + info, _ = analyze_point_target(out, nchip//2, nchip//2, nov=upsample_factor, + chipsize=nchip//2) + tofloatvals(info) + + # print point target info + print(json.dumps(info, indent=2)) + + # range resolution (m) + range_res = c / (2. * B) + + # range position error & -3 dB main lobe width (m) + range_err = dr * info["range"]["offset"] + range_width = dr * info["range"]["resolution"] + + # azimuth position error & -3 dB main lobe width (m) + _, vel = orbit.interpolate(target_azimuth) + azimuth_err = dt * info["azimuth"]["offset"] * np.linalg.norm(vel) + azimuth_width = dt * info["azimuth"]["resolution"] * np.linalg.norm(vel) + + # require positioning error < resolution/128 + assert(range_err < range_res / 128.) + assert(azimuth_err < azimuth_res / 128.) + + # require 3dB width in range to be <= range resolution + assert(range_width <= range_res) + + # azimuth response is spread slightly by the antenna pattern so the + # threshold is slightly higher - see + # https://github.jpl.nasa.gov/bhawkins/nisar-notebooks/blob/master/Azimuth%20Resolution.ipynb + assert(azimuth_width <= 6.62) + + +def test_polar_grid_hash(): + LS = isce.core.Linspace + args = (0.0, 1.0, [0, 0, 0], [1, 0, 0], LS(0, 1, 2), LS(0, 1, 2), + isce.core.LookSide.Left) + # Different instances with same parameters should be equal. + grid1 = isce.focus.PolarGrid(*args) + grid2 = isce.focus.PolarGrid(*args) + assert hash(grid1) == hash(grid2) + assert grid1 == grid2 diff --git a/tests/python/extensions/pybind/geogrid/geogrid_ecef_coords.py b/tests/python/extensions/pybind/geogrid/geogrid_ecef_coords.py new file mode 100644 index 000000000..463fe7d93 --- /dev/null +++ b/tests/python/extensions/pybind/geogrid/geogrid_ecef_coords.py @@ -0,0 +1,39 @@ +import isce3.ext.isce3 as isce3 +import numpy as np +import numpy.testing as npt +from pyproj import Proj, CRS, Transformer + +def test_ecef_coords(): + x0, y0 = 402538, 3851590 + dx, dy = 1, -1 + m, n = 5, 3 + epsg = 32611 + + grid = isce3.product.GeoGridParameters(x0, y0, dx, dy, n, m, epsg) + + dem = isce3.geometry.DEMInterpolator() + href = 1000. + dem.ref_height = href + + xyz = isce3.geogrid.get_geogrid_ecef_coords(grid, dem) + + assert xyz.shape == (m, n, 3) + assert xyz.dtype == np.float64 + + # check against pyproj transform + xform = Transformer.from_crs( + CRS.from_epsg(epsg).to_3d(), + CRS.from_epsg(4978).to_3d()) + + xyz_ref = np.zeros_like(xyz) + for i in range(m): + v = y0 + dy * i + for j in range(n): + u = x0 + dx * j + xyz_ref[i, j, :] = xform.transform(u, v, href) + if not np.allclose(xyz_ref[i, j], xyz[i, j]): + print(f"mismatch xyz[{i},{j}]:") + for k in range(3): + print(f" {xyz[i,j,k]:.3f} vs {xyz_ref[i,j,k]:.3f}") + + npt.assert_allclose(xyz, xyz_ref) diff --git a/tests/python/extensions/pybind/signal/nfft2d.py b/tests/python/extensions/pybind/signal/nfft2d.py new file mode 100644 index 000000000..c62999a01 --- /dev/null +++ b/tests/python/extensions/pybind/signal/nfft2d.py @@ -0,0 +1,54 @@ +import isce3.ext.isce3 as isce3 +import numpy as np +import numpy.testing as npt + + +def make_delta_test_image(shape): + z = np.zeros(shape, "c8") + z[0, 0] = 1.0 + return z + + +def test_make_nfft2d(): + # Input is Kronecker delta. + image = make_delta_test_image((128, 365)) + + # NFFT parameters for both dimensions. + params = {'rows': {'m': 4, 's': 4.0}, 'cols': {'m': 4, 's': 4.0}} + + # Use convenience function to create an interpolator. + itp = isce3.signal.make_image_nfft2d(image, params) + + # Interp at origin. + z0 = itp.interp((0.0, 0.0)) + + npt.assert_allclose(z0, 1.0+0j, atol=1e-4) + + +def test_nfft2d(): + # Input is Kronecker delta. + image = make_delta_test_image((128, 365)) + + # Compute its spectrum. + image_spectrum = np.fft.fft2(image) + shape = image_spectrum.shape + + # NFFT parameters for both dimensions. + m = s = 4 + + # Plan 2D NFFT + nfft = isce3.signal.NFFT2dF32((m, m), shape, [n * s for n in shape]) + + # Execute to get time-domain image ready for interpolation. + itp = nfft.transform(image_spectrum) + + # Interpolate at origin. + z0 = itp.interp((0.0, 0.0)) + + # Check result + npt.assert_allclose(z0, 1.0+0j, atol=1e-4) + + # Check bindings + npt.assert_equal(nfft.sizes, shape) + npt.assert_equal(nfft.fft_sizes, [s * n for n in shape]) + npt.assert_equal(nfft.spectrum.shape, nfft.fft_sizes) \ No newline at end of file diff --git a/tests/python/packages/nisar/workflows/focus.py b/tests/python/packages/nisar/workflows/focus.py index a4e63e6e4..96efdf9ea 100644 --- a/tests/python/packages/nisar/workflows/focus.py +++ b/tests/python/packages/nisar/workflows/focus.py @@ -5,6 +5,7 @@ from nisar.workflows.point_target_analysis import slc_pt_performance import nisar from pathlib import Path +import pytest import numpy as np import numpy.testing as npt import os @@ -33,8 +34,14 @@ def slc_is_baseband(filename: str, tol=2*np.pi/100, frequency="A", polarization= return abs(np.angle(dz.sum())) < tol -def test_focus(): +direct_bp = [focus.Struct({"size": 1})] +factorized_bp_64 = [focus.Struct({"size": 64})] + +@pytest.mark.parametrize("factorization", (direct_bp, factorized_bp_64)) +def test_focus(factorization): cfg = get_test_cfg() + cfg.runconfig.groups.processing.azcomp.factorization = factorization + focus.focus(cfg) filename = cfg.runconfig.groups.product_path_group.sas_output_file