Merge pull request #141 from MagneticResonanceImaging/SNR_opt_reco

Noise decorrelation implementation for multiCoilMultiEcho reconstruction
MagneticResonanceImaging · Jul 12, 2023 · 0f33494 · 0f33494
2 parents 8a60db9 + 72520af
commit 0f33494
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 3 deletions.
diff --git a/src/Reconstruction/IterativeReconstruction.jl b/src/Reconstruction/IterativeReconstruction.jl
@@ -232,6 +232,7 @@ function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T}
                               , reg::Vector{Regularization}
                               , sparseTrafo
                               , weights::Vector{Vector{Complex{T}}}
+                              , L_inv::Union{LowerTriangular{Complex{T}, Matrix{Complex{T}}}, Nothing}
                               , solvername::String
                               , senseMaps::Array{Complex{T}}
                               , normalize::Bool=false
@@ -246,6 +247,9 @@ function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T}
   numContr, numChan, numSl, numRep = numContrasts(acqData), numChannels(acqData), numSlices(acqData), numRepetitions(acqData)
   encParams = getEncodingOperatorParams(;params...)
 
+  # noise decorrelation
+  senseMapsUnCorr = decorrelateSenseMaps(L_inv, senseMaps, numChan)
+
   # set sparse trafo in reg
   reg[1].params[:sparseTrafo] = DiagOp( repeat([sparseTrafo],numContr)... )
 
@@ -256,10 +260,13 @@ function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData{T}
     if encodingOps != nothing
       E = encodingOps[i]
     else
-      E = encodingOp_multiEcho_parallel(acqData, reconSize, senseMaps; slice=i, encParams...)
+      E = encodingOp_multiEcho_parallel(acqData, reconSize, senseMapsUnCorr; slice=i, encParams...)
     end
 
     kdata = multiCoilMultiEchoData(acqData, i) .* repeat(vcat(weights...), numChan)
+    if !isnothing(L_inv)
+      kdata = vec(reshape(kdata, :, numChan) * L_inv')
+    end
 
     EFull = ∘(W, E)#, isWeighting=true)
     EFullᴴEFull = normalOperator(EFull)

diff --git a/src/Reconstruction/RecoParameters.jl b/src/Reconstruction/RecoParameters.jl
@@ -128,7 +128,8 @@ function setupIterativeReco(acqData::AcquisitionData{T}, recoParams::Dict) where
   if isempty(noiseData)
     L_inv = nothing
   else
-    L = cholesky(covariance(noiseData), check = true)
+    psi = convert(typeof(noiseData), covariance(noiseData))
+    L = cholesky(psi, check = true)
     L_inv = inv(L.L) #noise decorrelation matrix
   end
 

diff --git a/src/Reconstruction/Reconstruction.jl b/src/Reconstruction/Reconstruction.jl
@@ -46,7 +46,7 @@ function reconstruction(acqData::AcquisitionData, recoParams::Dict)
     elseif recoParams[:reco] == "multiCoil"
         return reconstruction_multiCoil(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, L_inv, solvername, senseMaps, normalize, encOps, recoParams)
     elseif recoParams[:reco] == "multiCoilMultiEcho"
-        return reconstruction_multiCoilMultiEcho(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, solvername, senseMaps, normalize, encOps, recoParams)
+        return reconstruction_multiCoilMultiEcho(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, L_inv, solvername, senseMaps, normalize, encOps, recoParams)
     else
         @error "reco modell $(recoParams[:reco]) not found"
     end