NUFFT forward returns NaN for certain near-edge coordinates at default `oversamp` (Kaiser–Bessel gridding)

[sabbir-ahmed12]

Summary

sigpy.nufft (and therefore sigpy.linop.NUFFT) returns NaN for certain
valid, finite, in-range k-space coordinates when run on the GPU (cupy) at the
default oversamp=1.25. The input image and the coordinate are both finite,
and the coordinate is well inside the sampled region, yet a single output sample
comes back NaN. Increasing oversamp (e.g. to 2.0) makes the same call
finite.

This is a numerical edge case in the Kaiser–Bessel gridding kernel that triggers
when a coordinate maps essentially onto an integer line of the oversampled
grid. In a 3-D radial/FLORET acquisition it reliably poisons the affected
readout sample, which then propagates (e.g. into an iterative reconstruction).

Environment

• OS: CentOS 7
• sigpy 0.1.26 (verified: sigpy/interp.py is identical in the
  latest release 0.1.27 and on main, so this is not fixed upstream)
• cupy 10.6.0 (cupy-cuda114), CUDA 11.4
• Python 3.9.18, NumPy 1.24.4
• GPU: NVIDIA A100

Minimal reproducer (requires a CUDA GPU + cupy)

import numpy as np
import sigpy as sp

gpu = sp.Device(0)          # CUDA device; the bug is on the cupy/CUDA path
xp = gpu.xp
N = 256

# A real, finite, in-range coordinate from a 3-D FLORET trajectory.
# radius = 113.46 < N/2 = 128, so it is a valid sample for a 256^3 grid.
coord = np.array([[-108.42059147846787, 21.51759054102146, -25.60000000000001]])

with gpu:
    img = xp.ones((N, N, N), dtype=xp.complex64)
    cg = sp.to_device(coord, gpu)
    for oversamp in (1.25, 2.0):
        k = sp.nufft(img, cg, oversamp=oversamp)
        print(f"oversamp={oversamp}: output={complex(k.ravel()[0]):.4g}  "
              f"all_finite={bool(xp.isfinite(k).all())}")

Actual output

Confirmed on an NVIDIA A100 (sigpy 0.1.26, cupy-cuda114 10.6.0, CUDA 11.4):

oversamp=1.25: output=nan+nanj  all_finite=False
oversamp=2.0: output=0.0005427-0.002819j  all_finite=True
versions: sigpy 0.1.26, numpy 1.24.4

Expected output

A finite value at both oversampling ratios (the coordinate is valid at either).

Additional Analysis

Decomposing the forward NUFFT for this coordinate shows the apodization and the
oversampled FFT are both finite; the NaN is created in the interpolation gridding step, and it appears even when gridding an all-ones grid — so it is a
property of the kernel weights, not of the data being interpolated.

The tell is the coordinate's position on the oversampled grid. With N=256 and
oversamp=1.25 the oversampled size is 320, and _scale_coord maps the
-25.6 component to -25.6 * (320/256) + 160 = 128.0 — i.e. exactly on an
integer grid line. Changing oversamp moves the coordinate off that line and
the result becomes finite.

The Kaiser–Bessel kernel (sigpy/interp.py) does guard the boundary, in both
the numba and the CUDA variants — so this is not a naive unguarded sqrt:

// _kaiser_bessel_kernel_cuda
__device__ inline S kernel(S x, S beta) {
    if (fabsf(x) > 1) return 0;        // boundary guard
    x = beta * sqrt(1 - x * x);        // sqrt of a negative value would give NaN
    ...
}

Hypothesis (not proven):

The guard fabsf(x) > 1 is evaluated in single precision, while the coordinates and 1 - x*x are effectively double precision for a float64 coord. A value by a rounding-error-sized amount 1 can therefore pass the float guard yet still make 1 - x*x slightly negative, so sqrt(1 - x*x) returns NaN. This is consistent with all observations (only coordinates landing on grid lines fail; changing oversamp fixes it), but I have not instrumented the compiled kernel to confirm the exact operation, so I'm flagging it as a strong hypothesis rather than a certainty.

Suggested fix

Make the kernel robust regardless of which coordinates land on grid lines, e.g.:

• clamp the argument before the square root: sqrt(fmax(0, 1 - x*x)) (and the
  numba equivalent), and/or
• make the boundary guard use the same precision as the coordinate/x so the
  guard and the sqrt cannot disagree.