Factorized backprojection

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()

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

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

cxx/isce3/core/Interp2d.icc

@@ -33,8 +33,13 @@ DataType interp2d(const Kernel<KernelType>& 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<long>(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(

cxx/isce3/core/Kernels.h

@@ -103,6 +103,10 @@ class NFFTKernel : public Kernel<T> {
 
     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;

cxx/isce3/core/Linspace.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "Common.h"
+#include <array>
 
 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<T, 2> bounds() const;
+
 private:
     T _first = {};
     T _spacing = {};

cxx/isce3/core/Linspace.icc

@@ -106,4 +106,13 @@ bool operator!=(const Linspace<T> & lhs, const Linspace<U> & rhs)
     return !(lhs == rhs);
 }
 
+template<typename T>
+CUDA_HOSTDEV
+constexpr
+std::array<T, 2UL> Linspace<T>::bounds() const
+{
+    const T half = spacing() / 2;
+    return std::array<T, 2UL>{first() - half, last() + half};
+}
+
 }}

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<long>(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 {

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

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
 )

cxx/isce3/cuda/Sources.cmake

@@ -53,4 +53,5 @@ signal/gpuFilter.cu
 signal/gpuLooks.cu
 signal/gpuRangeFilter.cu
 signal/gpuSignal.cu
+signal/NFFT2d.cu
 )

cxx/isce3/cuda/core/Interp2d.h

@@ -0,0 +1,37 @@
+#pragma once
+
+#include "forward.h"
+#include <isce3/core/Common.h>
+
+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<typename KernelX, typename KernelY, typename DataType>
+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"

cxx/isce3/cuda/core/Interp2d.icc

@@ -0,0 +1,139 @@
+#include <isce3/math/complexOperations.h>
+
+#include "Interp1d.h"
+
+namespace isce3::cuda::core {
+
+namespace detail {
+    template<typename Kernel, typename KT>
+    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<long>(ceil(t));
+        } else {
+            i0 = static_cast<long>(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<typename DataType>
+    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<long>(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<DataType>(0);
+                }
+            }
+        }
+        return block;
+    }
+
+    template<typename TX, typename TY>
+    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<TX, TY>::type;
+        TO sum = 0;
+
+        for (int i = 0; i < width; ++i) {
+            sum += x[i] * y[i];
+        }
+        return sum;
+    }
+}
+
+template<typename KernelX, typename KernelY, typename DataType>
+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<int>(ceil(kernelx.width()));
+    const int widthy = static_cast<int>(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<TX, TY>::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<long>(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);
+}
+
+}

cxx/isce3/cuda/core/Kernels.h

@@ -125,6 +125,60 @@ class KnabKernel : public Kernel<T, KnabKernel<T>> {
     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<typename T>
+class NFFTKernel : public Kernel<T, NFFTKernel<T>> {
+    using Base = Kernel<T, NFFTKernel<T>>;
+    friend Base;
+
+public:
+    /** A non-owning kernel view type that can be passed to device code */
+    using view_type = NFFTKernel<T>;
+
+    /**
+     * 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<T>& 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<T>() 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<typename T>
 class TabulatedKernelView : public Kernel<T, TabulatedKernelView<T>> {