Skip to content

Commit

Permalink
Merge pull request #209 from MagneticResonanceImaging/floop_along_con…
Browse files Browse the repository at this point in the history
…trast

* floop along contrast in MultiCoil reco

* fix remaining end
  • Loading branch information
aTrotier authored Dec 11, 2024
2 parents 4548fd8 + 3d3d735 commit 29af5ed
Showing 1 changed file with 13 additions and 15 deletions.
28 changes: 13 additions & 15 deletions src/Reconstruction/IterativeReconstruction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -208,30 +208,28 @@ function reconstruction_multiCoil(acqData::AcquisitionData{T}
# solve optimization problem
Ireco = zeros(Complex{T}, prod(reconSize), numSl, numContr, numRep)
let reg = reg # Fix @floop warning due to conditional/multiple assignment to reg
@floop executor(arrayType) for l = 1:numRep, k = 1:numSl
@floop executor(arrayType) for l = 1:numRep, k = 1:numSl, j = 1:numContr
if encodingOps != nothing
E = encodingOps[:,k]
else
E = encodingOps_parallel(acqData, reconSize, senseMapsUnCorr; slice=k, encParams...)
end

for j = 1:numContr
W = WeightingOp(Complex{T}; weights=weights[j], rep=numChan)
kdata = arrayType((multiCoilData(acqData, j, k, rep=l))) .* repeat(weights[j], numChan)
if !isnothing(L_inv)
kdata = vec(reshape(kdata, :, numChan) * L_inv')
end
W = WeightingOp(Complex{T}; weights=weights[j], rep=numChan)
kdata = arrayType((multiCoilData(acqData, j, k, rep=l))) .* repeat(weights[j], numChan)
if !isnothing(L_inv)
kdata = vec(reshape(kdata, :, numChan) * L_inv')
end

EFull = (W, E[j])
EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...)
solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...)
I = solve!(solv, kdata)
EFull = (W, E[j])
EFullᴴEFull = normalOperator(EFull; normalOpParams(arrayType)...)
solv = createLinearSolver(solver, EFull; AHA=EFullᴴEFull, reg=reg, params...)
I = solve!(solv, kdata)

if isCircular( trajectory(acqData, j) )
circularShutter!(reshape(I, reconSize), 1.0)
end
Ireco[:,k,j,l] = Array(I)
if isCircular( trajectory(acqData, j) )
circularShutter!(reshape(I, reconSize), 1.0)
end
Ireco[:,k,j,l] = Array(I)
end
end
Ireco_ = reshape(Ireco, volumeSize(reconSize, numSl)..., numContr, 1,numRep)
Expand Down

0 comments on commit 29af5ed

Please sign in to comment.