`pygrappa.grappa` may leave k-space holes unfilled with NumPy 2.3.0 and later

[Arctic-Xiangjian]

Summary

The pygrappa.grappa implementation can silently leave many missing k-space locations unfilled when used with NumPy 2.3.0 and later.

This is caused by an upstream NumPy regression reported in:

• BUG: unravel_index returns incorrect coordinates for large (N, 1) index arrays across an iterator buffer boundary numpy/numpy#31980

In the affected NumPy versions, np.unravel_index can return incorrect and duplicated coordinates when given a sufficiently large integer array with shape (N, 1). The current grappa implementation passes exactly this type of array to np.unravel_index:

idx = np.unravel_index(
    np.argwhere(iidx == ii),
    Psh[:2],
)

For large sampling-pattern groups, some target coordinates are duplicated while other targets are never visited. The reconstruction completes without an exception, but some missing k-space samples remain zero.

Affected versions

The upstream issue has been reproduced with every tested stable NumPy release beginning with NumPy 2.3.0.

Last known good NumPy release: 2.2.6
First known bad NumPy release: 2.3.0

The issue is not specific to one Python version.

Impact on pygrappa.grappa

The relevant code is currently:

idx = np.unravel_index(
    np.argwhere(iidx == ii),
    Psh[:2],
)

x, y = idx[0] + kx2, idx[1] + ky2
x = np.atleast_1d(x.squeeze())
y = np.atleast_1d(y.squeeze())

For a one-dimensional iidx array:

np.argwhere(iidx == ii)

returns an array with shape:

(N, 1)

On affected NumPy versions, sufficiently large (N, 1) inputs can produce incorrect coordinates after an internal iterator-buffer boundary.

As a result, the following reconstruction loop may write repeatedly to duplicated coordinates:

for xx, yy in zip(x, y):
    ...
    recon[xx, yy, :] = ...

while some intended target coordinates are never written.

Observed behavior includes:

• the reconstruction returns normally;
• acquired k-space samples remain unchanged;
• some missing samples remain exactly zero;
• the resulting image may contain residual aliasing or other reconstruction artifacts.

Suggested fix

The indices should be converted to a one-dimensional flat-index array before calling np.unravel_index.

A direct fix is:

flat_indices = np.flatnonzero(iidx == ii)

x, y = np.unravel_index(
    flat_indices,
    Psh[:2],
)

x = x + kx2
y = y + ky2

The surrounding code would become:

flat_indices = np.flatnonzero(iidx == ii)

x, y = np.unravel_index(
    flat_indices,
    Psh[:2],
)

x = x + kx2
y = y + ky2

for xx, yy in zip(x, y):
    S = kspace[
        xx-kx2:xx+kx2+adjx,
        yy-ky2:yy+ky2+adjy,
        :,
    ]
    S = S[P[ii, ...]]
    recon[xx, yy, :] = (W @ S[:, None]).squeeze()

This also removes the need for:

x = np.atleast_1d(x.squeeze())
y = np.atleast_1d(y.squeeze())

because np.flatnonzero always returns a one-dimensional array.

An equivalent minimal change would be:

idx = np.unravel_index(
    np.argwhere(iidx == ii).ravel(),
    Psh[:2],
)

However, np.flatnonzero(iidx == ii) more clearly expresses that these values are flat indices.

Why this change is appropriate independently of the NumPy fix

Although the incorrect coordinate values are caused by the upstream NumPy regression, np.flatnonzero is also the more appropriate API for this operation.

The code needs the one-dimensional positions where:

iidx == ii

is true. np.flatnonzero returns exactly those flat positions, while np.argwhere returns a coordinate table with shape (N, iidx.ndim).

Using np.flatnonzero therefore:

• avoids the affected NumPy code path;
• makes the intended flat-index semantics explicit;
• simplifies the coordinate-handling code;
• remains compatible with older NumPy versions.

Suggested regression test

A regression test should create a sampling-pattern group containing more than 8193 targets and verify that every expected target coordinate is visited exactly once.

The strongest test would track target assignment counts directly rather than inferring assignment from reconstructed values:

assignment_count = np.zeros(Psh[:2], dtype=np.int64)

flat_indices = np.flatnonzero(iidx == ii)
x, y = np.unravel_index(flat_indices, Psh[:2])

assignment_count[x, y] += 1

The test should verify that:

assert np.all(assignment_count[expected_targets] == 1)

A full-size Cartesian mask, such as a 384 × 384 mask with a fixed central ACS region, can produce pattern groups large enough to exercise this case.

Upstream reference

NumPy issue:

• BUG: unravel_index returns incorrect coordinates for large (N, 1) index arrays across an iterator buffer boundary numpy/numpy#31980

The NumPy-only reproducer demonstrates that identical flat-index values produce different coordinates depending on whether they are supplied with shape (N,) or (N, 1).