diff --git a/Project.toml b/Project.toml index 04e409fb..c7706c75 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MRIReco" uuid = "bdf86e05-2d2b-5731-a332-f3fe1f9e047f" authors = ["Tobias Knopp "] -version = "0.4" +version = "0.4.1" [deps] Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -36,7 +36,7 @@ Graphics = "0.4, 1.0" HDF5 = "0.15" ImageUtils = "0.2.5" LightXML = "0.8.1, 0.9" -LinearOperators = "1.0" +LinearOperators = "2.0" LowRankApprox = "0.4.2,0.5" NFFT = "0.7" NIfTI = "0.4, 0.5" @@ -44,7 +44,7 @@ Polyester = "0.6" ProgressMeter = "1.2" Reexport = "0.2, 1" RegularizedLeastSquares = "0.7.1, 0.8" -SparsityOperators = "0.3" +SparsityOperators = "0.4" StatsBase = "0.32, 0.33" ThreadPools = "1.1.3, 2.0" Unitful = "0.17, 1.2" diff --git a/src/MRIReco.jl b/src/MRIReco.jl index 139047a2..6de42a67 100644 --- a/src/MRIReco.jl +++ b/src/MRIReco.jl @@ -7,6 +7,7 @@ using ProgressMeter using NFFT using Distributions using LinearOperators +using SparsityOperators @reexport using RegularizedLeastSquares using StatsBase using LinearAlgebra diff --git a/src/Operators/Composition.jl b/src/Operators/Composition.jl index 8d16a204..f25b3647 100644 --- a/src/Operators/Composition.jl +++ b/src/Operators/Composition.jl @@ -11,12 +11,17 @@ mutable struct CompositeOp{T} <: AbstractLinearOperator{T} ncol :: Int symmetric :: Bool hermitian :: Bool - prod :: Function - tprod :: Function - ctprod :: Function + prod! :: Function + tprod! :: Function + ctprod! :: Function nprod :: Int ntprod :: Int nctprod :: Int + args5 :: Bool + use_prod5! :: Bool + allocated5 :: Bool + Mv5 :: Vector{T} + Mtu5 :: Vector{T} isWeighting :: Bool A B @@ -32,22 +37,24 @@ function CompositeOp(A,B;isWeighting=false) ncol=B.ncol S = eltype(A) - function produ(x::Vector{T}) where T<:Union{Real,Complex} - return A*(B*x) + function produ!(res, x::Vector{T}) where T<:Union{Real,Complex} + res .= A*(B*x) end - function tprodu(y::Vector{T}) where T<:Union{Real,Complex} - return transposed(B)*(transposed(A)*y) + function tprodu!(res, y::Vector{T}) where T<:Union{Real,Complex} + res .= transposed(B)*(transposed(A)*y) end - function ctprodu(y::Vector{T}) where T<:Union{Real,Complex} - return adjoint(B)*(adjoint(A)*y) + function ctprodu!(res, y::Vector{T}) where T<:Union{Real,Complex} + res .= adjoint(B)*(adjoint(A)*y) end Op = CompositeOp{S}( nrow, ncol, false, false, - produ, - tprodu, - ctprodu, 0, 0, 0, isWeighting, A, B ) + produ!, + tprodu!, + ctprodu!, + 0, 0, 0, false, false, false, S[], S[], + isWeighting, A, B ) return Op end diff --git a/src/Operators/ExplicitOp.jl b/src/Operators/ExplicitOp.jl index 8c5e9ac9..768ad0bf 100644 --- a/src/Operators/ExplicitOp.jl +++ b/src/Operators/ExplicitOp.jl @@ -5,12 +5,17 @@ mutable struct ExplicitOp{T,F1,F2} <: AbstractLinearOperator{T} ncol :: Int symmetric :: Bool hermitian :: Bool - prod :: Function - tprod :: F1 - ctprod :: F2 + prod! :: Function + tprod! :: F1 + ctprod! :: F2 nprod :: Int ntprod :: Int nctprod :: Int + args5 :: Bool + use_prod5! :: Bool + allocated5 :: Bool + Mv5 :: Vector{T} + Mtu5 :: Vector{T} end """ @@ -42,9 +47,10 @@ function ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory, correctionmap::Array end return ExplicitOp{ComplexF64,Nothing,Function}(nrow, ncol, false, false - , x->produ(x, shape, nodes, times, echoOffset, correctionmap) + , (res,x)->(res .= produ(x, shape, nodes, times, echoOffset, correctionmap)) , nothing - , y->ctprodu(y, shape, nodes, times, echoOffset, correctionmap),0,0,0) + , (res,y)->(res .= ctprodu(y, shape, nodes, times, echoOffset, correctionmap)) + , 0,0,0, false, false, false, ComplexF64[], ComplexF64[]) end function produ(x::Vector{T}, shape::NTuple{2,Int64}, @@ -148,5 +154,5 @@ end function adjoint(op::ExplicitOp{T}) where T return LinearOperator{T}(op.ncol, op.nrow, op.symmetric, op.hermitian, - op.ctprod, nothing, op.prod) + op.ctprod!, nothing, op.prod!) end diff --git a/src/Operators/FieldmapNFFTOp.jl b/src/Operators/FieldmapNFFTOp.jl index 2524e21b..f9a9236a 100644 --- a/src/Operators/FieldmapNFFTOp.jl +++ b/src/Operators/FieldmapNFFTOp.jl @@ -15,12 +15,17 @@ mutable struct FieldmapNFFTOp{T,F1,F2,D} <:AbstractLinearOperator{T} ncol :: Int symmetric :: Bool hermitian :: Bool - prod :: Function - tprod :: F1 - ctprod :: F2 + prod! :: Function + tprod! :: F1 + ctprod! :: F2 nprod :: Int ntprod :: Int nctprod :: Int + args5 :: Bool + use_prod5! :: Bool + allocated5 :: Bool + Mv5 :: Vector{T} + Mtu5 :: Vector{T} plans idx::Vector{Vector{Int64}} circTraj::Bool @@ -89,17 +94,16 @@ function FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory, circTraj = isCircular(tr) - mul(x::Vector{T}) where T<:ComplexF64 = - produ(x,nrow,ncol,shape,plans,idx,cparam,circTraj,d) - ctmul(y::Vector{T}) where T<:ComplexF64 = - ctprodu(y,shape,plans,idx,cparam,circTraj,d) - inverse(y::Vector{T}) where T<:ComplexF64 = - inv(y,shape,plans,idx,cparam,circTraj,p,y,d) + mul!(res, x::Vector{T}) where T<:ComplexF64 = + (res .= produ(x,nrow,ncol,shape,plans,idx,cparam,circTraj,d)) + ctmul!(res, y::Vector{T}) where T<:ComplexF64 = + (res .= ctprodu(y,shape,plans,idx,cparam,circTraj,d)) return FieldmapNFFTOp{ComplexF64,Nothing,Function,D}(nrow, ncol, false, false - , mul + , mul! , nothing - , ctmul, 0, 0, 0, plans, idx, circTraj, shape, cparam) + , ctmul!, 0, 0, 0, false, false, false, ComplexF64[], ComplexF64[] + , plans, idx, circTraj, shape, cparam) end function Base.copy(S::FieldmapNFFTOp) @@ -111,18 +115,17 @@ function Base.copy(S::FieldmapNFFTOp) cparam = deepcopy(S.cparam) - mul(x::Vector{T}) where T<:ComplexF64 = - produ(x,S.nrow,S.ncol,S.shape,plans,idx,cparam,S.circTraj,d) - ctmul(y::Vector{T}) where T<:ComplexF64 = - ctprodu(y,S.shape,plans,idx,cparam,S.circTraj,d) - inverse(y::Vector{T}) where T<:ComplexF64 = - inv(y,S.shape,plans,idx,cparam,S.circTraj,p,y,d) + mul!(res, x::Vector{T}) where T<:ComplexF64 = + (res .= produ(x,S.nrow,S.ncol,S.shape,plans,idx,cparam,S.circTraj,d)) + ctmul!(res, y::Vector{T}) where T<:ComplexF64 = + (res .= ctprodu(y,S.shape,plans,idx,cparam,S.circTraj,d)) D = length(S.shape) return FieldmapNFFTOp{ComplexF64,Nothing,Function,D}(S.nrow, S.ncol, false, false - , mul + , mul! , nothing - , ctmul, 0, 0, 0, plans, idx, S.circTraj, S.shape, cparam) + , ctmul!, 0, 0, 0, false, false, false, ComplexF64[], ComplexF64[] + , plans, idx, S.circTraj, S.shape, cparam) end # function produ{T<:ComplexF64}(x::Vector{T}, numOfNodes::Int, numOfPixel::Int, shape::Tuple, plan::Vector{NFFTPlan{Float64,2,1}}, cparam::InhomogeneityData, density::Vector{Float64}, symmetrize::Bool) diff --git a/src/Operators/MapSliceOp.jl b/src/Operators/MapSliceOp.jl index 49fc927d..8dcfcfba 100644 --- a/src/Operators/MapSliceOp.jl +++ b/src/Operators/MapSliceOp.jl @@ -21,7 +21,7 @@ generates an operator that applies `trafo` to each slice of dimension `dim` of a """ function MapSliceOp(trafo, dim::Int64, size1::Tuple, size2::Tuple; T=ComplexF64) return LinearOperator(prod(size2), prod(size1), false, false - , x->mapSliceForeward(trafo, x, size1, dim) + , (res,x) -> ( res .= mapSliceForeward(trafo, x, size1, dim) ) , nothing - , y->mapSliceBackward(trafo, y, size2, dim) ) + , (res,y) -> ( res .= mapSliceBackward(trafo, y, size2, dim) ) ) end diff --git a/src/Operators/NFFTOp.jl b/src/Operators/NFFTOp.jl index a4286bba..b1e759d2 100644 --- a/src/Operators/NFFTOp.jl +++ b/src/Operators/NFFTOp.jl @@ -6,12 +6,17 @@ mutable struct NFFTOp{T} <: AbstractLinearOperator{T} ncol :: Int symmetric :: Bool hermitian :: Bool - prod :: Function - tprod :: Nothing - ctprod :: Function + prod! :: Function + tprod! :: Nothing + ctprod! :: Function nprod :: Int ntprod :: Int nctprod :: Int + args5 :: Bool + use_prod5! :: Bool + allocated5 :: Bool + Mv5 :: Vector{T} + Mtu5 :: Vector{T} plan toeplitz :: Bool end @@ -34,9 +39,11 @@ function NFFTOp(shape::Tuple, tr; nodes=nothing, toeplitz=false, plan = plan_nfft(nodes, shape, m=kernelSize, σ=oversamplingFactor, precompute=NFFT.FULL) return NFFTOp{ComplexF64}(size(nodes,2), prod(shape), false, false - , x->produ(plan,x) + , (res,x) -> (res .= produ(plan,x)) , nothing - , y->ctprodu(plan,y), 0, 0, 0, plan, toeplitz) + , (res,y) -> (res .= ctprodu(plan,y)) + , 0, 0, 0, false, false, false, ComplexF64[], ComplexF64[] + , plan, toeplitz) end function produ(plan::NFFT.NFFTPlan, x::Vector{T}) where T<:Union{Real,Complex} @@ -53,9 +60,11 @@ end function Base.copy(S::NFFTOp) plan = copy(S.plan) return NFFTOp{ComplexF64}(size(plan.x,2), prod(plan.N), false, false - , x->produ(plan,x) - , nothing - , y->ctprodu(plan,y), 0, 0, 0, plan, S.toeplitz) + , (res,x) -> (res .= produ(plan,x)) + , nothing + , (res,y) -> (res .= ctprodu(plan,y)) + , 0, 0, 0, false, false, false, ComplexF64[], ComplexF64[] + , plan, S.toeplitz) end ### Normal Matrix Code ### @@ -146,7 +155,7 @@ function diagonalizeOp(p::NFFT.NFFTPlan, weights=nothing) firstRow = nfft_adjoint(p2, nfft(p2,e1)[:].*weights) firstRow = reshape(firstRow, shape) - # construct FT matrix of the eigentvalues + # construct FT matrix of the eigenvalues # here a row and column in the center is filled with zeros. # this is done to ensure that the FFT can operate # on arrays with even size in each dimension diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index 1d5ceae8..7c2544c3 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -76,12 +76,17 @@ mutable struct DiagOp{T} <: AbstractLinearOperator{T} ncol :: Int symmetric :: Bool hermitian :: Bool - prod :: Function - tprod :: Function - ctprod :: Function + prod! :: Function + tprod! :: Function + ctprod! :: Function nprod :: Int ntprod :: Int nctprod :: Int + args5 :: Bool + use_prod5! :: Bool + allocated5 :: Bool + Mv5 :: Vector{T} + Mtu5 :: Vector{T} ops equalOps :: Bool xIdx :: Vector{Int} @@ -108,9 +113,10 @@ function diagOp(ops :: AbstractLinearOperator...) yIdx = cumsum(vcat(1,[ops[i].nrow for i=1:length(ops)])) Op = DiagOp{S}( nrow, ncol, false, false, - x->diagOpProd(x,nrow,xIdx,yIdx,ops...), - y->diagOpTProd(y,ncol,yIdx,xIdx,ops...), - y->diagOpCTProd(y,ncol,yIdx,xIdx,ops...), 0, 0, 0, + (res,x) -> (res .= diagOpProd(x,nrow,xIdx,yIdx,ops...)), + (res,y) -> (res .= diagOpTProd(y,ncol,yIdx,xIdx,ops...)), + (res,y) -> (res .= diagOpCTProd(y,ncol,yIdx,xIdx,ops...)), + 0, 0, 0, false, false, false, S[], S[], [ops...], false, xIdx, yIdx) return Op @@ -126,9 +132,10 @@ function diagOp(op::AbstractLinearOperator, N=1) yIdx = cumsum(vcat(1,[ops[i].nrow for i=1:length(ops)])) Op = DiagOp{S}( nrow, ncol, false, false, - x->diagOpProd(x,nrow,xIdx,yIdx,ops...), - y->diagOpTProd(y,ncol,yIdx,xIdx,ops...), - y->diagOpCTProd(y,ncol,yIdx,xIdx,ops...), 0, 0, 0, + (res,x) -> (res .= diagOpProd(x,nrow,xIdx,yIdx,ops...)), + (res,y) -> (res .= diagOpTProd(y,ncol,yIdx,xIdx,ops...)), + (res,y) -> (res .= diagOpCTProd(y,ncol,yIdx,xIdx,ops...)), + 0, 0, 0, false, false, false, S[], S[], ops, true, xIdx, yIdx ) return Op diff --git a/src/Operators/SamplingOp.jl b/src/Operators/SamplingOp.jl index ef437166..370fdc15 100644 --- a/src/Operators/SamplingOp.jl +++ b/src/Operators/SamplingOp.jl @@ -26,9 +26,9 @@ end function SamplingOp(pattern::Array{Bool}) return LinearOperator(length(pattern), length(pattern), false, false - , x->vec(pattern).*x + , (res,x) -> (res .= vec(pattern).*x) , nothing - , y->vec(pattern).*y ) + , (res,y) -> (res .= vec(pattern).*y) ) end function zeropad(x::Array{T}, pattern::Matrix{Bool}) where T diff --git a/src/Operators/SensitivityOp.jl b/src/Operators/SensitivityOp.jl index d9616e2e..18b859fb 100644 --- a/src/Operators/SensitivityOp.jl +++ b/src/Operators/SensitivityOp.jl @@ -46,10 +46,10 @@ the coil sensitivities specified in `sensMaps` function SensitivityOp(sensMaps::Matrix{ComplexF64}, numContr::Int=1) numVox, numChan = size(sensMaps) sensMapsC = conj.(sensMaps) - return LinearOperator(numVox*numContr*numChan,numVox*numContr,false,false, - x->prod_smap(sensMaps,x,numVox,numChan,numContr), + return LinearOperator{ComplexF64}(numVox*numContr*numChan, numVox*numContr, false, false, + (res,x) -> (res .= prod_smap(sensMaps,x,numVox,numChan,numContr)), nothing, - x->ctprod_smap(sensMapsC,x,numVox,numChan,numContr)) + (res,x) -> (res .= ctprod_smap(sensMapsC,x,numVox,numChan,numContr))) end """