Reconstructing an Undersampled Spiral

[nathanjbeaumont]

Hello everyone,
I have an undersampled spiral sequence that I want to reconstruct but I don't know how to go about it. The simulation signal looks as expected:

But reconstruction is not handling the under sampling:

I'm expecting to see an image of similar quality to figure 7 here https://onlinelibrary.wiley.com/doi/10.1002/mrm.25559:

What parameters do I need to set in the reconstruction function to do this?

Here's the code:

using KomaMRICore, KomaMRIPlots, KomaMRI, CUDA
obj = brain_phantom2D()
plot_phantom_map(obj, :ρ)
seq = read_seq("spiralNormal.seq")
plot_seq(seq)
sys = Scanner() # default hardware definition
raw = simulate(obj, seq, sys)
p4 = plot_signal(raw)
acq = AcquisitionData(raw)
reconParams = Dict{Symbol,Any}(:reconSize=>(255, 255))
image = reconstruction(acq, reconParams)
slice_abs = abs.(image[:, :, 1])
p5 = plot_image(slice_abs)

And here's the .seq:
spiralNormal.seq.txt

[cncastillo]

Hi! I believe the k-space FOV is too large for the sampling of the phantom, so you are getting aliasing artifacts. This happens because the spins are in a uniform rectangular grid, so they generate a repeating spectrum (separated in 1/dx). You can see a spurious echo at the end! (red arrow)

To increase the phantom sampling you can use the us keyword:

obj = brain_phantom2D(; us=2) # Up-sample by a factor of 2

This produces:

MRIReco, used for the reconstruction, uses a DCF based on the k-space samples, which could look slightly different from what you expect. If so, varying the DCF or reconstructing externally with the output raw data MRD would be advised.

BTW, KomaMRI contains KomaMRICore, KomaMRIFiles, and KomaMRIPlots. So using KomaMRI is the same as doing using KomaMRICore, KomaMRIFiles, KomaMRIPlots, the latter does not load the UI and recon code.

[nathanjbeaumont]

Thanks that looks more as expected.
I'm trying to recreate figure 7 here: https://onlinelibrary.wiley.com/doi/10.1002/mrm.25559

The authors say
"All acquired spiral data were reconstructed using nonuniform fast Fourier transform with a separately measured spiral trajectory that would correct the gradient imperfection, such as from eddy current or gradient delays. Images from individual coils were combined using the adaptive coil combination method, and coil-combined images were normalized to coil sensitivity map to mitigate the image intensity variation due to the coil sensitivity"

Does KomaMRI have an option to do NUFFT with a DCF?

[cncastillo]

This is more of an MRIReco thing. I looked into it, and MRIReco already uses a NUFFT with a DCF if the trajectory is not cartesian. I think the only "problems" were:

• There could be a DC offset, as the DCF is not 0 for k=0 (Based on https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.1910380117)

using MRIReco
density = samplingDensity(acq, reconParams[:reconSize])
density[1][1] = 0.0 # Fixes problem 1

• I believe the raw data is not pre-multiplied by sqrt(dcf) only the operators (the initial data b needs to be multiplied by it based on https://onlinelibrary.wiley.com/doi/epdf/10.1002/mrm.1241)

raw2 = copy(raw)
raw2.profiles[1].data .*= density[1] # Fixes problem 2

• The generated DCF is a little "noisy", I am not sure if a custom DCF can be used. Maybe @tknopp or @nHackel can help?

DCF from MRIReco.samplingDensity after modifications mentioned above:

Result (it should look sharper):

Full code:

using KomaMRI, CUDA
obj = brain_phantom2D(; us=2)
seq = read_seq("spiralNormal.seq")
sys = Scanner() # default hardware definition
raw = simulate(obj, seq, sys)
acq = AcquisitionData(raw)
acq.traj[1].circular = false #This is to remove the circular mask
reconParams = Dict{Symbol,Any}(:reconSize=>(255, 255))
image = reconstruction(acq, reconParams)
slice_abs = abs.(image[:, :, 1])
p5 = plot_image(slice_abs)

# DCF modifications
using MRIReco
density = samplingDensity(acq, reconParams[:reconSize]) # Gives sqrt(dcf)
density[1][1] = 0.0 # Fixes problem 1
raw2 = copy(raw)
raw2.profiles[1].data .*= density[1]  # Fixes problem 2
acq2 = AcquisitionData(raw2)
acq2.traj[1].circular = false #This is to remove the circular mask
image2 = reconstruction(acq2, reconParams)
slice_abs2 = abs.(image2[:, :, 1])
p6 = plot_image(slice_abs2)

[aTrotier]

@cncastillo the "noise" in the dcf is "normal" and won't impact a lot the image reconstruction.

The results seems ok to me, the spiral acquisition is especially subsampled (approximatively 256/8 turns = 16 times ?) I guess it is done like that for a MRF like acquisition.

@nathanjbeaumont the image on the siemens scanner was reconstructed with an Ice pipeline ?

[nathanjbeaumont]

I don't know what an Ice pipeline is, but I looked into it and the Siemens reconstruction was actually done using a window of the first 24 interleaves of a 48 interleaf sequence, so actually disregard the phantom image. However I do think a single interleaf simulation should be able to produce an image like this:

[cncastillo]

I have used a similar recon for our MRF implementation and have gotten images like the ones you show, but the contrast changes through the MRF sequence. These are in vivo acquisitions, btw.

MRF.mp4

For the first spiral readout, the image has a higher signal intensity, I would add a few more readouts. To reconstruct more images change Nz to the number of spiral readouts in the sequence definition.

Having said that, the reconstruction differences shouldn't affect the MRF quantification after dictionary matching as they are consistently different.

[cncastillo]

I forgot to say, but you can add readouts by scaling and rotating the single-shot acquisition.

seq = read_seq(...) # single-shot undersampled spiral
seq_mrf = Sequence()
for i in 1:Nreadouts
    seq_mrf += (rf_scaling[i] * cis(rf_phase[i])) * (rotz(spiral_angle[i]) * seq)
end

Gradients are rotated with a rotation matrix and RFs are scaled by a complex number.

[nathanjbeaumont]

Is each frame in the video you posted reconstruction of a single interleaf of a spiral? So the image quality varies drastically between TR indices? In this case figure 7 here is a cherry picked interleaf with one of the best looking reconstructions?

[cncastillo]

It could be a slightly different reconstruction, but yeah, the reconstructed interleafs can look quite different. The selected TR is probably the one that looked the best.