add CGNR-based regridding method

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

src/Tools/Tools.jl

@@ -1,3 +1,4 @@
 include("Shutter.jl")
 include("ErrorMeasures.jl")
 include("CoilSensitivity.jl")
+include("Regridding.jl")

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