Merge remote-tracking branch 'origin/subspace' into subspace

MRIOperators/src/FieldmapNFFTOp.jl

@@ -33,6 +33,7 @@ mutable struct FieldmapNFFTOp{T,F1,F2,D} <:AbstractLinearOperator{Complex{T}}
  circTraj::Bool
  shape::NTuple{D,Int64}
  cparam::InhomogeneityData{T}
+
 end
 
 LinearOperators.storage_type(op::FieldmapNFFTOp) = typeof(op.Mv5)
@@ -81,6 +82,9 @@ function FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory,
  nrow = size(nodes,2)
  ncol = prod(shape)
 
+ x_tmp = zeros(Complex{T}, ncol)
+ y_tmp = zeros(Complex{T}, nrow)
+
  # create and truncate low-rank expansion
  cparam = createInhomogeneityData_(vec(times), correctionmap; K=K, alpha=alpha, m=m, method=method, K_tol=K_tol, numSamp=numSamplingPerProfile(tr),step=step)
  K = size(cparam.A_k,2)
@@ -93,45 +97,51 @@ function FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory,
  idx[κ] = findall(x->x!=0.0, cparam.A_k[:,κ])
  plans[κ] = plan_nfft(nodes[:,idx[κ]], shape, m=3, σ=1.25, precompute = NFFT.POLYNOMIAL)
  end
-
- d = [zeros(ComplexF64, length(idx[κ])) for κ=1:K ]
- p = [zeros(Complex{T}, shape) for κ=1:K]
- x_tmp = zeros(Complex{T}, ncol)
- y_tmp = zeros(Complex{T}, nrow)
 
- circTraj = isCircular(tr)
+ d = [zeros(Complex{T}, length(idx[κ])) for κ=1:K ]
+ p = [zeros(Complex{T}, shape) for κ=1:K]
 
- mul! = (res,x) -> produ!(res,x,x_tmp,shape,plans,idx,cparam,circTraj,d,p)
- ctmul! = (res,y) -> ctprodu!(res,y,y_tmp,shape,plans,idx,cparam,circTraj,d,p)
+ circTraj = isCircular(tr)
 
  return FieldmapNFFTOp{T,Nothing,Function,D}(nrow, ncol, false, false
- , mul!
+ , (res,x) -> produ!(res,x,x_tmp,shape,plans,idx,cparam,circTraj,d,p)
  , nothing
- , ctmul!, 0, 0, 0, false, false, false, ComplexF64[], ComplexF64[]
+ , (res,y) -> ctprodu!(res,y,y_tmp,shape,plans,idx,cparam,circTraj,d,p), 0, 0, 0, false, false, false, Complex{T}[], Complex{T}[]
  , plans, idx, circTraj, shape, cparam)
 end
 
+function Base.copy(cparam::InhomogeneityData{T}) where T
+
+ return cparam
+
+end
+
 function Base.copy(S::FieldmapNFFTOp{T,Nothing,Function,D}) where {T,D}
+
+ shape = S.shape
+
  K=length(S.plans)
  plans = [copy(S.plans[i]) for i=1:K]
- idx = deepcopy(S.idx)
+ idx = copy(S.idx)
 
- d = [zeros(Complex{T}, length(idx[κ])) for κ=1:K ]
- p = [zeros(Complex{T}, S.shape) for κ=1:K]
  x_tmp = zeros(Complex{T}, S.ncol)
  y_tmp = zeros(Complex{T}, S.nrow)
 
- cparam = deepcopy(S.cparam)
+ cparam = copy(S.cparam)
+ d = [zeros(Complex{T}, length(idx[κ])) for κ=1:K]
+ p = [zeros(Complex{T}, shape) for κ=1:K]
+
+ D_ = length(shape)
+ circTraj = S.circTraj
 
- mul! = (res,x) -> produ!(res,x,x_tmp,S.shape,plans,idx,cparam,S.circTraj,d,p)
- ctmul! = (res,y) -> ctprodu!(res,y,y_tmp,S.shape,plans,idx,cparam,S.circTraj,d,p)
+ mul! = (res,x) -> produ!(res,x,x_tmp,shape,plans,idx,cparam,circTraj,d,p)
+ ctmul! = (res,y) -> ctprodu!(res,y,y_tmp,shape,plans,idx,cparam,circTraj,d,p)
 
- D_ = length(S.shape)
  return FieldmapNFFTOp{T,Nothing,Function,D_}(S.nrow, S.ncol, false, false
  , mul!
  , nothing
  , ctmul!, 0, 0, 0, false, false, false, Complex{T}[], Complex{T}[]
- , plans, idx, S.circTraj, S.shape, cparam)
+ , plans, idx, circTraj, shape, cparam)
 end
 
 function produ!(s::AbstractVector{T}, x::AbstractVector{T}, x_tmp::Vector{T},shape::Tuple, plan,

MRIOperators/test/testOperators.jl

@@ -169,6 +169,34 @@ function testSparseOp(T::Type,shape)
  @test (norm(xapprox-x)/norm(x)) < 1e-3
 end
 
+## test FieldmapNFFTOp
+function testCopySizes(N=16)
+ # random image
+ x = zeros(ComplexF64,N,N)
+ for i=1:N,j=1:N
+ x[i,j] = rand()
+ end
+
+ tr = CartesianTrajectory(Float64,N,N;TE=0.0,AQ=0.01)
+ times = readoutTimes(tr)
+ nodes = kspaceNodes(tr)
+ cmap = im*quadraticFieldmap(N,N)[:,:,1]
+
+ # FourierMatrix
+ idx = CartesianIndices((N,N))[collect(1:N^2)]
+
+ # Operators
+
+ F_nfft = NFFTOp((N,N),tr,symmetrize=false)
+ F_fmap_nfft = FieldmapNFFTOp((N,N),tr,cmap,symmetrize=false)
+
+ # Copy the FieldmapNFFTOp operator and change the plans field of the new operator to empty 
+ F_fmap_nfft_copy = copy(F_fmap_nfft)
+ F_nfft_copy = copy(F_nfft)
+
+ @test sizeof(F_fmap_nfft_copy) == sizeof(F_fmap_nfft)
+
+end
 
 function testOperators()
  @testset "Linear Operator" begin

Project.toml

@@ -9,11 +9,11 @@ FLoops = "cc61a311-1640-44b5-9fba-1b764f453329"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MRIBase = "f7771a9a-6e57-4e71-863b-6e4b6a2f17df"
 MRIOperators = "fb1137e3-90a6-46ce-a672-6e1e53d120f2"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RegularizedLeastSquares = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23"
-PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
@@ -25,10 +25,10 @@ MRIBase = "0.3.0"
 MRIOperators = "0.1"
 MRISampling = "0.1"
 MRISimulation = "0.1"
+PrecompileTools = "1"
 ProgressMeter = "1.2"
 Reexport = "0.2, 1"
-RegularizedLeastSquares = "0.9.0"
-PrecompileTools = "1"
+RegularizedLeastSquares = "0.9.0, 0.10"
 Unitful = "1.2"
 julia = "1.6"
 

src/Reconstruction/IterativeReconstruction.jl

@@ -235,6 +235,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
@@ -249,6 +250,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
  if isnothing(sparseTrafo)
  sparseTrafo = SparseOp(Complex{T},"nothing", reconSize; params...)
@@ -262,10 +266,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)

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
 

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)
  elseif recoParams[:reco] == "multiCoilMultiEchoSubspace"
  return reconstruction_multiCoilMultiEcho_subspace(acqData, reconSize[1:encodingDims], reg, sparseTrafo, weights, solvername, senseMaps, normalize, encOps, recoParams)
  else