forked from JuliaMath/NFFT.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
testToeplitz.jl
85 lines (67 loc) · 2.61 KB
/
testToeplitz.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
using Test
using NFFT
using FFTW
using BenchmarkTools
@testset "Toeplitz" begin
Nx = 32
## calculateToeplitzKernel vs calculateToeplitzKernel_explicit (2D, Float64)
for Nx ∈ [32, 33, 64]
trj = rand(2, 1000) .- 0.5
Kx = NFFTTools.calculateToeplitzKernel_explicit((Nx, Nx), trj)
Ka = calculateToeplitzKernel((Nx, Nx), trj, m=4, σ=2)
@test typeof(Kx) === Matrix{ComplexF64}
@test typeof(Ka) === Matrix{ComplexF64}
@test Ka ≈ Kx rtol = 1e-6
end
## calculateToeplitzKernel with less-allocating constructor vs calculateToeplitzKernel_explicit (2D, Float64)
for Nx ∈ [32, 33, 64]
trj = rand(2, 1000) .- 0.5
Kx = NFFTTools.calculateToeplitzKernel_explicit((Nx, Nx), trj)
p = NFFTPlan(trj, (2Nx, 2Nx))
Ka = similar(Kx)
fftplan = plan_fft(Ka; flags = FFTW.ESTIMATE)
calculateToeplitzKernel!(Ka, p, trj, fftplan)
@test typeof(Kx) === Matrix{ComplexF64}
@test typeof(Ka) === Matrix{ComplexF64}
@test Ka ≈ Kx rtol = 1e-6
end
## calculateToeplitzKernel vs calculateToeplitzKernel_explicit (2D, Float32)
Nx = 32
trj = Float32.(rand(2, 10000) .- 0.5)
Kx = NFFTTools.calculateToeplitzKernel_explicit((Nx, Nx), trj)
Ka = calculateToeplitzKernel((Nx, Nx), trj, m=4, σ=2)
@test typeof(Kx) === Matrix{ComplexF32}
@test typeof(Ka) === Matrix{ComplexF32}
@test Ka ≈ Kx rtol = 1e-5
## compare to the NFFT
x = randn(ComplexF32, Nx, Nx)
xN = nfft_adjoint(trj, (Nx, Nx), nfft(trj, x))
xOS1 = similar(Ka)
xOS2 = similar(Ka)
FFTW.set_num_threads(1)
fftplan = plan_fft(xOS1; flags = FFTW.MEASURE)
ifftplan = plan_ifft(xOS1; flags = FFTW.MEASURE)
convolveToeplitzKernel!(x, Ka, fftplan, ifftplan, xOS1, xOS2)
@test x ≈ xN rtol = 1e-5
## test that convolveToeplitzKernel! with all arguments is non-allocating (only true with FFTW.set_num_threads(1))
bm = @benchmark convolveToeplitzKernel!($x, $Ka, $fftplan, $ifftplan, $xOS1, $xOS2)
@test bm.allocs == 0
## calculateToeplitzKernel vs calculateToeplitzKernel_explicit (2D-rectangular, Float32)
Nx = 32
Ny = 33
trj = Float32.(rand(2, 1000) .- 0.5)
Kx = NFFTTools.calculateToeplitzKernel_explicit((Nx, Ny), trj)
Ka = calculateToeplitzKernel((Nx, Ny), trj, m=4, σ=2)
@test typeof(Kx) === Matrix{ComplexF32}
@test typeof(Ka) === Matrix{ComplexF32}
@test size(Ka) == (2Nx, 2Ny)
@test Ka ≈ Kx rtol = 1e-5
## calculateToeplitzKernel vs calculateToeplitzKernel_explicit 3D, Float32
Nx = 32
trj = Float32.(rand(3, 1000) .- 0.5)
Kx = NFFTTools.calculateToeplitzKernel_explicit((Nx, Nx, Nx), trj)
Ka = calculateToeplitzKernel((Nx, Nx, Nx), trj, m=4, σ=2)
@test typeof(Kx) === Array{ComplexF32,3}
@test typeof(Ka) === Array{ComplexF32,3}
@test Ka ≈ Kx rtol = 1e-5
end