Merge pull request #148 from MagneticResonanceImaging/subspace

Draft for discussion : Subspace implementation
MagneticResonanceImaging · Aug 29, 2023 · 5a8e4e1 · 5a8e4e1
2 parents 62e2077 + 9380e1c
commit 5a8e4e1
Show file tree

Hide file tree

Showing 12 changed files with 631 additions and 46 deletions.
diff --git a/MRIBase/src/Datatypes/AcqData.jl b/MRIBase/src/Datatypes/AcqData.jl
@@ -1,7 +1,7 @@
 export AcquisitionData, kData, kDataCart, kdataSingleSlice, convertUndersampledData,
        changeEncodingSize2D, convert3dTo2d, samplingDensity,
        numContrasts, numChannels, numSlices, numRepetitions,
-       encodingSize, fieldOfView, multiCoilData
+       encodingSize, fieldOfView, multiCoilData, multiCoilMultiEchoData, multiEchoData
 
 """
 struct describing MRI acquisition data.
@@ -112,8 +112,9 @@ function AcquisitionData(kspace::Array{Complex{T},6}) where T
         tr = MRIBase.CartesianTrajectory3D(T, sy, sx, numSlices=sz, TE=T(0), AQ=T(0))
     end
 
-    kdata = [reshape(kspace,:,nCh) for i=1:nEchos, j=1:1, k=1:nReps]
-    acq = AcquisitionData(tr, kdata)
+    kdata = [reshape(kspace[:,:,:,:,i,k],:,nCh) for i=1:nEchos, j=1:1, k=1:nReps]
+    traj = [tr for i=1:nEchos]
+    acq = AcquisitionData(traj, kdata, encodingSize=ntuple(d->0, ndims(tr)))
 
     acq.encodingSize = ndims(tr) == 3 ? (sx,sy,sz) : (sx,sy)
 
@@ -185,7 +186,7 @@ end
 returns the k-space contained in `acqData` for all echoes and given `coil`, `slice` and `rep`(etition).
 """
 function multiEchoData(acqData::AcquisitionData{T}, coil::Int64, slice::Int64;rep::Int64=1) where T
-  kdata = T[]
+  kdata = Complex{T}[]
   for echo=1:numContrasts(acqData)
     append!(kdata,acqData.kdata[echo,slice,rep][:,coil])
   end
@@ -208,10 +209,8 @@ returns the k-space contained in `acqData` for all coils, echoes and given `slic
 """
 function multiCoilMultiEchoData(acqData::AcquisitionData{T},slice::Int64;rep=1) where T
   kdata = Complex{T}[]
-  for coil=1:numChannels(acqData)
-    for echo=1:numContrasts(acqData)
-      append!(kdata, acqData.kdata[echo,slice,rep][:,coil])
-    end
+  for echo=1:numContrasts(acqData)
+      append!(kdata, vec(acqData.kdata[echo,slice,rep]))
   end
   return kdata
 end

diff --git a/MRIOperators/src/EncodingOp.jl b/MRIOperators/src/EncodingOp.jl
@@ -90,7 +90,8 @@ function encodingOp_multiEcho_parallel(acqData::AcquisitionData{T,D}, shape::NTu
   # fourier operators
   ft = encodingOps_simple(acqData, shape; kargs...)
   S = SensitivityOp(reshape(smaps,:,numChan),numContrasts(acqData))
-  return DiagOp(repeat(ft, outer=numChan)...) ∘ S
+  ops2 = [copy(ft[n]) for j=1:numChan,n=eachindex(ft)]
+  return DiagOp(ops2...) ∘ S
 end
 
 ###################################

diff --git a/MRIOperators/src/MRIOperators.jl b/MRIOperators/src/MRIOperators.jl
@@ -23,6 +23,7 @@ include("MapSliceOp.jl")
 include("FieldmapNFFTOp.jl")
 include("EncodingOp.jl")
 include("SparseOp.jl")
+include("SubspaceOp.jl")
 
 """
     hcat(A::AbstractLinearOperator, n::Int)

diff --git a/MRIOperators/src/SensitivityOp.jl b/MRIOperators/src/SensitivityOp.jl
@@ -2,25 +2,25 @@ export SensitivityOp
 
 function prod_smap!(y::AbstractVector{T}, smaps::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numChan, numContr=1) where T
   x_ = reshape(x,numVox,numContr)
-  y_ = reshape(y,numVox,numContr,numChan)
+  y_ = reshape(y,numVox,numChan,numContr)
 
-  @assert size(smaps) == (size(y_,1), size(y_,3))
+  @assert size(smaps) == (size(y_,1), size(y_,2))
 
   @inbounds for i ∈ CartesianIndices(y_)
-    y_[i] = x_[i[1],i[2]] * smaps[i[1],i[3]]
+    y_[i] = x_[i[1],i[3]] * smaps[i[1],i[2]]
   end
   return y
 end
 
 function ctprod_smap!(y::AbstractVector{T}, smapsC::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numChan, numContr=1) where T
-  x_ = reshape(x,numVox,numContr,numChan)
+  x_ = reshape(x,numVox,numChan,numContr)
   y_ = reshape(y,numVox,numContr)
 
-  @assert size(smapsC) == (size(x_,1), size(x_,3))
+  @assert size(smapsC) == (size(x_,1), size(x_,2))
 
   y_ .= 0
   @inbounds for i ∈ CartesianIndices(x_)
-    y_[i[1],i[2]] += x_[i] * smapsC[i[1],i[3]]
+    y_[i[1],i[3]] += x_[i] * smapsC[i[1],i[2]]
   end
   return y
 end

diff --git a/MRIOperators/src/SubspaceOp.jl b/MRIOperators/src/SubspaceOp.jl
@@ -0,0 +1,37 @@
+export SubspaceOp
+
+function prod_subspace!(y::Vector{T}, basis::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numContr, numBasis) where T
+  x_ = reshape(x,numVox,numBasis)
+  y_ = reshape(y,numVox,numContr)
+
+  y_ .= x_ * basis'
+  return y
+end
+
+function ctprod_subspace!(y::Vector{T},basis::AbstractMatrix{T}, x::AbstractVector{T}, numVox, numContr, numBasis) where T
+  x_ = reshape(x,numVox,numContr)
+  y_ = reshape(y,numVox,numBasis)
+
+  y_ .= x_ * basis
+  return y
+end
+
+
+"""
+  SubspaceOp(basis::AbstractMatrix{T},shape::NTuple{D,Int64},numContr) where {T,D}
+
+basis correspond to svd_object.V cropped to a certain level
+
+basis type should correspond to the type of the rawdata
+"""
+function SubspaceOp(basis::AbstractMatrix{T},shape::NTuple{D,Int64},numContr) where {T,D}
+    numVox = prod(shape)
+    numBasis = size(basis,2)
+
+    return LinearOperator{T}(numVox*numContr, numVox*numBasis, false, false,
+                         (res,x) -> prod_subspace!(res,basis,x,numVox,numContr,numBasis),
+                         nothing,
+                         (res,x) -> ctprod_subspace!(res,basis,x,numVox,numContr,numBasis))
+end
+
+