espirit for 3D data

[aTrotier]

If I wan to use espirit for 3D acquisition
I have to first extract the calibration data in an Array of size : kx × ky [× kz] × coils

I guess I can first extract the k-space with kdata(acquisitionData) but for what I see the results is not "ordered".
Am I missing something ?

Should I wrote a function for that ?

[aTrotier]

Ok, I have something working

kspace = extract3DKSpace(acq)

calibSize = parse.(Int,b["CenterMaskSize"]);
calibData = crop(kspace[:,:,:,2,:],(calibSize,calibSize,calibSize));


"""
Documentation extractKSpace
"""
function extract3DKSpace(acqData::AcquisitionData)
    if !MRIReco.isCartesian(trajectory(acqData, 1))
        @error "espirit 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], numEcho,coil] .= kData(acqData, echo, coil, 1)
        end
    end

    return reshape(kdata, nx, ny, nz, numEcho, numChan)
end

"""
Documentation : Crop the central area
"""
function crop(A::Array{T,4}, s::NTuple{3,Int64}) where {T}
    nx, ny, nz = size(A)
    idx_x = div(nx, 2)-div(s[1], 2)+1:div(nx, 2)-div(s[1], 2)+s[1]
    idx_y = div(ny, 2)-div(s[2], 2)+1:div(ny, 2)-div(s[2], 2)+s[2]
    idx_z = div(nz, 2)-div(s[3], 2)+1:div(nz, 2)-div(s[3], 2)+s[3]
    return A[idx_x, idx_y, idx_z,:]
  end

For now the espirit is taking too long to compute... I had to kill the process. I don't know why
sens_spirit = espirit(calibData, (size(kspace,1),size(kspace,2),size(kspace,3)),(6,6,6))

[aTrotier]

Output size of sens maps : 20^3 and 64^3 works

Julia ram takes 19 go when trying to do 128x128x64.

I am working on my mac and not on my lab server. Maybe I'm using some kind of swap memory

[aTrotier]

Ok I was able to get the maps It is just really long.
Pretty sure it is a swap memory problem :)

[aTrotier]

Results versus BART:

But BART takes 1 second to compute the maps.
My guess is that the C implementation in BART is far more efficient than the one in the initial Matlab package. @JakobAsslaender

[JakobAsslaender]

Is extract3DKSpace a new function you wrote? My apologies, I never used the extract functions myself and, therefore, did not convert them to ND. Quick comment w/o having looked at the code in too much detail: I generally tried to write one set of code for ND rather than duplicating code with 2D and 3D versions. For this, CartesianIndices are tremendously helpful. You can check out some of the ESPIRiT functions for examples. Though, here it seems that all it takes is to replace nx, ny, nz = acqData.encodingSize[1:3] by ima_shape = acqData.encodingSize[1:end-1] or something similar. and

Another side comment, as this thread is about speed: With the line zeros(ComplexF64, nx * ny * nz, numEcho,numChan), you are filling the entire array with zeros, only to overwrite all of it in the for loop. The call Array{ComplexF64}(nx * ny * nz, numEcho,numChan) is much faster here, as it does not write anything in the arrays.

[aTrotier]

Is extract3DKSpace a new function you wrote?

Yeah just wanted something to both extract the k-space (in order to write the cfl and I also cropped the kspace to get the calib region.
For now I want to have something working, later on for sure I should just change some parts of espirit to directly handle 3D :)

[JakobAsslaender]

Hi,
Are you using the master branch or the latest release? If the latter, can you try it with the master branch? I made some major modifications that save quite a bit of memory and speed up the calculation by using power iterations instead of an SVD. Here some comments on the state of things:

• With version on the master branch, calculating 3D maps on a desktop/laptop should be ok: it has to allocate a complex valued array of the size Nx x Ny x Nz x Ncoils^2, so if you, e.g. use ComplexF64 with above matrix size and, say, 12 coils, you will need round about 2.5GB, plus a little extra. Can you check if this roughly the case?
• Did you make sure that you have multiple threads enabled for Julia? The most essential part of the code is multi threaded. You can check it with the call Threads.nthreads() and you start Julia multithreaded with julia --threads 4 (replace the number with the number of cores on your machine). But there is also the issue of Julia vs. BLAS threads. I implemented some switching between them, which should work fine, but it is not the most beautiful implementation in history.
• The latest version is fastest if you a) calculate only one set coil maps by setting nmaps = 1 and b) by using this flag: use_poweriterations = true. Both should be the default.
• Can you run a quick benchmark comparing the Julia implementation and BART? It would be helpful to move forward from here. Make sure the settings are the same. I think BART uses by default a calibration size of 24x24x24 and a kernel size of 6x6x6. Beyond that, I am not surprised that BART is faster than our implementation, it is written by the inventor of ESPIRiT, who is also a pretty good programmer ;). Yet I am surprised BART is THAT quick for a 3D map, in my use case (192^3 to 256^3 and 20 coils) on 40 cores, it takes a few minutes. And yes, the implementation in BART it is much different compared to their MATLAB implementation. I think they are approximating all their SVDs. I have done so for the last, per voxel one. But I think they are also approximating the first big SVD with some kind of power iterations. But as we need more than one eigenvector from this SVD, it is a bit more involved. I don't think I will have time to look into the BART implementation any time soon, but feel free to check it out and make PRs. We can also talk to @uecker and see if he is willing to support us in building up competition ;).

[aTrotier]

With master branch and a matrix of 128x128x64x4 + 1 thread + F32

 60.602676 seconds (5.24 M allocations: 1.073 GiB, 0.84% gc time)

julia> Threads.nthreads()
1

Indeed it was an issue with my memory

[uecker]

@JakobAsslaender If you mention me, you are already talking to me. We only use power iteration for the second step and not the first. But we use single precision as this is generally sufficient for MRI. This alone would make a huge difference. And there are some other steps one can do to make it faster which may not be implemented in the Matlab code. I also had discussions with @spinicist when he implemented it in C++, maybe he remembers exactly what those steps were.

[JakobAsslaender]

Haha, fair point, this my brain changing its mind while typing... Thanks for the feedback!

Ah, ok. So you are using a regular SVD, or an eigen decomposition for the first step? I was trying to understand the debug output of BART and got the impression that the number of singular/eigen values above the first threshold were much smaller compare to the Julia implementation. But when I compared the Julia implementation to your Matlab code, they matched precisely.

I was aiming to write precision-agnostic code in Julia, but it might be worth double (pun intended) checking if that is consistent throughout, @aTrotier . If everything works as intended, you can just call the function with ComplexF32 data and all operations will be preformed with single precision.

[aTrotier]

I have to check at what point my data are converted to F64 because from Bruker I get F32 :)

[aTrotier]

Just for reminder :
espirit can accept and return F32 data. However the sensitivity maps needs to be in F64 for reconstruction().

[uecker]

We use an eigen decomposition in the first step which reduces problem size a lot, but then there are also some other tricks... One thing none needs to watch out for with single precision is that one does not exceed the range due bad scaling of the data.

[JakobAsslaender]

Ah, great! We should switch to an eigendecomposition as well. Can you quickly tell me if the first threshold acts on the eigenvalues, or on their square root? If I am assuming that you are performing the eigendecomposition AᴴA = EΛEᴴ where A is system matrix as defined in your paper; this implies the relation Λᵢᵢ = Sᵢᵢ² (and E=V) given the SVD A = USVᴴ. Just asking because this might explain the discrepancy between the number of vectors kept.

[uecker]

The threshold acts on the eigenvalues (relative to the largest), i.e. on the square of the singular value of the original matrix. I think the debug output is the singular values. I just remembered the other thing we do which has a much bigger impact: We do not zeropad the kernels to full size before transforming them to the image domain. Instead we zeropad to only the double size, then FFT and multiply with the conjugate and only then only transform this much smaller matrix to full size using FFTs and zero padding. Left as an exercise to the reader why it works ;-)

[JakobAsslaender]

Ok, this sounds like black magic! I think this might indeed address the bottleneck of our current implementation.

But you are still computing the final eigendecomposition/power iterations on the full size image?

[uecker]

Mathematics == Magics

Yes, the final eigendecomposition is done on the full size (but if you are really limited in memory, there is another trick you can do...)

[JakobAsslaender]

Mathematics == Magic
or
Mathematics === Magic ?

[spinicist]

Hello,

Thanks @uecker for the introduction - @JakobAsslaender we haven't met, but literally today I was planning to e-mail you about how to get Proton Density estimates from finger-printing type reconstructions. But this isn't that place for that conversation.

As Martin mentioned, I have a C++ implementation of ESPIRiT here https://github.com/spinicist/riesling/blob/main/src/espirit.cpp. I think Martin has answered most of your questions already, and frankly it sounds like my implementation needs a good deal of work. The one thing that I think hasn't been discussed here explicitly is that, at Martin's suggestion, I do the upsampling and image-space PCA/eigendecomposition step slice-by-slice. I'm dealing with oversampled 3D acquisitions and doing the image-space analysis on fully upsampled images was very impractical. Slice-by-slice made it feasible, but my implementation is still slow.

[tknopp]

Hi @JakobAsslaender,

I had a look at optimizations of espirit and tested some things where it would be good if we could move this forward together. My work can be found here: https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/tk/espiritSpeedup/src/Tools/CoilSensitivity.jl
The idea is to implement Martin's suggestion #55 (comment).

Remarks / questions:

• In line https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/tk/espiritSpeedup/src/Tools/CoilSensitivity.jl#L175 I create the padded ksize kernel.
• It is in turn used in kern2, kernel_k and kernel_i: https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/tk/espiritSpeedup/src/Tools/CoilSensitivity.jl#L193
• What I absolutely do not understand is your line https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/master/src/Tools/CoilSensitivity.jl#L199. If I get it right, it only fills the "lower positive frequencies" while the "lower negative frequencies", which are in the second part of kernel_k are not filled. I tried to fix this here: https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/tk/espiritSpeedup/src/Tools/CoilSensitivity.jl#L202 but maybe I just do not understand your code and you apply the shift at some other part of the code.
• The next question is when to upsample. I tried to do it here: https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/tk/espiritSpeedup/src/Tools/CoilSensitivity.jl#L219, which is early on and gives a small speedup. Or we can do it later (https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/tk/espiritSpeedup/src/Tools/CoilSensitivity.jl#L253) on the eigenVecsPadded and eigenValsPadded.
• In my code I currently do the upsampling in image domain using imresize which is not what Martin originally suggested. I just did this because struggling with the fftshift issue discussed above. We can of course change this back.

This is just some initial brain dump. I have to say that I am not deeply in that paper and basically was just looking at the code.

[spinicist]

Hello,

I had the same question for Martin when I wrote my implementation. Sadly, you cannot upsample at the end once you have the eigenvalues and vectors. The problem is upsampling regions of fast varying phase is not well defined, and you likely end up introducing phase poles into your sensitivity maps.

Consider what happens if you have two neighbouring voxels in the low resolution eigenvector maps with values -1 +1. If you upsample to twice the resolution, then these values will be linearly interpolated and you get -1 0 +1 - i.e., you have introduced a voxel with zero magnitude in your sensitivity maps. Not good!

If instead you do the upsampling before the eigenanalysis, what happens in practice is the same voxels will either end up as -1 +i +1 or -1 -i +1, i.e. the magnitude will be constant but the phase varies. This is physically meaningful.

Maybe there is some clever way of upsampling magnitude and phase separately, instead of real and imaginary? I'm not aware of one.

[tknopp]

Thanks for the explanation. I was just curious to hear if there are more opportunities.

[spinicist]

I am definitely curious about that too!

[JakobAsslaender]

That is a great explanation, thanks!

[tknopp]

I made a quick test to apply the upsampling even later at this stage https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/master/src/Tools/CoilSensitivity.jl#L133
That is after the filtering step.
Without having a look at the image at least the error looks not completely bad. It is 0.045 instead of 0.004, which we get when doing the upsampling earlier.

[tknopp]

Crossref MagneticResonanceImaging/MRIRecoBenchmarks#6 for the related benchmark.