How to ensure the final results similar to the initial input after a type-1 transformation and a subsequent type-2 transformation?

[zongjg]

Hi All,

Does anyone have experiences in using FINUFFT in non-cartesian MRI?

I'm designing an algorithm for radial-MRI, and there is an essential step which is called Data-Consistency Constraint, it means that the results of an algorithm must be similar to the original input.

But when I use FINUFFT to perform it, I can't get the expected results, the results are very different from the original input.

So, does anyone have some experiences in it?
How should I scale the results or configure the finufft plan?

Thanks very much for your help!

---

Please review the following matlab code!
Version of FINUFFT: 2.2.0

clear;

% init param
nRO = 512; % samples
nSpk = 1000; % how many radial-lines in a plane
nRowCol = nRO/2; % row and column of the reconstructed images

% create an radial trajectory by a group of golden angles
golden_ratio = 0.61803398874989484820458683436564;
golden_angle = 2 * pi * (1 - golden_ratio);
angles = golden_angle .* (1:nSpk);
k_rad = linspace(-pi, pi, nRO);
[kr_grid, theta_grid] = meshgrid(k_rad, angles);
kx = single(kr_grid.' .* cos(theta_grid.'));
ky = single(kr_grid.' .* sin(theta_grid.'));

% create dcf for density compensation
dcf_all = sqrt(kx.^2 + ky.^2);
dcf_all(dcf_all < 1e-6) = 1;
dcf_all = dcf_all / max(dcf_all(:));

%create the initial images and radial sampled lines

noisectrl = 0.00;
imgs_init = single(phantom(nRowCol)+rand(nRowCol)*noisectrl)*1000;

spokes_init = finufft2d2(kx(:), ky(:), -1, 1e-15, imgs_init);
spokes_init = reshape(spokes_init./nSpk,[nRO, nSpk]);

% to perform the density compensation
c1_dcomp = spokes_init.*dcf_all;

% 1st type-1 transform
imgs_1 = finufft2d1(kx(:), ky(:), c1_dcomp, 1, 1e-15, nRowCol, nRowCol);
imgs_1 = imgs_1./nRowCol; %!!! Here I expect imgs_1 equal or similar to imgs_init

% 1st type-2 transform
c2 = finufft2d2(kx(:), ky(:), -1, 1e-15, imgs_1);
c2 = reshape(c2./nSpk,[nRO, nSpk]); %!!! Here I expect c2 similar to spokes_init

% 2st type-1 transform
c2_dcomp = c2.*dcf_all;

imgs_2 = finufft2d1(kx(:), ky(:), c2_dcomp, 1, 1e-15, nRowCol, nRowCol);
imgs_2 = imgs_2./nRowCol; %!!! Here I expect imgs_2 similar to imgs_init

% 2st type-2 transform
c3 = finufft2d2(kx(:), ky(:), -1, 1e-15, imgs_2);
c3 = reshape(c3./nSpk,[nRO, nSpk]); %!!! Here I expect c3 similar to spokes_init

% more iterations
% todo ...

[DiamonDinoia]

Hi @zongjg,

I think the difference arises from the difference between an adjoint nufft and an inverse fft. The Ajoint is not the inverse.
Here's there is an example of how to invert a nufft: https://finufft.readthedocs.io/en/latest/tutorial/inv1d2.html

To illustrate the difference, the following python applies the adjoint and plots it. I also provided Conjugate gradient inverse nufft and I plot the result of the inversion.

import numpy as np
from skimage.data import shepp_logan_phantom
from skimage.transform import resize
import finufft
import matplotlib.pyplot as plt


def nufft_inverse_cg(c1, kx_flat, ky_flat, nRowCol, nRO, nSpk,
                     eps=1e-6, tol=1e-4, max_iter=20):
    """
    Solve A x = c1 via Conjugate‑Gradient, where
      A(x) = reshape(finufft2d2(kx,ky,x)/nSpk, (nRO,nSpk))
    and c1 is the measured, density‑compensated k‑space of shape (nRO,nSpk).
    """
    M = nRO * nSpk
    N = nRowCol

    # Forward NUFFT operator
    def A(x):
        flat = finufft.nufft2d2(kx_flat, ky_flat, x, eps=eps, isign=-1)
        return flat.reshape((nRO, nSpk), order='F') / nSpk

    # Adjoint NUFFT† operator
    def AT(y):
        y_flat = y.reshape(M, order='F')
        img_flat = finufft.nufft2d1(
            kx_flat, ky_flat, y_flat, (N, N),
            eps=eps, isign=+1
        )
        return img_flat.reshape((N, N), order='F') / N

    # Initialize
    x = np.zeros((N, N), dtype=np.complex64)
    r = AT(c1 - A(x))
    p = r.copy()
    rr_old = np.vdot(r, r)  # complex, but we only use real part below
    norm_r0 = np.sqrt(rr_old)

    for k in range(1, max_iter+1):
        Ap       = AT(A(p))
        alpha    = rr_old / np.vdot(p, Ap)
        x       += alpha * p
        r       -= alpha * Ap

        rr_new   = np.vdot(r, r)
        res_norm = np.linalg.norm(r)     # now a real number
        rel_norm = (res_norm / norm_r0).real          # relative
        if rel_norm <= tol:
            print(f"Converged at iter {k}, ‖r‖₂ = {rel_norm:.2e}")
            return x, {'niters': k, 'res_norm': rel_norm}

        p      = r + (rr_new/rr_old)*p
        rr_old  = rr_new

    return x, {'niters': max_iter, 'res_norm': rel_norm}

# Params
nRO, nSpk = 512, 1000
nRowCol = nRO // 2
eps = 1e-6

# Trajectory
golden_ratio = 0.6180339887498948
golden_angle = 2*np.pi*(1-golden_ratio)
angles = golden_angle * np.arange(1, nSpk+1)
k_rad = np.linspace(-np.pi, np.pi, nRO)
kr_grid, theta_grid = np.meshgrid(k_rad, angles, indexing='xy')
kx = (kr_grid*np.cos(theta_grid)).T.astype(np.float32)
ky = (kr_grid*np.sin(theta_grid)).T.astype(np.float32)

M = nRO * nSpk
kx_flat = kx.reshape(M)   # 1-D, length M
ky_flat = ky.reshape(M)

# Density compensation
dcf = np.sqrt(kx**2 + ky**2)
dcf[dcf < 1e-5] = 1.0
dcf /= dcf.max()

# Initial image
phantom_img = resize(shepp_logan_phantom(), (nRowCol, nRowCol),
                     mode='reflect', anti_aliasing=True)
imgs_init = (phantom_img*1000).astype(np.complex64)

# ——— Type‑2 NUFFT (image → nonuniform samples) ———
# Pass imgs_init *as a 2-D array*, not flattened
spokes_flat = finufft.nufft2d2(
    kx_flat,                # length‑M
    ky_flat,                # length‑M
    imgs_init,              # shape (nRowCol, nRowCol)
    eps=eps,
    isign=-1
)
# reshape to (nRO, nSpk) in column‑major order
spokes_init = (spokes_flat / nSpk).reshape((nRO, nSpk), order='F')

# apply density compensation
c1 = spokes_init * dcf

# ——— Type‑1 NUFFT (nonuniform → image) ———
# Flatten c1 too:
c1_flat = c1.flatten(order='F')    # length M

imgs1_flat = finufft.nufft2d1(
    kx_flat,                       # float[M]
    ky_flat,                       # float[M]
    c1_flat,                       # complex[M]
    (nRowCol, nRowCol),            # output image size
    eps=eps,
    isign=+1
)
imgs_1 = imgs1_flat.reshape((nRowCol, nRowCol), order='F') / nRowCol

imgs_rec, info = nufft_inverse_cg(
    c1, kx_flat, ky_flat,
    nRowCol, nRO, nSpk,
    eps=eps, tol=1e-4, max_iter=50
)
print(f"CG converged in {info['niters']} iterations, residual norm = {info['res_norm']:.2e}")


fig, axs = plt.subplots(1, 3, figsize=(15, 5))
for ax in axs:
    ax.axis('off')

axs[0].imshow(np.abs(imgs_init))
axs[0].set_title("Initial Phantom")

axs[1].imshow(np.abs(imgs_1))
axs[1].set_title("1‑step Adjoint Reconstruction")

axs[2].imshow(np.abs(imgs_rec))
axs[2].set_title("CG‑based Reconstruction")

plt.tight_layout()
plt.show()