Bruker MP2RAGE + CS reconstruction : MRIReco.jl and BART

[aTrotier]

Hi,

I have reimplemented my reconstruction of a custom Bruker sequence : MP2RAGE + CS
https://github.com/aTrotier/EDUC_JULIA_CS_MP2RAGE/

you can find the corresponding pluto webpage here : https://atrotier.github.io/EDUC_JULIA_CS_MP2RAGE/

The acq,raw object are created with a function in utils_MP2RAGE.jl

Remarks

• In order to use FISTA I had to use the latest commit of RegularizedLeastSquares.jl (Fista bugfix)

• I was not able to get the same results as with BART with the FISTA implementation when I try to overregularized it will just put a lot of 0 in the image.

• but ADMM the results are really close

• In terms of reconstruction times, MRIReco is slower by a factor of at least 5 for the FISTA. (ADMM is really slow > 500 sec) maybe just switching the acq data to ComplexF32 will help.

Let me know if you have some thoughts about how I worked with MRIReco.jl for custom sequences

[tknopp]

Its really cool that you are pushing MRIReco forward. I cannot deeply work into your code right now but we should work on integrating as much as possible of your custom code to make this working with MRIReco out of the box.

Regarding speed: If we there are no big differences on the algorithmical level (like for ESPIRIT) than this is because of suboptimal implementation. With the Julia profiler (-> ProfileView) one can find the bottlenecks and fix them. This is how I speed up NFFT.jl. But this is tedious work and I right now don't have enough time to fix them alone.

[aTrotier]

Let's discuss what I can try to incorporate in MRIReco.

Reuse parts of acqDataFromMP2RAGE(b::BrukerFile)

In terms of custom functions the main one is function acqDataFromMP2RAGE(b::BrukerFile). Inside it I need to create the profiles which is very specific to my custom sequence.

Maybe the part that fills params can be shared between all the sequences ? I had to change some fields like : params["TI"] = parse.(Float32,b["EffectiveTE"]) but I think 90% will remains unchanged in most of the sequence.
I can implement something like : setRawParams(b::BrukerFile) ?

Extract Calibration data

I had to create a function to extract the 3D k-space then crop the center to get the calibration data.

The crop function already exists. I can enhanced it to 3D if you want.

I did not find a function that convert an acquisitionData to kspace -> extract3DKSpace(acqData::AcquisitionData)

I found 2 functions that does something similar / closely related :

• setcalibrationarea(A::Array{Int64,2},...)
• espirit(acqData::AcquisitionData, ...) is doing that in 2D before calling espirit(calibData::Array{ComplexF64,3}, ...)

Playing with BART

I can create a function that directly create a CFLfile from acquisitionData, it might be something like a a loop on slice/echo/NR with extract3DKSpace and then reshape/permute the data according to BART cfl format ?

convertToBart(acqData::AcquisitionData ;coils = :, echo = :,...)

Additional stuff

• Removing conversion from ComplexF32 to F64 in acquisitionData
• Lambda only accept Float -> need to cast Int to Float

[aTrotier]

When BartIO.jl will be more stable than my forked version, do you want to add more benchmarks vs BART ?

It won't really be possible to do something as clean as in NFFT but it might be interesting to see where we stand in term of performance.

[tknopp]

yes absolutely, actually we even have an independent repository for that: https://github.com/MagneticResonanceImaging/MRIRecoBenchmarks
The idea is to have one folder per benchmark.

[aTrotier]

I played a little bit more with ADMM and get comparable results with BART in 110sec.
I guess reducing the inner tolerance helps :)

I have to check the parameter in bart to do something close in terms of benchmark

[tknopp]

I does not matter too much if MRIReco.jl is slower at the moment. A proper benchmark suite will enable us to work on performance issues.

[aTrotier]

Sorry, I was thinking in term of parameters for the ADMM algorithm :)

[aTrotier]

With BartIO we might get rid of the matlab part :) Do you prefer to work on real 3D datasets for benchmark or a numerical phantom using the BART function ? Le lun. 7 févr. 2022 à 11:58, Tobias Knopp ***@***.***> a écrit :

…

yes absolutely, actually we even have an independent repository for that: https://github.com/MagneticResonanceImaging/MRIRecoBenchmarks The idea is to have one folder per benchmark. — Reply to this email directly, view it on GitHub <#65 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AC5P7O6WGIAF5CKUMWGBZ4DUZ6QW7ANCNFSM5NSJ4AKQ> . You are receiving this because you authored the thread.Message ID: <MagneticResonanceImaging/MRIReco.jl/repo-discussions/65/comments/2124411@ github.com>

-- ------------------------------------------------------------ TROTIER Aurélien 31 allée de la Biotte 33470 Gujan-Mestras 06.09.12.61.61 ***@***.***

[tknopp]

Difficult question. I always like real data. But in this case we could also look at deviation and calculate errors if we go for phantom data. So slight preference on phantom data.

By the way: the "Benchmark" in MRIRecoBenchmarks actually means comparing speed and accuracy.

[aTrotier]

I am starting another project where I will use BART.

I again need to extract the k-space.

I previously wrote this function :

function extract3DKSpace(acqData::AcquisitionData)
    if !MRIReco.isCartesian(trajectory(acqData, 1))
        @error "extract does not yet support non-cartesian sampling"
    end
    nx, ny, nz = acqData.encodingSize[1:3]
    numChan, numSl = MRIReco.numChannels(acqData), MRIReco.numSlices(acqData)
    numEcho = length(acqData.traj)
    kdata = zeros(ComplexF64, nx * ny * nz, numEcho,numChan)

    for echo = 1:numEcho
        for coil = 1:numChan
            kdata[acqData.subsampleIndices[numEcho], echo,coil] .= kData(acqData, echo, coil, 1)
        end
    end
    kdata = reshape(kdata, nx, ny, nz, numEcho, numChan)
    return kdata
end

@tknopp I can add the repetition to the function and add it to MRIReco ?