Help with simple gridding

[thomasediolausson]

Hi,

I'm trying to grid an image with non-Cartesian coordinates using the KbInterpAdjoint module. I was thinking that a simple test of performing an adjoint after the gridding then nufft step would give a visually similar image as to without the gridding. But I'm seeing these artifacts when my trajectory(omega) is non-Cartesian. Do you know what could be the cause of these?

I'm adding some code to reproduce...

import torch
import torchkbnufft as tkbn
import numpy as np
from skimage.data import shepp_logan_phantom

x = shepp_logan_phantom().astype(np.cfloat)
im_size = x.shape
# convert to tensor, unsqueeze batch and coil dimension
# output size: (1, 1, ny, nx)
x = torch.tensor(x).unsqueeze(0).unsqueeze(0).to(torch.complex64)

spokelength = x.shape[-1] * 2
grid_size = (spokelength, spokelength)
nspokes = 100

ga = np.deg2rad(180 / ((1 + np.sqrt(5)) / 2))
kx = np.zeros(shape=(spokelength, nspokes))
ky = np.zeros(shape=(spokelength, nspokes))
ky[:, 0] = np.linspace(-np.pi, np.pi, spokelength)
for i in range(1, nspokes):
    kx[:, i] = np.cos(ga) * kx[:, i - 1] - np.sin(ga) * ky[:, i - 1]
    ky[:, i] = np.sin(ga) * kx[:, i - 1] + np.cos(ga) * ky[:, i - 1]
    
ky = np.transpose(ky)
kx = np.transpose(kx)

ktraj = np.stack((ky.flatten(), kx.flatten()), axis=0)

ktraj = torch.tensor(ktraj)

X, Y = np.meshgrid(np.linspace(0, 2*np.pi, im_size[0]), np.linspace(0, 2*np.pi, im_size[1]))
grid = torch.tensor(np.stack((Y,X), axis=0)).to(torch.float32)

grid_ob = tkbn.KbInterpAdjoint(im_size=im_size, grid_size=im_size)

x_gridded = grid_ob(x.view(1,1,-1), omega = grid.view(2,-1))

nufft_ob = tkbn.KbNufft(im_size=im_size)
# outputs a (1, 1, klength) vector of k-space data
kdata = nufft_ob(x, ktraj)
kdata_prior_gridded = nufft_ob(x_gridded, ktraj)

nufft_adjoint_ob = tkbn.KbNufftAdjoint(im_size=im_size)
adjoint = nufft_adjoint_ob(kdata.to(torch.cdouble), omega=ktraj)
adjoint_prior_gridded = nufft_adjoint_ob(kdata_prior_gridded.to(torch.cdouble), omega=ktraj)

adjoint:

adjoint_prior_gridded:

[mmuckley]

Hello @thomasediolausson, I apologize for the late response as I was out of office and on vacation the last few weeks.

I think I'm a little bit confused. After you apply the KbInterpAdjoint, the data should be gridded and the final step should be an inverse FFT, not a NUFFT adjoint. Am I missing something here?

[mmuckley]

@thomasediolausson actually I think I'm still confused. Why would you adjoint grid directly the image data? The sequence should always be like this for forward:

x -> FFT -> Interp -> k-space data

And this for backward:

k-space data -> adjInterp -> iFFT -> x_est

But it looks like you're doing

x -> adjInterp -> FFT -> Interp -> AdjInterp -> iFFT -> x_est

which doesn't make sense to me.

[thomasediolausson]

Thank you for your reply @mmuckley.

I have deformations of my image which I want to compute non-cartesian k-space for.

My goal was to try to warp the image with your KbInterAdjoint function, is this possible?

So as a dumb test, the code is doing
x (cartesian) -> nufft -> adjnufft (top image)
x (cartesian) -> adjinterp (cartesian) -> nufft -> adjnufft (bottom image)

Besides the modulation from the kernel, I expected the results to be visually the same when the adjinterp is given cartesian coordinates. Am I misunderstanding something?

[mmuckley]

Ah I see, yes I don't think adjinterp can be used this way. I guess you have some sort of deformation field? You may want to use bilinear interpolation instead of the Kaiser-Bessel kernel. The NUFFT kernel is designed so that the image-space error is minimized with interpolations in k-space, but it's not necessarily built for deformations.

If the KB kernel is used in the image deformation literature feel free to share a paper.

[thomasediolausson]

Thank you for the fast replies @mmuckley.

I'm currently playing around with some interpolation methods, so I thought that your function would work for images too.

Some work relating deformation fields directly to k-space is here, where they use the type-3 nufft from the finufft package to go from non cartesian spatial points to non cartesian k-space points. So I'm basically trying to find differences by using a similar interpolation function in image space that finufft would make.

Do you know if the finufft uses a kaiser-bessel kernel in the type-3 forward step in image space?

[mmuckley]

Hello @thomasediolausson, I don't think they use the Kaiser-Bessel kernel. I thought they used a sequence of Gaussian approximations. They have a few papers detailing their approach and documenting its performance vs. other NUFFT packages.