Conversation
This reverts commit 5ca283a.
This reverts commit 733479f.
Change accumulatePolarImagesToRadarGrid and mergePolarImages to use std::vector<const NFFT2dResult<float>*> instead of std::vector<NFFT2dResult<float>> to match the GPU implementation and avoid unnecessary copies. Updated Python bindings to accept py::sequence and convert to vector of pointers, consistent with the GPU bindings. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
|
All the CI failures are in Seems unrelated to the PR and the test passes locally. |
bhawkins
left a comment
bhawkins
left a comment
There was a problem hiding this comment.
Comments from today's review meeting.
Documented all functions and classes in serialization.py with NumPy format docstrings including parameter types, return values, and attribute descriptions. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
|
Fix for the CI failure is in #348 |
|
I just pushed a change making FBP the default algorithm in order to simplify PCM/SDS testing. @hfattahi @jshimada47 |
|
I'm getting a unit test failure while doing the build on dev-5 in |
Tyler-g-hudson
left a comment
Tyler-g-hudson
left a comment
There was a problem hiding this comment.
I'm slowly chewing my way through the PR, here are a very basic review item and question
| # List of aperture factorization parameters. | ||
| # Default to two-stage with 128 pulses per subaperture. | ||
| factorization: | ||
| - size: 128 |
There was a problem hiding this comment.
To clarify, do we intend to modify this list to be hierarchical later on for more efficient processing or no?
Also, for a typical input L0B ~40 seconds per RSLC containing ~76k pulses, wouldn't 256 be the nearest power of 2 to the optimum per Yegulalp? Have you had a chance to check this case? Are my numbers wrong here? (I'm basing them off of a single full-length RSLC I have, and assuming a 1.5 second pad before and after in the L0B at the L0B native PRF)
Do we need additional scaled testing for hierarchical runs? I could produce some scaled workflows on a g6e if you have one or more RSLCs you want me to reproduce and a handful of factorization lists.
There was a problem hiding this comment.
Recall that we're using stripmap mode (not spotlight), so there are only around 6000 pulses in a NISAR synthentic aperture. Then Yegulalp's sqrt(num_aperture) suggestion for a two-stage algorithm would work out to 77 pulses.
However, in the version here (at least in the default config) the first stage uses a 9-point 1D interpolator while the second (and higher) stage(s) use a 5x5=25 point 2D interpolator. So if interpolation dominates the runtime, then later stages are around f=2.7 times more expensive than the first stage. To split the work evenly among a two-stage algorithm, then you'd want to solve num_aperture = (f * n) * n which results in around f * n = 129 pulses in the first stage. In the general case you'd want to solve
num_aperture = f**(num_stages - 1) * n**num_stagesfor n. So for a NISAR 3-stage algorithm you might try the config
factorization:
- size: 70
- size: 9However, there are other considerations like the geometry calculations (including the fact that the geometry is in float64 while the interpolation is float32), the FFT overhead, GPU memory bandwidth, cache sizes, and so-on that complicate the picture. I measured runtimes for two-stage cases using size in {64, 128, 256} and didn't see very large differences between them on g6/g6e instances. It might be worth repeating that experiment since I've hacked on the code a fair amount since then.
py_image_interpolators elements were only borrowed via a transient py::object during std::transform, leaving no owning reference once the GIL was released in accumulate_polar_images_to_radar_grid. Another thread could drop the last reference to an element while the C++ code held a raw pointer into it. Materialize an owning std::vector<py::object> before releasing the GIL, and apply the same pattern in merge_polar_images for consistency. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Now that image_interpolators elements are held via owning py::object references, it's safe to release the GIL during the potentially long-running mergePolarImages call, matching the pattern already used in accumulate_polar_images_to_radar_grid. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Tyler-g-hudson
left a comment
Tyler-g-hudson
left a comment
There was a problem hiding this comment.
Some questions and naming/documentation recommendations
| # 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. |
There was a problem hiding this comment.
Would Python-level parallelism even help? The individual tasks are parallel at the C++/OMP level so, at least for CPU, they should saturate the available computation resources, no?
This would probably be more complicated to figure out for GPU, but you'd need to see if the GPU usage is low enough at any step of the task to justify parallel task requests and associated complexity. Maybe worth a test later.
| _, img, _ = bp_to_polar_grid(fdata, sr, x, v, | ||
| polar_grid, dem, fc, kernel, "nodelay", rdr2geo_params) | ||
| if debugfile is not None: | ||
| import h5py |
There was a problem hiding this comment.
Recommend moving this to the top of the file, I don't think this is a heavy dependency relative to ISCE3 and should only be performed once as the file is uploaded if it has not already been imported (and it likely has)
| template <typename T> | ||
| class NFFT2dResult { |
There was a problem hiding this comment.
Readability suggestion:
| template <typename T> | |
| class NFFT2dResult { | |
| /** | |
| * Prepared interpolator for 2D image queries at arbitrary positions. | |
| * Holds oversampled grid enabling rapid Kaiser-Bessel interpolation | |
| * queries. | |
| */ | |
| template <typename T> | |
| class KaiserBesselInterpolator2d { |
I think this could use a rename and/or or some nearby documentation. My reason for suggesting this is that the data it holds is in the polar domain (Seeing NFFT in the name might make a reader mistake at a glance that this as some sort of frequency spectrum) and that, in principle, an oversampler of any kind could produce one of these interpolators (presumably in any domain?) so long as it can also produce the correct kernels.
In an ideal world I'd recommend moving this to a separate file and renaming it - failing those, a quick comment here would probably help future maintenance.
I'm interested in this method in general - do you think it's practical for InSAR resampling? Do you know any pros and cons to doing so off the top of your head?
| /** | ||
| * @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<T> object for interpolating the image. | ||
| */ | ||
| template<typename T> | ||
| NFFT2dResult<T> makeImageNFFT2d( |
There was a problem hiding this comment.
Same comment as above - this no longer returns an NFFT2d so I recommend a rename for it. Docstring parameters also look to have been changed.
| /** | |
| * @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<T> object for interpolating the image. | |
| */ | |
| template<typename T> | |
| NFFT2dResult<T> makeImageNFFT2d( | |
| /** | |
| * @brief Create an rapid interpolator using an intermediate NFFT2d. | |
| * | |
| * @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 params parameters for the intermediate NFFT | |
| * @param pad_input Whether to also zero-pad input data to an efficient | |
| * forward transform size. Requires extra memory. | |
| * | |
| * @return KaiserBesselInterpolator2d<T> object for interpolating the image. | |
| */ | |
| template<typename T> | |
| KaiserBesselInterpolator2d<T> makeFastImageInterpolator( |
4 participants
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
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
which means sum groups of 128 pulses into polar grids before summing those into the output image. Similarly,
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