Benchmark vs Bart

[aTrotier]

Hi,

I have created a repository that perform a benchmark against Bart during in a Github actions : https://github.com/aTrotier/MRIRecoVsBART_Benchmark

It is done in a pluto notebook that is published : https://atrotier.github.io/MRIRecoVsBART_Benchmark/test_bart.html

The notebook test_bart.jl perform a CS reconstruction of a shepp logan phantom with Fista.

Results are :

Reconstruction time julia = 256.0 sec vs t_bart = 19.432334866 sec

Bart is 13.19 faster

Precision :
RMSE:

julia = 0.0182

bart = 0.0167

SSIM:

julia = 0.875

bart = 0.962

I can't find the right parameters to get a accuracy close to BART. (It is also something I pointed out here : https://atrotier.github.io/EDUC_JULIA_CS_MP2RAGE/)
We can create other notebooks that will be published in order to benchmark others functions, algorithm.

@tknopp we can merge it into your Benchmark repository at some points

[tknopp]

Thanks, this is very useful. Basically we seem to have two problems here, the speed and the accuracy.
One question on the speed. Could it be that BART runs multi-threaded? If you look at this line: https://github.com/MagneticResonanceImaging/MRIRecoBenchmarks/blob/master/benchmark1/runBenchmark.sh#L7
we did explicitly set the number of cores to make things comparable.

[aTrotier]

Reconstruction is now performed in F32 which gives :

Reconstruction time julia = 69.0 sec vs t_bart = 22.506465278 sec
number of threads for julia = 2 and bart = 2
Bart is 3.07 faster

RMSE:
julia = 0.0174f0
bart = 0.0167

SSIM:
julia = 0.888
bart = 0.967

[tknopp]

Julia is likely not using the 2 threads right now. I suggest to set ENV["OMP_NUM_THREADS"]=1 and retry.
@migrosser and I have, however, found additional performance issues in FFTOp wich we have a look at.

[aTrotier]

Now we're close during CI. On my mac, I bart is faster between 2 to 3 times.

Reconstruction time julia = 60.0 sec vs t_bart = 47.455016155 sec
number of threads for julia = 2 and bart = 1
Bart is 1.26 faster

Should I run the the reco twice to remove the compilation ?

[migrosser]

Concerning the accuracy, the proper choice of the regularization factor also seems to be an issue because bart performs a different form of data rescaling, which step. When using your benchmark script I obtained the best results using λ=3.e-2 and normalizeReg=true. Finally, it seems that FISTA did not fully converge after 30 iterations. When increasing the number of iterations to 40 or 50. The images improved further. I am not yet sure where exactly the difference to bart comes from. First i thought that it might be the choice of ρ but this seems to be consistent with bart.

[aTrotier]

rescaling -> related to #92

My experience with Bart make me think that scaling the data before is a good idea to easily choose the regularization parameter (always close to 0.01 :D ). Do you think it is a good idea to implement that ?

Maybe for a first test I can scale the data by hand accordingly to bart output : Inverse scaling of the data: 193895.437500

An other possibility is the wavelet implementation.
I asked @uecker about the wavelet detph level for non-squared matrix :

Hi Aurélien,
Am Dienstag, den 14.06.2022, 12:42 +0200 schrieb Aurélien:

Hi,
I have a question about the implementation of the DWT for CS.

When performing the DWT on a non-squared matrix. Are you zero-padding the
image to the next 2^n size in order to maximize the level of wavelet
transforms ? Or just taking the smallest possible level among all the
dimensions ? (for example with a matrix 220x160 the DWT level is only 2 for
the dimension 220).

neither, it will do more levels for the larger dimensions.

Martin

In our case 128x128x128 it is not an issue but we should keep that in mind for benchmark on real datasets like the one in example.jl (matrix 220x160)

[tknopp]

I have now integrated your example in MRIRecoBenchmarks @aTrotier. Thanks for preparing!

https://github.com/MagneticResonanceImaging/MRIRecoBenchmarks/tree/master/benchmark2

Our timings are actually quite ok:

We will have a deeper look at the accuracy differences once we find some time. For now it was more important for me to improve the infrastructure of MRIRecoBenchmarks. It should now be very easy to integrate new benchmarks since benchmark1 and benchmark2 are good templates.

[aTrotier]

I did not expexted MRIReco to be as fast as bart. That's great :)
I guess the benchmark has been done with your last commits on master branch ? Should I change my CI benchmark to use the same number of threads for both ?

Maybe we should keep an historic about thus timing somewhere for tagged release of MRIReco ?

Something interesting will be to also benchmark the memory consumption, we are working with multi-echo/repetition/channels data non-cartesian acquisition in my labs and the memory is one of the limitation. However I have no clue how it can be benchmark... we can add that for a futur online meeting discussion.

@mrikasper I don't know how it is implemented in MRIreco but integrating that in Bart will required too much work I guess. Another alternative might be to implement that in sigpy but I don't really know how fast the toolbox is.

[uecker]

It should be fairly simple to add off-resonance correction to BART. It is on my TODO list. If someone needs it, I think I could add a simple time-segmented nuFFT fairly quickly.

[tknopp]

I did not want to push you there Martin. Our feature set is actually much(!) smaller than BART's. In that particular case we will start benchmarking against ourselves. What's also import that a "benchmark" is meant to be not just timing but also accuracy. And regarding accuracy it's of course simple to compare to a ground truth.

By the way: I am very excited that it is possible to mix BART and MRIReco.jl so nicely because of BartIO.jl. This could be a very valuable combination for users.

[tknopp]

@mrikasper: Do you combine (in your research) off-resonance correction with parallel imaging or not?

@aTrotier: Memory consumption is a tricky one. I failed at that when I tried to do that with NFFT.jl. But I am open to suggestions. Things may change with Julia 1.8 since it includes a memory profiler.

[mrikasper]

Yes, it's (nearly) all undersampled spirals with long readouts (10s of ms). So off resonance correction plus parallel imaging.

Like this publication: https://doi.org/10.1016/j.neuroimage.2021.118738

And recon with MRIReco.jl:
https://github.com/mrikasper/julia-recon-advances-in-spiral-fmri