quick fix for #106

Author: tknopp

Project.toml

@@ -1,7 +1,7 @@
 name = "NFFT"
 uuid = "efe261a4-0d2b-5849-be55-fc731d526b0d"
 authors = ["Tobias Knopp <tobias@knoppweb.de>"]
-version = "0.13.1"
+version = "0.13.2"
 
 [deps]
 AbstractNFFTs = "7f219486-4aa7-41d6-80a7-e08ef20ceed7"
@@ -21,7 +21,7 @@ SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 [compat]
 AbstractNFFTs = "0.8"
 BasicInterpolators = "0.6.5"
-DataFrames = "1.3.1"
+DataFrames = "1.3.1, 1.4.1"
 FFTW = "1.5"
 FLoops = "0.2"
 Reexport = "1.0"

src/precomputation.jl

@@ -180,26 +180,32 @@ function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::Nothing,
   return (y, tmpWin)
 end
 
-@generated function shiftedWindowEntries(winLin::Vector, idx, scale, d, L::Val{Z}) where {Z}
+@generated function shiftedWindowEntries(winLin::Vector, idx::T, scale, d, L::Val{Z}) where {T,Z}
   quote
-    idxL = floor(Int,idx)
-    idxInt = Int(idxL)
-    α = ( idx-idxL )
-
-    tmpWin = @ntuple $(Z) l -> begin
-      # Uncommented code: This is the version where we pull in l into the abs.
-      # We pulled this out of the iteration.
-      # idx = abs((kscale - (l-1)  - off)*LUTSize)/(m)
-
-      # The second +1 is because Julia has 1-based indexing
-      # The first +1 is part of the index calculation and needs(!)
-      # to be part of the abs. The abs is shifting to the positive interval
-      # and this +1 matches the `floor` above, which rounds down. In turn
-      # for positive and negative indices a different neighbour is calculated
-      idxInt1 = abs( idxInt - (l-1)*scale ) +1
-      idxInt2 = abs( idxInt - (l-1)*scale +1) +1
-
-      (winLin[idxInt1] + α * (winLin[idxInt2] - winLin[idxInt1]))
+    idxInt = floor(Int,idx)
+    α = ( idx-idxInt )
+
+    if α != zero(T) 
+      tmpWin = @ntuple $(Z) l -> begin
+        # Uncommented code: This is the version where we pull in l into the abs.
+        # We pulled this out of the iteration.
+        # idx = abs((kscale - (l-1)  - off)*LUTSize)/(m)
+
+        # The second +1 is because Julia has 1-based indexing
+        # The first +1 is part of the index calculation and needs(!)
+        # to be part of the abs. The abs is shifting to the positive interval
+        # and this +1 matches the `floor` above, which rounds down. In turn
+        # for positive and negative indices a different neighbour is calculated
+        idxInt1 = abs( idxInt - (l-1)*scale ) +1
+        idxInt2 = abs( idxInt - (l-1)*scale +1) +1
+
+        (winLin[idxInt1] + α * (winLin[idxInt2] - winLin[idxInt1]))
+      end
+    else
+      tmpWin = @ntuple $(Z) l -> begin
+        idxInt1 = abs( idxInt - (l-1)*scale ) + 1
+        winLin[idxInt1]
+      end
     end
     return tmpWin
   end

test/issues.jl

@@ -0,0 +1,18 @@
+@testset "Issues Encountered during Development" begin
+
+  # https://github.com/JuliaMath/NFFT.jl/issues/106
+  # The issue was an out of bounce in the linear interpolation code.
+  T = Float32
+  Nz = 120
+  img_shape_os = (2Nz,)
+  λ = Array{Complex{T}}(undef, img_shape_os)
+  
+  trj = Matrix{T}(undef, 1, 2)
+  nfftplan = plan_nfft(trj, img_shape_os; precompute = LINEAR, blocking = false, fftflags = FFTW.ESTIMATE)
+  
+  trj[1,:] .= 0.008333333 # throws error
+  nfftplan = plan_nfft(trj, img_shape_os; precompute = LINEAR, blocking = false, fftflags = FFTW.ESTIMATE)
+  mul!(λ, adjoint(nfftplan), ones(Complex{T}, size(trj,2)))
+  
+
+end

test/runtests.jl

@@ -7,6 +7,7 @@ using NFFTTools
 
 Random.seed!(123)
 
+include("issues.jl")
 include("accuracy.jl")
 include("constructors.jl")
 include("performance.jl")