WIP : Implement labels

[aTrotier]

This PR will read the labels section of PR and in the end fill each profile with the corresponding MDH.

It could be used to reconstruct multi-contrast / slice / rep and also to use the MDH for cartesian position of the data in kspace.

Related to : #308

Remaining things to do

• Manage extension in sequence structure
• Link ADC with the MDH
• fill the ismrmrd MDH
• Implement a way to create label in the seq designer

[aTrotier]

@cncastillo can you take a look at the initial implementation before going further.

[aTrotier]

@cncastillo forget to push the rest !

[cncastillo]

Sure! Thanks for implementing this, I will take a look later today.

[cncastillo]

Hi! This looks awesome. I just have a quick nitpick, I added comments to some code blocks that could be benefited by the use of multiple dispatch. I don't have a clear idea yet, so it is only a comment. But feel free to continue developing this.

[aTrotier]

Labels are now implemented

The label field is first empty and filled on request (if not empty). The ADC is not directly linked but the user can use the index of each adc block and look at seq.label[i].

[aTrotier]

I am able to use the label for multi-slice reconstruction with KomaMRI but we are limited by something in MRIReco.jl :
The multi-slice reconstruction from MRIReco expects that all slices has the same trajectory (when I want to have different trajectory I manually label the raw.profiles along the rep or contrast dimensions which handles different trajectory).

In my local example the trajectory is the same so I just "cut" the trajectory at 1/3 of the total traj

acq = AcquisitionData(raw)
acq.traj[1].nodes = acq.traj[1].nodes[:,1:4096]

params = Dict{Symbol,Any}()
params[:reco] = "direct"
params[:encodedSize] = acq.encodingSize[1:2]
params[:densityWeighting] = true

im_reco = reconstruction(acq,params)

Reconstruction with cartesian labels

I was able to reconstruct the same dataset with the FFT (using the cartesian label LIN and PAR) :

It requires the following modifications to the reconstruction pipeline :

• encodedSize was estimated by Koma = 63,63 -> I changed that with Nx,Ny from seq.DEF
• the subsampleIndices are wrong because the UInt16(0), #center_sample uint16: Sample at the center of k-space is equal to 0.
• With the flag estimateProfileCenter = true it works (for centered acquisition). I don't know if we can handle that with the calculated ktraj

raw2 = deepcopy(raw)
raw2.params["trajectory"] = "cartesian"
raw2.params["encodedSize"] = [seq.DEF["Nx"],seq.DEF["Ny"]]
#raw2.params["reconSize"] = [seq.DEF["Nx"],seq.DEF["Ny"]]

acq2 = AcquisitionData(raw2,estimateProfileCenter=true)

params = Dict{Symbol,Any}()
params[:reco] = "direct"

im_reco_2 = reconstruction(acq2,params)

The reconstruction part might cause some issues between labeled "cartesian" sequence etc @cncastillo

[aTrotier]

I removed the label from Sequence struct. It is easier to handle that way and is not directly connected to the sequence (issue with subset of sequence, or if the user already extract the label and continue to add block etc)

now the label can be extracted with : get_label(seq) and returns a Vector{AdcLabels} (and maximum(label))

[aTrotier]

I now estimates the center directly from ktraj it works for both cases with my simple example (I need a partial-fourier acquisition along the readout in order to see if it really works)

2 points remains to be added in the reconstruction part :
raw2.params["trajectory"] = "cartesian" -> nothing to do here except adding that in the documentation

raw2.params["encodedSize"] = [seq.DEF["Nx"],seq.DEF["Ny"]] -> The Nx/Ny from seq.DEF (calculated during the conversion from Pulseq) are different to the one calculated here in ISMRMRD.jl :

    if haskey(seq.DEF, "FOV")
        FOVx, FOVy, _ = seq.DEF["FOV"] #[m]
        if FOVx > 1 FOVx *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm
        if FOVy > 1 FOVy *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm
        Nx = round(Int64, FOVx / Δx[1])
        Ny = round(Int64, FOVy / Δx[2])
    else
        FOVx = Nx * Δx[1]
        FOVy = Ny * Δx[2]
    end

64,64 vs 63,63 (from the FOV). Should I put a ceil instead of round ? @cncastillo

BTW : I opened an issue to implement all the dimensions independently (SLC, ECO...) in MRIReco : MagneticResonanceImaging/MRIReco.jl#246

[cncastillo]

Oh sorry for the late reply.

64,64 vs 63,63 (from the FOV). Should I put a ceil instead of round ? @cncastillo

I don't have a strong preference, ceil could make more sense.

So can this be used in the UI already? Or is there something that needs to be modified for the reconstruction to incorporate the labels?

[aTrotier]

My example with 3 slices with SLC / LIN (/ PAR) set correctly can be reconstructed by KomaUI() with obj_ui[] = brain_phantom3D() and I can slide along the 3 slices in the UI with and without defining the LIN (SLC are still defined)

You can give it a try with the attached sequences :
Archive.zip

I am pretty sure it will failed at some point if the SLC label are not defined or if it is a 3D acquisition without label.
We need to be careful and use more test sequences.

The estimated matrix is still 63 (estimated from the FOV and k-space trajectory) it is not an issue for the reconstruction and not related to my modifications

---

MRIReco only handles SLC / ECO / REP as label and only the eco can have a different trajectory (I think it is one of my PR that added this enhancement).
Maybe a google GSoC could be to enhance MRIReco with all the label that would be directly beneficial for Koma -> see issue
Seems it is too late for that :)

[aTrotier]

I am not a user of the seq designer, do you have an idea about what should be done to implement the labels ?

[aTrotier]

The seq designer is already compatible with label:

lSet = LabelSet(1,"LIN")

seq = Sequence()  # empty sequence
tmp_seq = Sequence([gr],[exc])       # adding RF-only block
tmp_seq.EXT = [[lSet]]
seq += tmp_seq

What should I add to the documentation ?

• Explanations/Sequence Events
• I need to describe somewhere the function 'get_label(seq)'
• Should we add an example to show how to use the label for reconstruction ? Like removing ADC lines with labels -> ECO = 1 or something like that

After the doc enhancement, the remaining issue is mostly related to MRIReco limitations.

It would be great if we could merge this feature before summer (I would like to start writing my MRI course with this possibility)

[cncastillo]

Hi @aTrotier Thanks for the enhancements.

Should we add an example to show how to use the label for reconstruction ? Like removing ADC lines with labels -> ECO = 1 or something like that

I think a literate tutorial would be great.

It would be great if we could merge this feature before summer (I would like to start writing my MRI course with this possibility)

Awesome! I have no problems merging this before summer. @Stockless can you make sure this is merged asap and tag new versions so it is usable by summer?

[aTrotier]

I'll wrote the literate example next week.
I think I also need to add more test to the sequence construction.

Should I wrote the example using a dummy .seq séquence or with the koma sequence designer ? (And just mention that it works with pulseq labels)

[JanWP]

Hello everyone,

Just stopping by to say a big thanks for this PR !
It's very timely: For the past few days I have been trying to implement an AFI sequence, which simply consists of two interleaved GRE acquisitions with different TRs. The k-space lines from the two images are interleaved, just like for multi-echo data, and I was struggling with the reconstruction.

With this PR, it worked perfectly. I have to use the ECO label for now (I think it's preferable to use SET, and that's what I tried initially, but it's not implemented on the reconstruction side yet). There was only one tiny issue: I had to change the "EXT.jl" to "Ext.jl" here.