Draft for discussion : Subspace implementation

[aTrotier]

Hi,

I am working on a subspace reconstruction for T2 mapping Multi-Echo Spin-Echo sequence. I am currently doing the reconstruction part with BART.

This is a draft to implement the reconstruction in Julia. It seems to work (at least on fully datasets) :

• for multi-coil multi-echo -> same fix than in export function and correct SparseOp for echoes #146
• for now the reconstruction is in a specific reconstruction pipeline params[:reco] = "multiCoilMultiEchoSubspace". Maybe we can directly put it in the multi-echoes operators ?

Let me know what you think I should change ?

I think @JakobAsslaender and @alexjaffray are also interested in / using a similar implementation

example :

b = BrukerFile("MESE_FULLY")
raw = RawAcquisitionData_MESE(b)
acq = AcquisitionData(raw,OffsetBruker=true)

sens = espirit(acq)


## generate basis

function createExpBasis(TE_vec::AbstractVector{T},T2_vec::AbstractVector{T}) where {T<:Real}
    nTE = length(TE_vec)
    nSimu = length(T2_vec)
    expSignal = zeros(T,nSimu,nTE)

    for (i,T2) in enumerate(T2_vec)
        expSignal[i,:] = exp.(-TE_vec/T2_vec[i])
    end

    return expSignal
end

T = Float32
TE_vec = 7:7:50*7
T2_vec = 1:1:2000

sDict = createExpBasis(T.(TE_vec),T.(T2_vec))
lines(TE_vec,sDict[200,:])

svd_obj = svd(sDict)
V = svd_obj.V

############### TRY RECO ###############
########################################


params = Dict{Symbol, Any}()
#params[:reco] = "direct"
params[:reconSize] = acq.encodingSize

params[:reco] = "multiCoilMultiEchoSubspace"
params[:regularization] = "L2"
params[:λ] = 1.e-3
params[:iterations] = 1
params[:solver] = "cgnr"
params[:senseMaps] = sens
params[:basis] = Complex.(V[:,1:5])

im_reco_sub = reconstruction(acq,params)

heatmap(abs.(im_reco_sub.data[:,:,16,5,1,1]),colormap = :greys)


## Apply subspace

I2 = applySubspace(im_reco_sub.data,Complex.(V[:,1:5]))
heatmap(abs.(I2[:,:,16,15,1,1]),colormap = :greys)

im2 = applySubspace(im_reco_sub, Complex.(V[:,1:5]))

heatmap(abs.(im2[:,:,16,15,1,1]),colormap = :greys)

[JakobAsslaender]

Hi,
I already implemented this for non-Cartesian trajectories in MRFingerprintingRecon.jl (not the best name, I know). @andrewwmao is currently working on a Cartesian implementation.

[JeffFessler]

I already implemented this for non-Cartesian trajectories in MRFingerprintingRecon.jl

Would it help discoverability to house it under MRIReco?
Or maybe MRIReco should at least have a "related packages" link somewhere that points to this?

[JakobAsslaender]

Probably :). I guess the best way might be to change ownership to MagneticResonanceImaging? I think I was just hesitant to do so because the package is still WIP. In terms of MRIReco.jl itself, I have no objections against wrappers, but I would be hesitant to move the into MRIReco as I personally prefer a more low-level interface.

[JeffFessler]

oops, i meant MagneticResonanceImaging not MRIReco for exactly the reason you mentioned. just like MRIFieldmaps is under MagneticResonanceImaging and has only low-level interfaces. up to you of course!
everything is a WIP so don't let that deter you :)

[codecov-commenter]

Codecov Report

Patch coverage: 92.30% and project coverage change: +18.59 🎉

Comparison is base (0f33494) 65.56% compared to head (9380e1c) 84.16%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files

@@             Coverage Diff             @@
##           master     #148       +/-   ##
===========================================
+ Coverage   65.56%   84.16%   +18.59%     
===========================================
  Files           9       10        +1     
  Lines         273      322       +49     
===========================================
+ Hits          179      271       +92     
+ Misses         94       51       -43     

Impacted Files	Coverage Δ
src/Reconstruction/IterativeReconstruction.jl	87.41% <88.88%> (+40.13%)	⬆️
src/Reconstruction/Reconstruction.jl	67.85% <100.00%> (+14.01%)	⬆️
src/Tools/Subspace.jl	100.00% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[JakobAsslaender]

Ha, fair point! OK, I will move it.

[aTrotier]

I have some kind of implementation + test ready.
I have to try it on real accelerated k-space to see the results, I guess with fista the results won't be as good as with BART (issue #114 ).

@JakobAsslaender I suppose your implementation is faster than mine. Do you think we should keep the high level implementation ?

[aTrotier]

shape mismatch with wavelet regularization :)

[aTrotier]

I have some weird instability for the reconstruction with multiCoilMultiEcho acquisition.

I tried :

• reconstruction with only the low level functions -> same instability
• ComplexF64 -> same instability
• Real and simulated data -> same instability
• reco with "multiCoil" or "multiEcho" -> stable...

any idea ?

Somebody else can also try to reproduce the error with that branch (if needed I can share the dataset). @alexjaffray maybe you can also try with one of your MESE datasets :/

High and Low level reco :

using CairoMakie
using ImageUtils: shepp_logan
using LinearAlgebra
using MRIReco, MRISimulation, MRICoilSensitivities, MRISampling


T = ComplexF32
N=128
TEnum = collect(7:7:30*7)
x = T.(shepp_logan(N))
rmap = 5.0*abs.(x)

coilsens = T.(birdcageSensitivity(N, 4, 1.5))

TEnum = collect(2.e-2:2.e-2:50.e-2)
# simulation
params = Dict{Symbol, Any}()
params[:simulation] = "fast"
params[:trajName] = "Cartesian"
params[:numProfiles] = floor(Int64, N)
params[:numSamplingPerProfile] = N
params[:r2map] = rmap
params[:T_echo] = TEnum
params[:seqName] = "ME"
params[:refocusingAngles] = Float64.(repeat([pi],length(TEnum)))
params[:senseMaps] = coilsens

acq1 = simulation( real(x), params)
sens = espirit(acq1,eigThresh_2=0.0)

#= REAL DATA

b=BrukerFile("/Users/aurelien/Downloads/MESE_QUEU/")
raw = RawAcquisitionData(b)
acq=AcquisitionData(raw)

# extract 1 slice and take only 10 first echoes

acq1 = deepcopy(acq)
acq1.kdata = acq.kdata[1:30,5:5,:]
acq1.subsampleIndices = acq.subsampleIndices[1:30]
acq1.traj=acq.traj[1:30]

sens = espirit(acq1,eigThresh_2=0.0)
=#

begin
params = Dict{Symbol, Any}()
params[:reconSize] = acq1.encodingSize
params[:reco] = "multiCoilMultiEcho"
#params[:sparseTrafo] = "nothing" #"nothing" #sparse trafo
params[:regularization] = "L2"
params[:λ] = 1.e-3
params[:iterations] = 1
params[:solver] = "cgnr"
params[:senseMaps] = sens

x_approx= reconstruction(acq1,params);


p1 = (64,64)
f=Figure()
ax=Axis(f[1,1])
heatmap!(ax,abs.(x_approx[:,:,1,10,1,1]))
scatter!(ax,p1,color=:red)
ax=Axis(f[1,2])
lines!(ax,abs.(x_approx[p1...,1,:,1,1]))
f
end



## try again with low level
begin
numContr = MRIReco.numContrasts(acq1)
numChan = MRIReco.numChannels(acq1)
shape_ = acq1.encodingSize

#tr = [trajectory(acq_copy,i) for i=1:numContr]

tr = acq1.traj
idx = acq1.subsampleIndices
ftOp = [fourierEncodingOp(acq1.encodingSize, tr[i], "fast",subsampleIdx=idx[i]) for i=1:numContr]

S = SensitivityOp(reshape(sens,:,numChan),length(ftOp))
E = DiagOp(repeat(ftOp, inner=numChan)...) ∘ S


kdata = multiCoilMultiEchoData(acq1, 1);


im_res = E' * vec(kdata);
im_res = abs.(reshape(im_res,shape_...,:));

p1 = (64,64)
f=Figure()
ax=Axis(f[1,1])
heatmap!(ax,im_res[:,:,10])
scatter!(ax,p1,color=:red)
ax=Axis(f[1,2])
lines!(ax,im_res[p1...,:])
#ylims!(ax, 0,60)
f
end

weirdly for real dataset if I run all this code the instability is gone. But if I run multiple times only from the lines im_res = E' * vec(kdata); I can again observe the instability (for simulated data, the instability are present in both case

[aTrotier]

@JakobAsslaender I am not able to get a working "multiCoil" reconstruction with the diagonal operator build like that :

E = DiagOp(repeat(ftOp, inner=numChan)...) ∘ S # do not work

but it works (and is stable) with this line (which correspond to the high level implementation of "multiCoil"-

E = DiagOp(ftOp[1], numChan) ∘ S

It is an issue with "repeat" it works if I replace it by

  ops2 = [copy(ftOp[1]) for n=1:numChan]
  E= DiagOp(ops2...)

I have to check how to change that in the high level interface

[alexjaffray]

Thanks @aTrotier for pushing this forward. As far as the question about repeat() vs. copy(), my intuition is that operators which maintain any computed or mutating internal variables must be copied.

I will pull this and try it out with some of my data

[alexjaffray]

@aTrotier PR #141 was merged, which had a competing change to reconstruction_multiCoilmultiEcho (the addition of the L_inv argument, which already existed for multiCoil recon). Latest commit resolves the conflict, if it causes any problems downstream feel free to roll it back

[aTrotier]

Thank @alexjaffray, I don't see any reason to broke my branch :)
I'll try to finish the literate example in the doc today and then change the PR to review

[aTrotier]

Seems good on my side.

We will have to benchmark the speed / memory / precision against BART

[aTrotier]

Can I merge It ?