From 891e5df002b0118ff3bcbfdfd25e324cca3b6930 Mon Sep 17 00:00:00 2001
From: migrosser <mi.grosser@uke.de>
Date: Wed, 3 Jun 2020 14:23:08 +0200
Subject: [PATCH] add CGNR-based regridding method

---
 src/Tools/Regridding.jl    | 40 +++++++++++++++++++++++++++++++
 src/Tools/Tools.jl         |  1 +
 test/testReconstruction.jl | 48 ++++++++++++++++++++++++++++++++++++--
 3 files changed, 87 insertions(+), 2 deletions(-)
 create mode 100644 src/Tools/Regridding.jl

diff --git a/src/Tools/Regridding.jl b/src/Tools/Regridding.jl
new file mode 100644
index 00000000..4b8f4cf2
--- /dev/null
+++ b/src/Tools/Regridding.jl
@@ -0,0 +1,40 @@
+export regrid2d
+
+"""
+    `regrid2d(acqData::AcquisitionData, kspaceSize::NTuple{2,Int64}
+              ; cgnr_iter::Int64=3, correctionMap::Array{ComplexF64}=ComplexF64[]) where D`
+
+Regrid non-cartesian k-space data in `acqData` to a cartesian grid of size `kspaceSize`.
+Uses the CGNR method to invert the non-cartesian Fourier encoding.
+
+# Arguments:
+* `acqData::AcquisitionData`        - AcquisitionData
+* `ksize::NTuple{2,Int64}`          - size of the cartesian k-space grid
+* `cgnr_iter::Int64=3`              - number of CGNR iterations
+* `correctionMao::Array{ComplexF64}`- relaxation/b0 map
+"""
+function regrid2d(acqData::AcquisitionData, kspaceSize::NTuple{2,Int64}; cgnr_iter::Int64=3, correctionMap::Array{ComplexF64}=ComplexF64[]) where D
+  dcf = samplingDensity(acqData,kspaceSize)
+
+  nx,ny = kspaceSize
+  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
+  Ireco = zeros(ComplexF64, prod(kspaceSize), numSl, numContr, numChan)
+
+  kdata_cart = [zeros(ComplexF64,nx*ny,numChan) for j=1:numContr, k=1:numSl, rep=1:1]
+  F = sqrt(prod(kspaceSize))*FFTOp(ComplexF64, kspaceSize)
+  img = zeros(ComplexF64,nx*ny)
+  for k = 1:numSl
+    E = encodingOps2d_simple(acqData, kspaceSize, slice=k, correctionMap=correctionMap)
+    for j = 1:numContr
+      W = WeightingOp(dcf[j])
+      solver = RegularizedLeastSquares.CGNR(W*E[j], iterations=cgnr_iter)
+      for i = 1:numChan
+        kdata = kData(acqData,j,i,k).* dcf[j]
+        img .= solve(solver, kdata)
+        kdata_cart[j,k,1][:,i] .= F*img
+      end
+    end
+  end
+
+  return AcquisitionData(CartesianTrajectory(ny,nx), kdata_cart, encodingSize=[nx,ny,1], fov=acqData.fov)
+end
diff --git a/src/Tools/Tools.jl b/src/Tools/Tools.jl
index 4d609c16..a2443125 100644
--- a/src/Tools/Tools.jl
+++ b/src/Tools/Tools.jl
@@ -1,3 +1,4 @@
 include("Shutter.jl")
 include("ErrorMeasures.jl")
 include("CoilSensitivity.jl")
+include("Regridding.jl")
diff --git a/test/testReconstruction.jl b/test/testReconstruction.jl
index 8b0912fa..7159ca0f 100644
--- a/test/testReconstruction.jl
+++ b/test/testReconstruction.jl
@@ -209,8 +209,8 @@ function testSENSEReco(N = 64)
   params[:reco] = "multiCoil" #"standard"
   params[:reconSize] = (N,N)
   params[:regularization] = "L2"
-  params[:iterations] = 3
-  params[:solver] = "admm"
+  params[:iterations] = 50
+  params[:solver] = "cgnr"
   params[:senseMaps] = reshape(coilsens, N*N, numCoils, 1)
 
   Ireco = reconstruction(acqData, params)
@@ -359,6 +359,49 @@ function testCSSenseReco3d(N=64)
   @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 5e-2
 end
 
+function testRegridding(N=64)
+
+  # image
+  x = shepp_logan(N)
+
+  # simulation
+  params = Dict{Symbol, Any}()
+  params[:simulation] = "fast"
+  params[:trajName] = "Radial"
+  params[:numProfiles] = round(Int64,1.25*N)
+  params[:numSamplingPerProfile] = round(Int64,1.25*N)
+
+  acqDataRad = simulation(x, params)
+
+  params[:trajName] = "Cartesian"
+  params[:numProfiles] = N
+  params[:numSamplingPerProfile] = N
+  acqDataCart = simulation(x, params)
+
+  # regridding
+  acqDataReg = regrid2d(acqDataRad,(N,N))
+  # cartesian reconstruction
+  params[:reco] = "standard"
+  params[:reconSize] = (N,N)
+  params[:solver] = "cgnr"
+  params[:regularization] = "L2"
+  params[:λ] = 0.0
+  params[:iterations] = 3
+
+  x_reg = collect(reshape(reconstruction(acqDataReg, params),N,N))
+  circularShutter!(x_reg)
+  x_rad = collect(reshape(reconstruction(acqDataRad, params),N,N))
+  circularShutter!(x_rad)
+
+  k_reg = reshape(acqDataReg.kdata[1],N,N)
+  k_cart = reshape(acqDataCart.kdata[1],N,N)
+  circularShutter!(k_reg)
+  circularShutter!(k_cart)
+
+  @test (norm(vec(k_cart)-vec(k_reg))/norm(vec(k_cart))) < 1e-1
+  @test (norm(vec(x_rad).-vec(x_reg))/norm(vec(x_rad))) < 1e-2
+end
+
 function testReco(N=32)
   @testset "Reconstructions" begin
     testGriddingReco()
@@ -378,6 +421,7 @@ function testReco(N=32)
     testSENSEReco()
     testOffresonanceSENSEReco()
     testDirectRecoMultiEcho()
+    testRegridding()
   end
 end