Merge pull request #47 from MagneticResonanceImaging/LinearOperators2

Linear operators2

Project.toml

@@ -1,7 +1,7 @@
 name = "MRIReco"
 uuid = "bdf86e05-2d2b-5731-a332-f3fe1f9e047f"
 authors = ["Tobias Knopp <[email protected]>"]
-version = "0.4"
+version = "0.4.1"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -36,15 +36,15 @@ 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"
 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"

src/MRIReco.jl

@@ -7,6 +7,7 @@ using ProgressMeter
 using NFFT
 using Distributions
 using LinearOperators
+using SparsityOperators
 @reexport using RegularizedLeastSquares
 using StatsBase
 using LinearAlgebra

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

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

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)

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

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

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

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

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
 
 """