Clean tests (#6)

* clean tests

* fix lit

* DOCUMENTER_KEY

* docstring

* adjoint rotate

* more test work

.github/workflows/Documentation.yml

@@ -10,17 +10,12 @@ on:
 
 jobs:
   build:
-    runs-on: ${{ matrix.os }}
-    strategy:
-      matrix:
-        julia-version: ["1"]
-        julia-arch: [x86]
-        os: [ubuntu-latest]
+    runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: ${{ matrix.julia-version }}
+          version: '1.6'
       - name: Cache artifacts
         uses: actions/cache@v1
         env:
@@ -35,10 +30,14 @@ jobs:
       - name: Install dependencies
         run: |
           julia --project=docs/ -e '
+          ENV["JULIA_PKG_SERVER"] = ""
           using Pkg
           Pkg.develop(PackageSpec(path=pwd()))
           Pkg.instantiate()'
       - name: BuildAndDeploy
         env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+# https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#Authentication:-GITHUB_TOKEN
+#         GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          ssh: ${{ secrets.DOCUMENTER_KEY }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
         run: julia --project=docs/ docs/make.jl

src/spectplan.jl

@@ -17,14 +17,15 @@ Struct for storing key factors for a SPECT system model
 - `dy` voxel size in y direction (dx is the same value)
 - `imgsize{nx, ny, nz}` number of voxels in {x,y,z} direction of the image, must be integer
 - `psfsize{nx_psf, nz_psf}` number of voxels in {x, z} direction of the psf, must be integer
-- `pad_fft{padu_fft, pad_fft, padl_fft, padr_fft} pixels padded for {left,right,up,down} direction for convolution with psfs, must be integer
+- `pad_fft{padu_fft, pad_fft, padl_fft, padr_fft}` pixels padded for {left,right,up,down} direction for convolution with psfs, must be integer
 - `pad_rot{padu_rot, padd_rot, padl_rot, padr_rot}` padded pixels for {left,right,up,down} direction for image rotation
 - `ncore` number of CPU cores used to process data, must be integer
-Currently code assumes each of the `nview` projection views is `[nx,nz]`
-Currently code assumes `nx = ny`
-Currently code assumes uniform angular sampling
-Currently code uses multiprocessing using # of cores specified by Threads.nthreads() in Julia
-Currently code assumes psf is symmetric
+Currently code assumes the following:
+* each of the `nview` projection views is `[nx,nz]`
+* `nx = ny`
+* uniform angular sampling
+* `psf` is symmetric
+* multiprocessing using # of cores specified by `Threads.nthreads()`
 """
 struct SPECTplan
     T::DataType # default type for work arrays etc.
@@ -51,20 +52,19 @@ struct SPECTplan
         interpidx::Int = 2, # 1 is for 1d interpolation, 2 is for 2d interpolation
         T::DataType = promote_type(eltype(mumap), Float32),
     )
-        # check nx = ny ? typically 128 x 128 x 81
-        (nx, ny, nz) = size(mumap)
-        @assert isequal(nx, ny)
-        @assert iseven(nx) && iseven(ny)
+        (nx, ny, nz) = size(mumap) # typically 128 x 128 x 81
+        isequal(nx, ny) || throw("nx != ny")
+        (iseven(nx) && iseven(ny)) || throw("nx odd")
 
         # check psf
         nx_psf = size(psfs, 1)
         nz_psf = size(psfs, 2)
         nview = size(psfs, 4)
-        @assert isodd(nx_psf) && isodd(nz_psf)
-        @assert all(mapslices(x -> x == reverse(x), psfs, dims = [1, 2]))
+        (isodd(nx_psf) && isodd(nz_psf)) || throw("non-odd size psfs")
+        all(mapslices(x -> x == reverse(x), psfs, dims = [1, 2])) || throw("asym. psf")
 
         # check interpidx
-        @assert interpidx == 1 || interpidx == 2
+        (interpidx == 1 || interpidx == 2) || throw("bad interpidx")
 
         padu_fft = _padup(mumap, psfs)
         padd_fft = _paddown(mumap, psfs)
@@ -101,6 +101,24 @@ struct SPECTplan
 end
 
 
+"""
+    SPECTplan(mumap, psfs, dy; viewangle, interpidx, T)
+
+Constructor for `SPECTplan`
+
+In
+* `mumap::AbstractArray{<:RealU, 3}` 3D attenuation map
+* `psfs::AbstractArray{<:RealU, 4}` 4D PSF array
+* `dy::RealU` pixel size
+
+Option
+* `viewangle::AbstractVector{<:RealU}` default 0 to almost 2π
+* `interpidx::Int = 2` 1 is for 3-pass 1D interpolator; 2 is for 2D interpolator
+* `T::DataType = promote_type(eltype(mumap), Float32)`
+"""
+function SPECTplan() end
+
+
 """
     Workarray
 Struct for storing keys of the work array for a single thread

test/Project.toml

@@ -1,4 +1,17 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearInterpolators = "2fc109c4-9d77-11e9-32ff-e5c650334c16"
+LinearMapsAA = "599c1a8e-b958-11e9-0d14-b1e6b2ecea07"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+BenchmarkTools = "1.2"
+FFTW = "1.4.5"
+ImageFiltering = "0.7"
+LinearInterpolators = "0.1.5"
+OffsetArrays = "1.10"

test/adjoint-fftconv.jl

@@ -0,0 +1,42 @@
+# adjoint-fftconv.jl
+# test adjoint consistency for FFT convolution methods on very small case
+
+using SPECTrecon: fft_conv!, fft_conv_adj!, Power2
+using FFTW: plan_fft!, plan_ifft!
+using LinearMapsAA: LinearMapAA
+using Test: @test, @testset
+
+
+@testset "adjoint-fftconv" begin
+    M = 200 # todo: make smaller
+#   M = 20
+    N = 64
+    T = Float32
+    fftpadsize = (28, 28, 32, 32)
+    img_compl = zeros(Complex{T}, 256, 128)
+    ker_compl = similar(img_compl)
+    workmat = zeros(T, 256, 128)
+    workvec1 = zeros(T, 128)
+    workvec2 = zeros(T, 256)
+    fft_plan = plan_fft!(img_compl)
+    ifft_plan = plan_ifft!(img_compl)
+
+#   x = randn(T, M, N)
+#   output_x = similar(x)
+#   y = randn(T, M, N)
+#   output_y = similar(y)
+    ker = rand(T, 11, 11) # change from 5 -> 11
+    ker /= sum(ker)
+    kerev = reverse(ker)
+
+    forw!(y,x) = fft_conv!(y, workmat, x,
+        ker, fftpadsize, img_compl, ker_compl, fft_plan, ifft_plan)
+
+    back!(x,y) = fft_conv_adj!(x, workmat, workvec1, workvec2, y,
+        kerev, fftpadsize, img_compl, ker_compl, fft_plan, ifft_plan)
+
+    idim = (M,N)
+    odim = (M,N)
+    A = LinearMapAA(forw!, back!, (prod(odim), prod(idim)); T, odim, idim)
+#   @test Matrix(A') ≈ Matrix(A)' # todo
+end

test/adjoint-project.jl

@@ -0,0 +1,50 @@
+# adjoint-project.jl
+# test adjoint consistency for SPECT projector/back-projector on very small case
+
+using SPECTrecon: project, backproject
+using SPECTrecon: project!, backproject!
+using SPECTrecon: SPECTplan, Workarray
+using LinearMapsAA: LinearMapAA
+using Test: @test, @testset
+
+
+#@testset "adjoint" begin
+    T = Float32
+    nx = 16
+    ny = 16
+    nz = 5
+    nview = 7
+
+    mumap = rand(T, nx, ny, nz)
+
+    nx_psf = 5
+    nz_psf = 3
+    psfs = rand(T, nx_psf, nz_psf, ny, nview)
+    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2]) # symmetrize
+    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])
+
+    dy = 4.7952
+
+#   for interpidx = 1:2
+    interpidx = 1
+        plan = SPECTplan(mumap, psfs, dy; interpidx)
+#       workarray = Workarray(plan.T, plan.imgsize, plan.pad_fft, plan.pad_rot)
+        workarray = [Workarray(plan.T, plan.imgsize, plan.pad_fft, plan.pad_rot) for i in 1:plan.ncore] # allocate
+#       workarray = [workarray]
+        forw = x -> project(x, mumap, psfs, dy; interpidx)
+        back = y -> backproject(y, mumap, psfs, dy; interpidx)
+        forw! = (y,x) -> project!(y, x, plan, workarray)
+        back! = (x,y) -> backproject!(x, y, plan, workarray)
+        idim = (nx,ny,nz)
+        odim = (nx,nz,nview)
+        A0 = LinearMapAA(forw, back, (prod(odim),prod(idim)); T, idim, odim)
+        A! = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, idim, odim)
+        x = 0 * rand(T, idim)
+        y = Array{T}(undef, odim)
+        @show extrema(forw(x))
+        @show extrema(forw!(y,x)) # NaN values even for 0 input :(
+#       @test forw!(y,x) ≈ forw(x)
+#       @test Matrix(A0)' ≈ Matrix(A0')
+#       @test Matrix(A!)' ≈ Matrix(A!')
+#    end
+#end

test/adjoint-rotate.jl

@@ -0,0 +1,45 @@
+# adjoint-rotate.jl
+# test adjoint consistency for rotate methods on very small case
+
+using SPECTrecon: linearinterp!, rotate_x!, rotate_y!
+using SPECTrecon: rotate_x_adj!, rotate_y_adj!
+using SPECTrecon: rotl90!, rotr90!, rot180!
+using SPECTrecon: imrotate3!, imrotate3_adj!
+using LinearAlgebra: dot
+using LinearInterpolators: SparseInterpolator, LinearSpline
+using LinearMapsAA: LinearMapAA
+using Test: @test, @testset
+using Random: seed!
+
+
+@testset "adjoint-rotate" begin
+    Ntest = 7
+    θ_list = (0:Ntest-1) / Ntest * 2π
+    M,N = 16, 16 # todo: test non-square?
+    T = Float32
+    pad_x = ceil(Int, 1 + M * sqrt(2)/2 - M / 2)
+    pad_y = ceil(Int, 1 + N * sqrt(2)/2 - N / 2)
+    workvec_x = zeros(T, M + 2 * pad_x)
+    workvec_y = zeros(T, N + 2 * pad_y)
+    workmat2_x = Array{T}(undef, M + 2 * pad_x, N + 2 * pad_y)
+    workmat2_y = Array{T}(undef, M + 2 * pad_x, N + 2 * pad_y)
+    workmat1_x = Array{T}(undef, M + 2 * pad_x, N + 2 * pad_y)
+    workmat1_y = Array{T}(undef, M + 2 * pad_x, N + 2 * pad_y)
+    A_x = SparseInterpolator(LinearSpline(T), workvec_x, length(workvec_x))
+    A_y = SparseInterpolator(LinearSpline(T), workvec_y, length(workvec_y))
+
+    for θ in θ_list
+		idim = (M,N)
+		odim = (M,N)
+
+        forw1!(y, x) = imrotate3!(y, workmat1_x, workmat2_x, x, θ)
+        back1!(x, y) = imrotate3_adj!(x, workmat1_y, workmat2_y, y, θ)
+		A = LinearMapAA(forw1!, back1!, (prod(odim), prod(idim)); T, odim, idim)
+        @test Matrix(A') ≈ Matrix(A)' # 3-pass 1D version
+
+        forw2!(y, x) = imrotate3!(y, workmat1_x, workmat2_x, x, θ, A_x, A_y, workvec_x, workvec_y)
+        back2!(x, y) = imrotate3_adj!(x, workmat1_y, workmat2_y, y, θ, A_x, A_y, workvec_x, workvec_y)
+		A = LinearMapAA(forw2!, back2!, (prod(odim), prod(idim)); T, odim, idim)
+        @test Matrix(A') ≈ Matrix(A)' # 2D version
+    end
+end

test/fft_convolve.jl

@@ -1,11 +1,12 @@
 # fft_convolve.jl
 
-using Main.SPECTrecon: imfilter3!
-using Main.SPECTrecon: fft_conv!, fft_conv_adj!
-using LazyAlgebra:vdot
+using SPECTrecon: imfilter3!
+using SPECTrecon: fft_conv!, fft_conv_adj!
+using LinearAlgebra: dot
 using FFTW: plan_fft!, plan_ifft!
+using Random: seed!
 using ImageFiltering: centered, imfilter
-using Test: @test, @testset, detect_ambiguities
+using Test: @test, @testset
 
 
 @testset "imfilter3!" begin
@@ -39,6 +40,7 @@ end
     workvec2 = zeros(T, 256)
     fft_plan = plan_fft!(img_compl)
     ifft_plan = plan_ifft!(img_compl)
+    seed!(0) # todo: fails sometimes
     for t = 1:testnum
         x = randn(T, M, N)
         output_x = similar(x)
@@ -48,10 +50,10 @@ end
         ker /= sum(ker)
         kerev = reverse(ker)
         fft_conv!(output_x, workmat, x, ker, fftpadsize,
-                                img_compl, ker_compl, fft_plan, ifft_plan)
+            img_compl, ker_compl, fft_plan, ifft_plan)
 
         fft_conv_adj!(output_y, workmat, workvec1, workvec2, y, kerev, fftpadsize,
-                                    img_compl, ker_compl, fft_plan, ifft_plan)
-        @test isapprox(vdot(y, output_x), vdot(x, output_y))
+            img_compl, ker_compl, fft_plan, ifft_plan)
+        @test isapprox(dot(y, output_x), dot(x, output_y)) # todo tol?
     end
 end

test/helper.jl

@@ -1,13 +1,13 @@
 # helper.jl
 
-using Main.SPECTrecon: padzero!, padrepl!, pad2sizezero!, pad_it!
-using Main.SPECTrecon: fftshift!, ifftshift!, fftshift2!
-using Main.SPECTrecon: plus1di!, plus1dj!, plus2di!, plus2dj!
-using Main.SPECTrecon: plus3di!, plus3dj!, plus3dk!, scale3dj!, mul3dj!
-using Main.SPECTrecon: copy3dj!
+using SPECTrecon: padzero!, padrepl!, pad2sizezero!, pad_it!
+using SPECTrecon: fftshift!, ifftshift!, fftshift2!
+using SPECTrecon: plus1di!, plus1dj!, plus2di!, plus2dj!
+using SPECTrecon: plus3di!, plus3dj!, plus3dk!, scale3dj!, mul3dj!
+using SPECTrecon: copy3dj!
 using ImageFiltering: Fill, Pad, BorderArray
-using OffsetArrays
-using Test: @test, @testset, detect_ambiguities
+import OffsetArrays
+using Test: @test, @testset
 
 
 @testset "padzero!" begin

test/project-mem.jl

@@ -16,7 +16,7 @@ xtrue = rand(T, nx, nx, nz)
 
 dy = T(4.7952)
 nview = size(psfs, 4)
-plan = SPECTplan(mumap, psfs, nview, dy; interpidx = 1)
+plan = SPECTplan(mumap, psfs, dy; interpidx = 1)
 workarray = Vector{Workarray}(undef, plan.ncore)
 for i = 1:plan.ncore
     workarray[i] = Workarray(plan.T, plan.imgsize, plan.pad_fft, plan.pad_rot) # allocate

test/project.jl

@@ -1,8 +1,8 @@
 # project.jl
 
-using Main.SPECTrecon: project, backproject
-using LazyAlgebra: vdot
-using Test: @test, @testset, detect_ambiguities
+using SPECTrecon: project, backproject
+using LinearAlgebra: dot
+using Test: @test, @testset
 
 
 @testset "proj-adj-test" begin
@@ -26,9 +26,9 @@ using Test: @test, @testset, detect_ambiguities
 
     output_x = project(x, mumap, psfs, dy; interpidx = 1)
     output_y = backproject(y, mumap, psfs, dy; interpidx = 1)
-    @test isapprox(vdot(y, output_x), vdot(x, output_y))
+    @test isapprox(dot(y, output_x), dot(x, output_y))
 
     output_x = project(x, mumap, psfs, dy; interpidx = 2)
     output_y = backproject(y, mumap, psfs, dy; interpidx = 2)
-    @test isapprox(vdot(y, output_x), vdot(x, output_y))
+    @test isapprox(dot(y, output_x), dot(x, output_y))
 end