Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/tknopp/NFFT.jl into Toepl…
Browse files Browse the repository at this point in the history 
…itzGram
  • Loading branch information
@JakobAsslaender
JakobAsslaender committed Dec 31, 2021
2 parents f6c7ea4 + e5c18d5 commit 11dd913
Show file tree
Hide file tree
Showing 17 changed files with 690 additions and 229 deletions.
12 changes: 11 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,20 +9,30 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
Graphics = "a2bd30eb-e257-5431-a919-1863eab51364"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
ThreadedSparseCSR = "39591389-4884-462a-89e7-9edfce40b6eb"

[compat]
CUDA = "1.3.3, 2.3, 3.1"
DataFrames = "1.3.1"
FFTW = "0.2, 1"
Graphics = "0.4, 1.0"
LoopVectorization = "0.12"
Polyester = "0.5.5, 0.6"
SpecialFunctions = "0.8, 0.10, 1, 2"
ThreadedSparseCSR = "0.1.1"
julia = "1.6"

[extras]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
NFFT3 = "53104703-03e8-40a5-ab01-812303a44cae"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]
test = ["Test", "BenchmarkTools", "NFFT3", "DataFrames"]
35 changes: 0 additions & 35 deletions appveyor.yml
View file

This file was deleted.

28 changes: 6 additions & 22 deletions src/Abstract.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,17 +26,17 @@ For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function nfft(p::AbstractNFFTPlan{D,0,T}, f::AbstractArray{U,D}, args...) where {D,T,U}
function nfft(p::AbstractNFFTPlan{D,0,T}, f::AbstractArray{U,D}, args...; kargs...) where {D,T,U}
fHat = similar(f,Complex{T}, p.M)
nfft!(p, f, fHat, args...)
nfft!(p, f, fHat, args...; kargs...)
return fHat
end

function nfft(p::AbstractNFFTPlan{D,DIM,T}, f::AbstractArray{U,D}, args...) where {D,DIM,T,U}
function nfft(p::AbstractNFFTPlan{D,DIM,T}, f::AbstractArray{U,D}, args...; kargs...) where {D,DIM,T,U}
sz = [p.N...]
sz[DIM] = p.M
fHat = similar(f, Complex{T}, Tuple(sz))
nfft!(p, f, fHat, args...)
nfft!(p, f, fHat, args...; kargs...)
return fHat
end

Expand All @@ -51,25 +51,9 @@ For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function nfft_adjoint(p::AbstractNFFTPlan{D,DIM,T}, fHat::AbstractArray{U}, args...) where {D,DIM,T,U}
function nfft_adjoint(p::AbstractNFFTPlan{D,DIM,T}, fHat::AbstractArray{U}, args...; kargs...) where {D,DIM,T,U}
f = similar(fHat, Complex{T}, p.N)
nfft_adjoint!(p, fHat, f, args...)
nfft_adjoint!(p, fHat, f, args...; kargs...)
return f
end

@generated function consistencyCheck(p::AbstractNFFTPlan{D,DIM,T}, f::AbstractArray{U,D},
fHat::AbstractArray{Y}) where {D,DIM,T,U,Y}
quote
M = numFourierSamples(p)
N = size(p)
if $DIM == 0
fHat_test = (M == length(fHat))
elseif $DIM > 0
fHat_test = @nall $D d -> ( d == $DIM ? size(fHat,d) == M : size(fHat,d) == N[d] )
end

if N != size(f) || !fHat_test
throw(DimensionMismatch("Data is not consistent with NFFTPlan"))
end
end
end
80 changes: 52 additions & 28 deletions src/CpuNFFT.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,14 +75,19 @@ end
# constructors
##############
function NFFTPlan(x::Matrix{T}, N::NTuple{D,Int}, m = 4, sigma = 2.0,
window = :kaiser_bessel, K = 2000;
precompute::PrecomputeFlags = LUT, kwargs...) where {D,T}
window=:kaiser_bessel, K=2000;
precompute::PrecomputeFlags=LUT, sortNodes=false, kwargs...) where {D,T}

if D != size(x, 1)
throw(ArgumentError("The number of dimensions of x and N do not match"))
if D != size(x,1)
throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $D!"))
end

n = ntuple(d -> round(Int, sigma * N[d]), D)
if any(isodd.(N))
throw(ArgumentError("N = $N needs to consist of even integers!"))
end

n = ntuple(d->(ceil(Int,sigma*N[d])÷2)*2, D) # ensure that n is an even integer
sigma = n[1] / N[1]

tmpVec = Array{Complex{T}}(undef, n)

Expand All @@ -91,19 +96,29 @@ function NFFTPlan(x::Matrix{T}, N::NTuple{D,Int}, m = 4, sigma = 2.0,
FP = plan_fft!(tmpVec; kwargs...)
BP = plan_bfft!(tmpVec; kwargs...)

# Sort nodes in lexicographic way
if sortNodes
x .= sortslices(x,dims=2)
end

windowLUT, windowHatInvLUT, B = calcLookUpTable(x, N, n, m, sigma, window, K; precompute=precompute)

NFFTPlan{D,0,T}(N, M, x, m, sigma, n, K, windowLUT, windowHatInvLUT, FP, BP, tmpVec, B)
end

function NFFTPlan!(p::AbstractNFFTPlan{D,DIM,T}, x::Matrix{T}, window = :kaiser_bessel) where {D,DIM,T}
function NFFTPlan!(p::AbstractNFFTPlan{D,DIM,T}, x::Matrix{T}, window = :kaiser_bessel; sortNodes=false) where {D,DIM,T}

if isempty(p.B)
precompute = LUT
else
precompute = FULL
end

# Sort nodes in lexicographic way
if sortNodes
x .= sortslices(x,dims=2)
end

windowLUT, windowHatInvLUT, B = calcLookUpTable(x, p.N, p.n, p.m, p.sigma, window, p.K; precompute=precompute)

p.M = size(x, 2)
Expand All @@ -125,9 +140,15 @@ A 1D plan that is applied along dimension `d` of a `D` dimensional array of size
It takes as optional keywords all the keywords supported by `plan_fft` function (like
`flags` and `timelimit`). See documentation of `plan_fft` for reference.
"""
function NFFTPlan(x::AbstractVector{T}, dim::Integer, N::NTuple{D,Int64}, m = 4, sigma = 2.0, window = :kaiser_bessel, K = 2000; kwargs...) where {D,T}
function NFFTPlan(x::AbstractVector{T}, dim::Integer, N::NTuple{D,Int64}, m=4,
sigma=2.0, window=:kaiser_bessel, K=2000; kwargs...) where {D,T}

n = ntuple(d -> round(Int, sigma * N[d]), D)
if any(isodd.(N))
throw(ArgumentError("N = $N needs to consist of even integers!"))
end

n = ntuple(d->(ceil(Int,sigma*N[d])÷2)*2, D) # ensure that n is an even integer
sigma = n[1] / N[1]

sz = [N...]
sz[dim] = n[dim]
Expand Down Expand Up @@ -162,6 +183,7 @@ numFourierSamples(p::NFFTPlan) = p.M
function calcLookUpTable(x::Union{Matrix{T},Vector{T}}, N::NTuple{D,Int}, n, m = 4, sigma = 2.0, window = :kaiser_bessel, K = 2000; precompute::PrecomputeFlags = LUT) where {T,D}

win, win_hat = getWindow(window)
M = size(x, 2)

windowLUT = Vector{Vector{T}}(undef, D)
windowHatInvLUT = Vector{Vector{T}}(undef, D)
Expand All @@ -173,25 +195,15 @@ function calcLookUpTable(x::Union{Matrix{T},Vector{T}}, N::NTuple{D,Int}, n, m =
end

if precompute == LUT
Z = round(Int, 3 * K / 2)
for d = 1:D
windowLUT[d] = Vector{T}(undef, Z)
@batch for l = 1:Z
y = ((l - 1) / (K - 1)) * m / n[d]
windowLUT[d][l] = win(y, n[d], m, sigma)
end
end

B = sparse([], [], Float64[])
precomputeLUT(win, windowLUT, n, m, sigma, K, T)
B = sparse([],[],T[])
elseif precompute == FULL
B = precomputeB(win, x, N, n, m, M, sigma, K, T)
elseif precompute == FULL_LUT
precomputeLUT(win, windowLUT, n, m, sigma, K, T)
B = precomputeB(windowLUT, x, N, n, m, M, sigma, K, T)
else
if typeof(x) <: Vector
error("For directional NFFTs, only `precompute = LUT` is supported")
end

M = size(x, 2)
U1 = ntuple(d -> (d == 1) ? 1 : (2 * m + 1)^(d - 1), D)
U2 = ntuple(d -> (d == 1) ? 1 : prod(n[1:(d-1)]), D)
B = precomputeB(win, x, n, m, M, sigma, T, U1, U2)
error("precompute = $precompute not supported by NFFT.jl!")
end
return (windowLUT, windowHatInvLUT, B)
end
Expand All @@ -208,7 +220,8 @@ Calculate the NFFT of `f` with plan `p` and store the result in `fHat`.
Both `f` and `fHat` must be complex arrays.
"""
function nfft!(p::NFFTPlan{D,DIM,T}, f::AbstractArray, fHat::StridedArray, verbose = false) where {D,DIM,T}
function nfft!(p::NFFTPlan{D,DIM,T}, f::AbstractArray, fHat::StridedArray;
verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {D,DIM,T}
consistencyCheck(p, f, fHat)

fill!(p.tmpVec, zero(Complex{T}))
Expand All @@ -222,6 +235,11 @@ function nfft!(p::NFFTPlan{D,DIM,T}, f::AbstractArray, fHat::StridedArray, verbo
if verbose
@info "Timing: apod=$t1 fft=$t2 conv=$t3"
end
if timing != nothing
timing.conv = t3
timing.fft = t2
timing.apod = t1
end
return fHat
end

Expand All @@ -232,7 +250,8 @@ Calculate the adjoint NFFT of `fHat` and store the result in `f`.
Both `f` and `fHat` must be complex arrays.
"""
function nfft_adjoint!(p::NFFTPlan, fHat::AbstractArray, f::StridedArray, verbose = false)
function nfft_adjoint!(p::NFFTPlan, fHat::AbstractArray, f::StridedArray;
verbose=false, timing::Union{Nothing,TimingStats} = nothing)
consistencyCheck(p, f, fHat)

t1 = @elapsed @inbounds convolve_adjoint!(p, fHat, p.tmpVec)
Expand All @@ -245,5 +264,10 @@ function nfft_adjoint!(p::NFFTPlan, fHat::AbstractArray, f::StridedArray, verbos
if verbose
@info "Timing: conv=$t1 fft=$t2 apod=$t3"
end
if timing != nothing
timing.conv_adjoint = t1
timing.fft_adjoint = t2
timing.apod_adjoint = t3
end
return f
end
11 changes: 10 additions & 1 deletion src/CuNFFT.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,16 @@ end
function CuNFFTPlan(x::Matrix{T}, N::NTuple{D,Int}, m=4, sigma=2.0,
window=:kaiser_bessel, K=2000; kwargs...) where {T,D}

n = ntuple(d->round(Int,sigma*N[d]), D)
if D != size(x,1)
throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $D!"))
end

if any(isodd.(N))
throw(ArgumentError("N = $N needs to consist of even integers!"))
end

n = ntuple(d->(ceil(Int,sigma*N[d])÷2)*2, D) # ensure that n is an even integer
sigma = n[1] / N[1]

tmpVec = CuArray(zeros(Complex{T},n))
tmpVec2 = CuArray(zeros(Complex{T},N))
Expand Down
8 changes: 4 additions & 4 deletions src/NDFT.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 11dd913

Please sign in to comment.