Merge pull request #44 from IvoMaatman/master

Fix for lower triangle elements of kernel

src/NFFTNormalOpBasisFunc.jl

@@ -100,11 +100,9 @@ function calculateToeplitzKernelBasis(img_shape_os, trj::AbstractVector{<:Abstra
     @assert all(kmask_indcs .<= prod(img_shape_os))
 
     Ncoeff = size(U, 2)
-    Nt = size(U,1)
 
     λ  = Array{Complex{T}}(undef, img_shape_os)
     λ2 = similar(λ)
-    λ3 = similar(λ)
     Λ  = Array{Complex{T}}(undef, Ncoeff, Ncoeff, length(kmask_indcs))
 
     trj_l = [size(trj[it],2) for it in eachindex(trj)]
@@ -113,6 +111,9 @@ function calculateToeplitzKernelBasis(img_shape_os, trj::AbstractVector{<:Abstra
     fftplan  = plan_fft(λ; flags = FFTW.MEASURE, num_threads=Threads.nthreads())
     nfftplan = plan_nfft(reduce(hcat, trj), img_shape_os; precompute = TENSOR, blocking = true, fftflags = FFTW.MEASURE, m=5, σ=2)
 
+    # Evaluating only the upper triangular matrix assumes that the PSF from the rightmost voxel to the leftmost voxel is the adjoint of the PSF in the opposite direction. 
+    # For the outmost voxel, this is not correct, but the resulting images are virtually identical in our test cases. 
+    # To avoid this error, remove the `if ic2 >= ic1` and the `Λ[ic2,ic1,it] = conj.(λ[kmask_indcs[it]])` statements at the cost of computation time.  
     for ic2 ∈ axes(Λ, 2), ic1 ∈ axes(Λ, 1)
         if ic2 >= ic1 # eval. only upper triangular matrix
             t = @elapsed begin
@@ -125,12 +126,10 @@ function calculateToeplitzKernelBasis(img_shape_os, trj::AbstractVector{<:Abstra
                 mul!(λ, adjoint(nfftplan), vec(S))
                 fftshift!(λ2, λ)
                 mul!(λ, fftplan, λ2)
-                λ2 .= conj.(λ2)
-                mul!(λ3, fftplan, λ2)
 
                 Threads.@threads for it ∈ eachindex(kmask_indcs)
-                    @inbounds Λ[ic2,ic1,it] = λ3[kmask_indcs[it]]
-                    @inbounds Λ[ic1,ic2,it] =  λ[kmask_indcs[it]]
+                    @inbounds Λ[ic2,ic1,it] = conj.(λ[kmask_indcs[it]]) 
+                    @inbounds Λ[ic1,ic2,it] =       λ[kmask_indcs[it]]
                 end
             end
             verbose && println("ic = ($ic1, $ic2): t = $t s"); flush(stdout)

test/data_removal.jl

@@ -40,7 +40,6 @@ for i ∈ CartesianIndices(@view cmaps[:,:,1])
 end
 cmaps = [cmaps[:,:,ic] for ic=1:Ncoil]
 
-
 ## set up trajectory
 α_g = 2π / (1+√5)
 phi = Float32.(α_g * (1:Nt*Ncyc))
@@ -95,20 +94,6 @@ b = calculateBackProjection(data, trj, cmaps; U=U)
 ## construct forward operator
 A = NFFTNormalOp((Nx,Nx), trj, U, cmaps=cmaps)
 
-## test forward operator
-λ = zeros(Complex{T}, Nc, Nc, 2Nx*2Nx)
-for i ∈ eachindex(A.prod!.A.kmask_indcs)
-    λ[:,:,A.prod!.A.kmask_indcs[i]] .= A.prod!.A.Λ[:,:,i]
-end
-λ = reshape(λ, Nc, Nc, 2Nx, 2Nx)
-
-for i = 1:Nc, j = 1:Nc
-    l1 = conj.(λ[i,j,:,:])
-    l2 = λ[j,i,:,:]
-    l2 = conj.(fft(conj.(ifft(l2))))
-    @test l1 ≈ l2 rtol = 1e-4
-end
-
 ## reconstruct
 xr = cg(A, vec(b), maxiter=20)
 xr = reshape(xr, Nx, Nx, Nc)

test/reconstruct_radial.jl

@@ -75,20 +75,6 @@ b = calculateBackProjection(data, trj, cmaps; U)
 ## construct forward operator
 A = NFFTNormalOp((Nx,Nx), trj, U, cmaps=cmaps)
 
-## test forward operator
-λ = zeros(Complex{T}, Nc, Nc, 2Nx*2Nx)
-for i ∈ eachindex(A.prod!.A.kmask_indcs)
-    λ[:,:,A.prod!.A.kmask_indcs[i]] .= A.prod!.A.Λ[:,:,i]
-end
-λ = reshape(λ, Nc, Nc, 2Nx, 2Nx)
-
-for i = 1:Nc, j = 1:Nc
-    l1 = conj.(λ[i,j,:,:])
-    l2 = λ[j,i,:,:]
-    l2 = conj.(fft(conj.(ifft(l2))))
-    @test l1 ≈ l2 rtol = 1e-4
-end
-
 ## reconstruct
 xr = cg(A, vec(b), maxiter=20)
 xr = reshape(xr, Nx, Nx, Nc)