diff --git a/Project.toml b/Project.toml index 918e3bb..a591ec8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NFFT" uuid = "efe261a4-0d2b-5849-be55-fc731d526b0d" authors = ["Tobias Knopp "] -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" diff --git a/src/precomputation.jl b/src/precomputation.jl index 339d532..5f592bb 100644 --- a/src/precomputation.jl +++ b/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 diff --git a/test/issues.jl b/test/issues.jl new file mode 100644 index 0000000..8789d8a --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 615e051..8a7ba29 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ using NFFTTools Random.seed!(123) +include("issues.jl") include("accuracy.jl") include("constructors.jl") include("performance.jl")