MagneticResonanceImaging / MRIReco.jl

Julia Package for MRI Reconstruction
https://magneticresonanceimaging.github.io/MRIReco.jl/latest/
Other
85 stars 22 forks source link

export function and correct SparseOp for echoes #146

Closed aTrotier closed 1 year ago

aTrotier commented 1 year ago

Correct some errors for multiCoilMultiEcho reconstruction.

Instability in the reconstruction

Now the reconstruction returns images but I have a stability / reproductibility issue when the number of coils > 1

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

N = 64
nCoil = 4
T = Complex{Float32}

x = T.(shepp_logan(N))

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

# simulation
params = Dict{Symbol, Any}()
params[:simulation] = "fast"
params[:trajName] = "Cartesian"
params[:numProfiles] = floor(Int64, N)
params[:numSamplingPerProfile] = N
params[:senseMaps] = coilsens

acqData = simulation( real(x), params )

params = Dict{Symbol, Any}()
params[:reconSize] = (N,N)
params[:reco] = "multiCoil"
params[:regularization] = "L2"
params[:λ] = 1.e-3
params[:iterations] = 1
params[:solver] = "cgnr"
params[:senseMaps] = reshape(coilsens,N,N,1,nCoil)

x_approx = reshape(reconstruction(acqData,params),N,N,:)

# Reco multiCoilMultiEcho
params2 = Dict{Symbol, Any}()
params2[:reconSize] = (N,N)
params2[:reco] = "multiCoilMultiEcho"
#params[:sparseTrafo] = "nothing" #"nothing" #sparse trafo
params2[:regularization] = "L2"
params2[:λ] = 1.e-3
params2[:iterations] = 1
params2[:solver] = "cgnr"
params2[:senseMaps] = reshape(coilsens,N,N,1,nCoil)

x_approx2= reshape(reconstruction(acqData,params2),N,N,:)

f=Figure()
for i = 1:8
    x_approx = reshape(reconstruction(acqData,params),N,N,:)
    heatmap!(Axis(f[i,1],aspect = 1),abs.(x_approx[:,:,1]))
    hidedecorations!(current_axis())

    x_approx2= reshape(reconstruction(acqData,params2),N,N,:)
    heatmap!(Axis(f[i,2],aspect=1),abs.(x_approx2[:,:,1]))
    hidedecorations!(current_axis())
end
f

which returns (left : MultiCoil reco / right : multiCoilEchoes):

image

Decomposition of the operator

I also get that error with the code bellow (which only used the operators)

## test operator
kdata = multiCoilMultiEchoData(acqData, 1)
encParams = MRIReco.getEncodingOperatorParams(;params...)
E = encodingOp_multiEcho_parallel(acqData, (N,N), reshape(coilsens,N,N,1,nCoil); slice=1)

f = Figure()
for i = 1:4,j=1:4
    im_reco = E'*kdata
    im_reco = reshape(im_reco,N,N,:)

    heatmap!(Axis(f[i,j],aspect=1),abs.(im_reco[:,:,1]))
    hidedecorations!(current_axis())
end
f
image

I don't think it is an issue with the simulation because I also get the same kind of results on my real data

aTrotier commented 1 year ago

My reconstruction (on real data) works with the line :


  return diagOp(repeat(ft, inner=numChan)...) ∘ S

rather than outer

However I still have an issue with the data from simulation... I don't really know why.

aTrotier commented 1 year ago

I can share a dataset if someone else want to test if the multiCoilMultiEcho reconstruction works on their side.