Factorized backprojection

[bhawkins]

This PR adds a factorized backprojection (FBP) algorithm for SAR focusing as a faster alternative to the existing direct backprojection method.

Algorithm Design

The basic idea is from Yegulalp 1999 where local, low-resolution polar grids can efficiently store partial sums of subapertures. As in other papers, the idea is extended hierarchically, so that multiple local grids can be combined together into a higher resolution grid. This process can be repeated in successive stages until the desired resolution is met. In the limit one achieves a factor of $O(N / log(N))$ reduction in the number of sums per output pixel (where $N$ is the number of pulses in the synthetic aperture). In practice there's a fair amount of overhead, but the speedup is still significant.

The FBP algorithm is usually described for the spotlight case, and the one here has a few modifications for the stripmap case relevant to NISAR. First, the azimuth extent of the local grid is extended based on the duration of the subaperture, since the beam is sliding. Second, the azimuth resolution of the local grid is not allowed to exceed the L/2 limit. The radar travel time between transmit and receive ("bistatic" correction) is accounted for both in the initial back projection stage and in the selection of grid coordinates. The latter is necessary to maintain a baseband signal in azimuth.

I'm not sure I'd call this a "fast factorized" (FFBP) algorithm. In the interest of maximal accuracy, it does full 2D interpolation and makes no assumptions about the trajectory or nesting of coordinate grids. Non-uniform FFT algorithms are used for interpolation in order to maintain high accuracy with small kernels. The zero-padded transforms do occupy a fair amount of memory, though, which may be a bottleneck. Care has been taken to ensure that only the "active" subimages reside in memory at any given time.

You might notice a few new functions that are not used in the NISAR RSLC workflow. These have been used for processing UAVSAR data from circular tracks (CSAR).

Key Features

• Factorized backprojection algorithm: Multi-stage approach (initial → middle → final) that reduces computational complexity compared to direct backprojection
• 2D NFFT implementation: New CPU and GPU implementations of 2D Non-uniform FFT (NFFT2d) required by the FBP algorithm
• GPU acceleration: CUDA implementations for all FBP stages and NFFT2d operations
• Bug fixes in interp1d and interp2d
• Python interface: New azcomp_bp and azcomp_fbp functions exposed via isce3.focus.azcomp_bp module
• Configuration support: Added FBP parameters to focus workflow configuration schema

Usage

The algorithm is highly configurable. You can specify the number of stages, the size of each stage, and the interpolation parameters for each stage (NFFT kernel size and oversampling ratio in both dimensions). To use a two-stage algorithm, the RSLC configuration file should look like

runconfig:
    groups:
        processing:
            azcomp:
                factorization:
                - size: 128

which means sum groups of 128 pulses into polar grids before summing those into the output image. Similarly,

# runconfig.groups.processing.azcomp
                factorization:
                - size: 64
                - size: 2

means a three-stage algorithm where the initial grids comprise 64 pulses each, these get merged together two at a time, and then the merged grids get summed into the output image. The default factorization: [{"size": 1}] is understood to mean the original direct backprojection algorithm.

In my testing I've found the two-stage n=128 algorithm is a 10x speedup of the azcomp algorithm on a g6e instance in AWS EC2.

I included a script that can generate animations like the following from the FBP debug file:

factors.mp4

[bhawkins]

Hmm the GPU build is failing, though it builds fine in my dev environment. I notice that I'm on CUDA 12 while the CI is building with CUDA 13. It looks like there was a change regarding visibility of __global__ symbols (CUDA kernels) that might be relevant (link). It's a little tricky since I used a template to allow it to use Python containers without copying the contents, and I'd like to avoid instantiating Python junk in the C++ API. Will chew on this a bit more...

[bhawkins]

All the CI failures are in test.cxx.isce3.io.raster.raster with

[ RUN      ] RasterTest.addRasterToVRT
ERROR 3: Failed to read scanline 100.
ERROR 1: IReadBlock failed at X offset 0, Y offset 100: Failed to read scanline 100.
In isce3::io::Raster::getSetValue() - error in RasterIO.
/home/runner/work/isce3/isce3/tests/cxx/isce3/io/raster/raster.cpp:258: Failure
Expected equality of these values:
  std::isfinite(val)
    Which is: false
  true

[  FAILED  ] RasterTest.addRasterToVRT (55 ms)

Seems unrelated to the PR and the test passes locally.

[bhawkins]

Comments from today's review meeting.

[bhawkins]

Fix for the CI failure is in #348