@@ -183,125 +183,5 @@ function NFFT.nfft(k::AbstractMatrix, f::Reactant.AnyTracedRArray, args...; kwar
183183 return p * f
184184end
185185
186- function NFFT. initParams(
187- k:: AbstractMatrix{T} ,
188- N:: NTuple{D, Int} ,
189- dims:: Union{Integer, UnitRange{Int64}} = 1 : D;
190- kargs... ,
191- ) where {D, T}
192- # convert dims to a unit range
193- dims_ = (typeof(dims) <: Integer ) ? (dims: dims) : dims
194-
195- params = NFFTParams{T, D}(; kargs... )
196- m, σ, reltol = accuracyParams(; kargs... )
197- params. m = m
198- params. σ = σ
199- params. reltol = reltol
200-
201- # Taken from NFFT3
202- m2K = [1 , 3 , 7 , 9 , 14 , 17 , 20 , 23 , 24 ]
203- K = m2K[min(m + 1 , length(m2K))]
204- params. LUTSize = 2 ^ (K) * (m) # ensure that LUTSize is dividable by (m)
205-
206- if length(dims_) != size(k, 1 )
207- throw(ArgumentError(" Nodes x have dimension $(size(k, 1 )) != $(length(dims_)) " ))
208- end
209-
210- doTrafo = ntuple(d -> d ∈ dims_, Val(D))
211-
212- Ñ = ntuple(
213- d ->
214- doTrafo[d] ? (ceil(Int, params. σ * N[d]) ÷ 2 ) * 2 : # ensure that n is an even integer
215- N[d], Val(D)
216- )
217-
218- params. σ = Ñ[dims_[1 ]] / N[dims_[1 ]]
219-
220- # params.blockSize = ntuple(d-> Ñ[d] , D) # just one block
221- if haskey(kargs, :blockSize)
222- params. blockSize = kargs[:blockSize]
223- else
224- params. blockSize = ntuple(d -> NFFT. _blockSize(Ñ, d), Val(D))
225- end
226-
227- J = size(k, 2 )
228-
229- # calculate output size
230- NOut = Int[]
231- Mtaken = false
232- ntuple(Val(D)) do d
233- if ! doTrafo[d]
234- return N[d]
235- elseif ! Mtaken
236- return J
237- Mtaken = true
238- end
239- end
240- for d in 1 : D
241- if ! doTrafo[d]
242- push!(NOut, N[d])
243- elseif ! Mtaken
244- push!(NOut, J)
245- Mtaken = true
246- end
247- end
248- # Sort nodes in lexicographic way
249- if params. sortNodes
250- k .= sortslices(k; dims = 2 )
251- end
252- return params, N, Tuple(NOut), J, Ñ, dims_
253- end
254-
255- function NFFT. precomputation(k:: AbstractVecOrMat , N:: NTuple{D, Int} , Ñ, params) where {D}
256- m = params. m
257- σ = params. σ
258- window = params. window
259- LUTSize = params. LUTSize
260- precompute = params. precompute
261-
262- win, win_hat = getWindow(window) # highly type instable. But what should be do
263- J = size(k, 2 )
264-
265- windowHatInvLUT_ = Vector{Vector{T}}(undef, D)
266- precomputeWindowHatInvLUT(windowHatInvLUT_, win_hat, N, Ñ, m, σ, T)
267-
268- if params. storeDeconvolutionIdx
269- windowHatInvLUT = Vector{Vector{T}}(undef, 1 )
270- windowHatInvLUT[1 ], deconvolveIdx = precompWindowHatInvLUT(
271- params, N, Ñ, windowHatInvLUT_
272- )
273- else
274- windowHatInvLUT = windowHatInvLUT_
275- deconvolveIdx = Array{Int64, 1 }(undef, 0 )
276- end
277-
278- if precompute == LINEAR
279- windowLinInterp = precomputeLinInterp(win, m, σ, LUTSize, T)
280- windowPolyInterp = Matrix{T}(undef, 0 , 0 )
281- B = sparse([], [], T[])
282- elseif precompute == POLYNOMIAL
283- windowLinInterp = Vector{T}(undef, 0 )
284- windowPolyInterp = precomputePolyInterp(win, m, σ, T)
285- B = sparse([], [], T[])
286- elseif precompute == FULL
287- windowLinInterp = Vector{T}(undef, 0 )
288- windowPolyInterp = Matrix{T}(undef, 0 , 0 )
289- B = precomputeB(win, k, N, Ñ, m, J, σ, LUTSize, T)
290- # windowLinInterp = precomputeLinInterp(win, windowLinInterp, Ñ, m, σ, LUTSize, T) # These versions are for debugging
291- # B = precomputeB(windowLinInterp, k, N, Ñ, m, J, σ, LUTSize, T)
292- elseif precompute == TENSOR
293- windowLinInterp = Vector{T}(undef, 0 )
294- windowPolyInterp = Matrix{T}(undef, 0 , 0 )
295- B = sparse([], [], T[])
296- else
297- windowLinInterp = Vector{T}(undef, 0 )
298- windowPolyInterp = Matrix{T}(undef, 0 , 0 )
299- B = sparse([], [], T[])
300- error(" precompute = $precompute not supported by NFFT.jl!" )
301- end
302-
303- return (windowLinInterp, windowPolyInterp, windowHatInvLUT, deconvolveIdx, B)
304- end
305-
306186
307187end
0 commit comments