fix ndft and finalize accuracy comparison

JuliaMath · Dec 30, 2021 · e5c18d5 · e5c18d5
1 parent fd24ee6
commit e5c18d5
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 9 deletions.
diff --git a/src/NDFT.jl b/src/NDFT.jl
@@ -2,18 +2,18 @@
 ### ndft functions ###
 
 """
-    ndft!(plan::NFFTPlan{D}, g::AbstractArray{Tg,D}, f::AbstractArray{T,D})
+    ndft!(plan::NFFTPlan{D}, g::AbstractVector{Tg}, f::AbstractArray{T,D})
 
 Compute NDFT of input array `f`
 and store result in pre-allocated output array `g`.
 Both arrays must have the same size compatible with the NFFT `plan`.
 
 """
-function ndft!(plan::NFFTPlan{D}, g::AbstractArray{Tg,D}, f::AbstractArray{T,D}) where {D,T,Tg}
+function ndft!(plan::NFFTPlan{D}, g::AbstractArray{Tg}, f::AbstractArray{T,D}) where {D,T,Tg}
 
     plan.N == size(f) ||
         throw(DimensionMismatch("Data f is not consistent with NFFTPlan"))
-    plan.N == size(g) ||
+    plan.M == length(g) ||
         throw(DimensionMismatch("Output g is inconsistent with NFFTPlan"))
 
     g .= zero(Tg)
@@ -41,7 +41,7 @@ end
 Non pre-allocated versions of NDFT; see `ndft!`.
 """
 ndft(plan::NFFTPlan{D}, f::AbstractArray{T,D}) where {D,T} =
-    ndft!(plan, similar(f), f)
+    ndft!(plan, similar(f,plan.M), f)
 
 ndft(x, f::AbstractArray, rest...; kwargs...) =
     ndft(NFFTPlan(x, size(f), rest...; kwargs...), f)

diff --git a/test/comparison.jl b/test/comparison.jl
@@ -26,24 +26,28 @@ function nfft_accuracy_comparison()
           p = plan_nfft(x, NN, m, sigma; precompute=NFFT.FULL)
           f = ndft_adjoint(p, fHat)
           fApprox = nfft_adjoint(p, fHat)
-          etrafo = norm(f[:] - fApprox[:]) / norm(f[:])
+          eadjoint = norm(f[:] - fApprox[:]) / norm(f[:])
 
           gHat = ndft(p, f)
           gHatApprox = nfft(p, f)
-          eadjoint = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
+          etrafo = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
 
           push!(df, ("NFFT.jl", D, M, N[D], m, sigma, etrafo, eadjoint))
 
           pnfft3 = NFFT3.NFFT(NN, M, Int32.(p.n), m) 
           pnfft3.x = (D==1) ? vec(x) : x
           pnfft3.f = fHat
           NFFT3.nfft_adjoint(pnfft3)
-          fApprox = pnfft3.fhat
-          etrafo = norm(f[:] - fApprox[:]) / norm(f[:])
+          fApprox = reshape(pnfft3.fhat,reverse(NN)...)
+          # switch from column major to row major format
+          fApprox = (D==1) ? vec(fApprox) : vec(collect(permutedims(fApprox,D:-1:1)))   
+          eadjoint = norm(f[:] - fApprox[:]) / norm(f[:])
 
+          # switch from column major to row major format
+          pnfft3.fhat = (D==1) ? vec(f) : vec(collect(permutedims(f,D:-1:1))) 
           NFFT3.nfft_trafo(pnfft3)
           gHatApprox = pnfft3.f
-          eadjoint = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
+          etrafo = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
 
           push!(df, ("NFFT3", D, M, N[D], m, sigma, etrafo, eadjoint))
       end