Draft for discussion : Subspace implementation

aTrotier commented 1 year ago

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 #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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

Codecov Report

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

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

:exclamation: 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

```diff @@ 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](https://app.codecov.io/gh/MagneticResonanceImaging/MRIReco.jl/pull/148?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MagneticResonanceImaging) | Coverage Δ | | |---|---|---| | [src/Reconstruction/IterativeReconstruction.jl](https://app.codecov.io/gh/MagneticResonanceImaging/MRIReco.jl/pull/148?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MagneticResonanceImaging#diff-c3JjL1JlY29uc3RydWN0aW9uL0l0ZXJhdGl2ZVJlY29uc3RydWN0aW9uLmps) | `87.41% <88.88%> (+40.13%)` | :arrow_up: | | [src/Reconstruction/Reconstruction.jl](https://app.codecov.io/gh/MagneticResonanceImaging/MRIReco.jl/pull/148?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MagneticResonanceImaging#diff-c3JjL1JlY29uc3RydWN0aW9uL1JlY29uc3RydWN0aW9uLmps) | `67.85% <100.00%> (+14.01%)` | :arrow_up: | | [src/Tools/Subspace.jl](https://app.codecov.io/gh/MagneticResonanceImaging/MRIReco.jl/pull/148?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MagneticResonanceImaging#diff-c3JjL1Rvb2xzL1N1YnNwYWNlLmps) | `100.00% <100.00%> (ø)` | |

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

JakobAsslaender commented 1 year ago

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

aTrotier commented 1 year ago

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 commented 1 year ago

shape mismatch with wavelet regularization :)

aTrotier commented 1 year ago

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 commented 1 year ago

@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 commented 1 year ago

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 commented 1 year ago

@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 commented 1 year ago

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 commented 1 year ago

Seems good on my side.

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

aTrotier commented 10 months ago

Can I merge It ?

MagneticResonanceImaging / MRIReco.jl

Draft for discussion : Subspace implementation #148

Codecov Report