From 2d22ac86a53ce993ca97b42b69323fc990374545 Mon Sep 17 00:00:00 2001
From: Zongyu Li <54018575+ZongyuLi-umich@users.noreply.github.com>
Date: Thu, 19 May 2022 13:09:06 -0400
Subject: [PATCH 1/5] src code

---
 src/CuSPECTplan.jl    | 160 +++++++++++++++++++++
 src/CuSPECTrecon.jl   |  14 ++
 src/Cubackproject.jl  | 106 ++++++++++++++
 src/Cufft_convolve.jl | 145 +++++++++++++++++++
 src/Cuhelper.jl       | 316 ++++++++++++++++++++++++++++++++++++++++++
 src/Cuinterp.jl       |  40 ++++++
 src/Cuplan-fft.jl     | 122 ++++++++++++++++
 src/Cuplan-rotate.jl  |  89 ++++++++++++
 src/Cuproject.jl      | 108 +++++++++++++++
 src/Cupsf-gauss.jl    |  55 ++++++++
 src/Curotate.jl       |  52 +++++++
 11 files changed, 1207 insertions(+)
 create mode 100644 src/CuSPECTplan.jl
 create mode 100644 src/Cubackproject.jl
 create mode 100644 src/Cufft_convolve.jl
 create mode 100644 src/Cuhelper.jl
 create mode 100644 src/Cuinterp.jl
 create mode 100644 src/Cuplan-fft.jl
 create mode 100644 src/Cuplan-rotate.jl
 create mode 100644 src/Cuproject.jl
 create mode 100644 src/Cupsf-gauss.jl
 create mode 100644 src/Curotate.jl

diff --git a/src/CuSPECTplan.jl b/src/CuSPECTplan.jl
new file mode 100644
index 0000000..3bee563
--- /dev/null
+++ b/src/CuSPECTplan.jl
@@ -0,0 +1,160 @@
+# CuSPECTplan.jl
+
+export CuSPECTplan
+
+"""
+    CuSPECTplan
+Struct for storing key factors for a SPECT system model
+- `T` datatype of work arrays
+- `imgsize` size of image
+- `px,pz` psf dimension
+- `imgr [nx, ny, nz]` 3D rotated version of image
+- `add_img [nx, ny, nz]` 3D image for adding views and backprojection
+- `mumap [nx,ny,nz]` attenuation map, must be 3D, possibly zeros()
+- `mumapr [nx, ny, nz]` 3D rotated mumap
+- `exp_mumapr [nx, nz]` 2D exponential rotated mumap
+- `psfs [px,pz,ny,nview]` point spread function, must be 4D, with `px and `pz` odd, and symmetric for each slice
+- `nview` number of views, must be integer
+- `viewangle` set of view angles, must be from 0 to 2π
+- `dy` voxel size in y direction (dx is the same value)
+- `planrot` Vector of struct `CuPlanRotate`
+- `planpsf` Vector of struct `CuPlanFFT`
+Currently code assumes the following:
+* each of the `nview` projection views is `[nx,nz]`
+* `nx = ny`
+* uniform angular sampling
+* `psf` is symmetric
+"""
+struct CuSPECTplan{T,A2,A3,A4}
+    T::DataType # default type for work arrays etc.
+    imgsize::NTuple{3, Int}
+    px::Int
+    pz::Int
+    imgr::A3 # 3D rotated image, (nx, ny, nz)
+    add_img::A3
+    add_view::A2
+    mumap::A3 # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
+    mumapr::A3 # 3D rotated mumap, (nx, ny, nz)
+    exp_mumapr::A2 # 2D exponential rotated mumap, (nx, ny)
+    psfs::A4 # PSFs must be 4D, [px, pz, ny, nview], finally be centered psf
+    nview::Int # number of views
+    viewangle::StepRangeLen{T}
+    dy::T
+    planrot::CuPlanRotate # todo: make it concrete?
+    planpsf::CuPlanFFT
+
+    """
+        CuSPECTplan(mumap, psfs, dy; T, viewangle)
+    """
+    function CuSPECTplan(
+        mumap::AbstractArray{<:RealU,3},
+        psfs::AbstractArray{<:RealU,4},
+        dy::RealU;
+        T::DataType = promote_type(eltype(mumap), eltype(psfs), Float32),
+        viewangle::StepRangeLen{<:RealU} = (0:size(psfs, 4) - 1) / size(psfs, 4) * T(2π), # set of view angles
+        )
+
+        # convert to the same type
+        dy = convert(T, dy)
+        mumap .= T.(mumap)
+        psfs .= T.(psfs)
+
+        # convert to gpu if not
+        mumap = typeof(mumap) <: CuArray ? mumap : CuArray(mumap)
+        psfs = typeof(psfs) <: CuArray ? psfs : CuArray(psfs)
+
+        (nx, ny, nz) = size(mumap) # typically 128 x 128 x 81
+
+        isequal(nx, ny) || throw("nx != ny")
+        (iseven(nx) && iseven(ny)) || throw("nx odd")
+
+        imgsize = (nx, ny, nz)
+        # check psf
+        px, pz, _, nview = size(psfs)
+        (isodd(px) && isodd(pz)) || throw("non-odd size psfs")
+        psfs == reverse(reverse(psfs, dims = 1), dims = 2) ||
+            throw("asym. psf")
+        # all(mapslices(x -> x == reverse(x, dims=:), psfs, dims = [1, 2])) ||
+        #     throw("asym. psf")
+        # imgr stores 3D image in different view angles
+        imgr = CuArray{T, 3}(undef, nx, ny, nz)
+        # add_img stores 3d image for backprojection
+        add_img = CuArray{T, 3}(undef, nx, ny, nz)
+        # mumapr stores 3D mumap in different view angles
+        mumapr = CuArray{T, 3}(undef, nx, ny, nz)
+
+        A2 = typeof(mumap[:,:,1])
+        A3 = typeof(mumap)
+        A4 = typeof(psfs)
+
+        exp_mumapr = CuArray{T, 2}(undef, nx, nz)
+        add_view = CuArray{T, 2}(undef, nx, nz)
+
+        planrot = CuPlanRotate(nx; T)
+
+        planpsf = CuPlanFFT(; nx, nz, px, pz, T)
+
+        new{T, A2, A3, A4}(
+            T, # default type for work arrays etc.
+            imgsize,
+            px,
+            pz,
+            imgr, # 3D rotated image, (nx, ny, nz)
+            add_img,
+            add_view,
+            mumap, # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
+            mumapr, # 3D rotated mumap, (nx, ny, nz)
+            exp_mumapr,
+            psfs, # PSFs must be 4D, [px, pz, ny, nview], finally be centered psf
+            nview, # number of views
+            viewangle,
+            dy,
+            planrot,
+            planpsf,
+        )
+    end
+end
+
+
+"""
+    Cugen_attenuation!(plan, y)
+Generate depth-dependent attenuation map
+"""
+function Cugen_attenuation!(plan::CuSPECTplan, y::Int)
+    Cuscale3dj!(plan.exp_mumapr, plan.mumapr, y, -0.5)
+    for j = 1:y
+        Cuplus3dj!(plan.exp_mumapr, plan.mumapr, j)
+    end
+
+    broadcast!(*, plan.exp_mumapr, plan.exp_mumapr, - plan.dy)
+    broadcast!(exp, plan.exp_mumapr, plan.exp_mumapr)
+    return plan.exp_mumapr
+end
+
+
+"""
+    show(io::IO, ::MIME"text/plain", plan::CuSPECTplan)
+"""
+function Base.show(io::IO, ::MIME"text/plain", plan::CuSPECTplan{T,A2,A3,A4}) where {T,A2,A3,A4}
+    t = typeof(plan)
+    println(io, t)
+    for f in (:imgsize, :px, :pz, :nview, :viewangle, :dy)
+        p = getproperty(plan, f)
+        t = typeof(p)
+        println(io, " ", f, "::", t, " ", p)
+    end
+    for f in (:mumap, )
+        p = getproperty(plan, f)
+        println(io, " ", f, ":", " ", summary(p))
+    end
+    println(io, " (", sizeof(plan), " bytes)")
+end
+
+
+"""
+    sizeof(::CuSPECTplan)
+Show size in bytes of `CuSPECTplan` object.
+"""
+function Base.sizeof(ob::T) where {T <: Union{CuSPECTplan}}
+    sum(f -> sizeof(getfield(ob, f)), fieldnames(typeof(ob)))
+end
diff --git a/src/CuSPECTrecon.jl b/src/CuSPECTrecon.jl
index 875c2e7..84709ba 100644
--- a/src/CuSPECTrecon.jl
+++ b/src/CuSPECTrecon.jl
@@ -1,3 +1,17 @@
+"""
+    CuSPECTrecon
+GPU version of system matrix (forward and back-projector) for SPECT image reconstruction.
+"""
 module CuSPECTrecon
 
+    include("Cuhelper.jl")
+    include("Cuplan-rotate.jl")
+    include("Curotate.jl")
+    include("Cuplan-fft.jl")
+    include("Cupsf-gauss.jl")
+    include("Cufft_convolve.jl")
+    include("CuSPECTplan.jl")
+    include("Cuproject.jl")
+    include("Cubackproject.jl")
+
 end # module
diff --git a/src/Cubackproject.jl b/src/Cubackproject.jl
new file mode 100644
index 0000000..3e32bcb
--- /dev/null
+++ b/src/Cubackproject.jl
@@ -0,0 +1,106 @@
+# Cubackproject.jl
+
+export Cubackproject, Cubackproject!
+
+
+"""
+    Cubackproject!(image, view, plan, viewidx)
+Backproject a single view.
+"""
+function Cubackproject!(
+    image::AbstractArray{<:RealU, 3},
+    view::AbstractMatrix{<:RealU},
+    plan::CuSPECTplan,
+    viewidx::Int,
+)
+
+    for z = 1:plan.imgsize[3] # 1:nz
+        # rotate mumap by "-angle"
+        Cuimrotate!((@view plan.mumapr[:, :, z]),
+                    (@view plan.mumap[:, :, z]),
+                    -plan.viewangle[viewidx],
+                    plan.planrot)
+    end
+
+    # adjoint of convolving img with psf and applying attenuation map
+    for y = 1:plan.imgsize[2] # 1:ny
+
+        Cugen_attenuation!(plan, y)
+
+        Cufft_conv_adj!((@view plan.imgr[:, y, :]),
+                         view,
+                        (@view plan.psfs[:, :, y, viewidx]),
+                         plan.planpsf)
+
+        Cumul3dj!(plan.imgr, plan.exp_mumapr, y)
+    end
+
+    # adjoint of rotating image, again, rotating by "-angle"
+    for z = 1:plan.imgsize[3] # 1:nz
+
+        Cuimrotate!((@view image[:, :, z]),
+                    (@view plan.imgr[:, :, z]),
+                    -plan.viewangle[viewidx],
+                     plan.planrot)
+    end
+
+    return image
+end
+
+
+"""
+    Cubackproject!(image, views, plan ; viewlist)
+Backproject multiple views into `image`.
+Users must initialize `image` to zero.
+"""
+function Cubackproject!(
+    image::AbstractArray{<:RealU, 3},
+    views::AbstractArray{<:RealU, 3},
+    plan::CuSPECTplan;
+    viewlist::AbstractVector{<:Int} = 1:plan.nview, # all views
+)
+
+    # loop over each view index
+    for viewidx in viewlist
+        Cubackproject!(plan.add_img, (@view views[:, :, viewidx]), plan, viewidx)
+        broadcast!(+, image, image, plan.add_img)
+    end
+    return image
+end
+
+
+"""
+    image = Cubackproject(views, plan ; kwargs...)
+SPECT backproject `views`; this allocates the returned 3D array.
+"""
+function Cubackproject(
+    views::AbstractArray{<:RealU, 3},
+    plan::CuSPECTplan;
+    kwargs...,
+)
+    image = CuArray(zeros(plan.T, plan.imgsize))
+    Cubackproject!(image, views, plan; kwargs...)
+    return image
+end
+
+
+"""
+    image = Cubackproject(views, mumap, psfs, dy; kwargs...)
+SPECT backproject `views` using attenuation map `mumap` and PSF array `psfs` for pixel size `dy`.
+This method initializes the `plan` as a convenience.
+Most users should use `Cubackproject!` instead after initializing those, for better efficiency.
+"""
+function Cubackproject(
+    views::AbstractArray{<:RealU, 3}, # [nx,nz,nview]
+    mumap::AbstractArray{<:RealU, 3}, # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
+    psfs::AbstractArray{<:RealU, 4},
+    dy::RealU;
+    kwargs...,
+)
+
+    size(mumap,1) == size(mumap,1) == size(views,1) ||
+        throw(DimensionMismatch("nx"))
+    size(mumap,3) == size(views,2) || throw(DimensionMismatch("nz"))
+    plan = CuSPECTplan(mumap, psfs, dy; kwargs...)
+    return Cubackproject(views, plan; kwargs...)
+end
diff --git a/src/Cufft_convolve.jl b/src/Cufft_convolve.jl
new file mode 100644
index 0000000..4ff00f5
--- /dev/null
+++ b/src/Cufft_convolve.jl
@@ -0,0 +1,145 @@
+# Cufft_convolve.jl
+
+export Cufft_conv!, Cufft_conv_adj!
+export Cufft_conv, Cufft_conv_adj
+
+"""
+    Cuimfilterz!(plan)
+GPU version of FFT-based convolution of `plan.img_compl`
+and kernel `plan.ker_compl` (not centered),
+storing result in `plan.workmat`.
+"""
+function Cuimfilterz!(plan::CuPlanFFT)
+    mul!(plan.img_compl, plan.fft_plan, plan.img_compl)
+    mul!(plan.ker_compl, plan.fft_plan, plan.ker_compl)
+    broadcast!(*, plan.img_compl, plan.img_compl, plan.ker_compl)
+    mul!(plan.img_compl, plan.ifft_plan, plan.img_compl)
+    fftshift!(plan.ker_compl, plan.img_compl)
+    plan.workmat .= real.(plan.ker_compl)
+    return plan.workmat
+end
+
+
+"""
+    Cufft_conv!(output, img, ker, plan)
+Convolve 2D image `img` with 2D (symmetric!) kernel `ker` using FFT,
+storing the result in `output`.
+"""
+function Cufft_conv!(
+    output::AbstractMatrix{<:RealU},
+    img::AbstractMatrix{<:RealU},
+    ker::AbstractMatrix{<:RealU},
+    plan::CuPlanFFT,
+)
+    @boundscheck size(output) == size(img) || throw("size output")
+    @boundscheck size(img) == (plan.nx, plan.nz) || throw("size img")
+    @boundscheck size(ker) == (plan.px, plan.pz) ||
+        throw("size ker $(size(ker)) $(plan.px) $(plan.pz)")
+
+    # filter image with a kernel, using replicate padding and fft convolution
+    Cupadrepl!(plan.img_compl, img, plan.padsize)
+
+    Cupad2sizezero!(plan.ker_compl, ker, size(plan.ker_compl)) # zero pad kernel
+
+    Cuimfilterz!(plan)
+
+    (M, N) = size(img)
+    copyto!(output, (@view plan.workmat[plan.padsize[1] .+ (1:M),
+                                        plan.padsize[3] .+ (1:N)]))
+    return output
+end
+
+
+"""
+    Cufft_conv(img, ker; T)
+GPU version of convolving 2D image `img` with 2D (symmetric!) kernel `ker` using FFT.
+"""
+function Cufft_conv(
+    img::AbstractMatrix{I},
+    ker::AbstractMatrix{K};
+    T::DataType = promote_type(I, K, Float32),
+) where {I <: Number, K <: Number}
+
+    ker ≈ reverse(ker, dims=:) || throw("asymmetric kernel")
+    nx, nz = size(img)
+    px, pz = size(ker)
+    plan = CuPlanFFT( ; nx, nz, px, pz, T)
+    output = similar(img)
+    Cufft_conv!(output, img, ker, plan)
+    return output
+end
+
+
+"""
+    Cufft_conv_adj!(output, img, ker, plan)
+GPU version of the adjoint of `fft_conv!`.
+"""
+function Cufft_conv_adj!(
+    output::AbstractMatrix{<:RealU},
+    img::AbstractMatrix{<:RealU},
+    ker::AbstractMatrix{<:RealU},
+    plan::CuPlanFFT{T},
+) where {T}
+
+    @boundscheck size(output) == size(img) || throw("size output")
+    @boundscheck size(img) == (plan.nx, plan.nz) || throw("size img")
+    @boundscheck size(ker) == (plan.px, plan.pz) ||
+        throw("size ker $(size(ker)) $(plan.px) $(plan.pz)")
+
+    Cupadzero!(plan.img_compl, img, plan.padsize) # pad the image with zeros
+    Cupad2sizezero!(plan.ker_compl, ker, size(plan.ker_compl)) # pad the kernel with zeros
+
+    Cuimfilterz!(plan)
+    (M, N) = size(img)
+    # adjoint of replicate padding
+    plan.workvecz .= zero(T)
+    for i = 1:plan.padsize[1]
+        Cuplus2di!(plan.workvecz, plan.workmat, i)
+    end
+    Cuplus1di!(plan.workmat, plan.workvecz, 1+plan.padsize[1])
+
+    plan.workvecz .= zero(T)
+    for i = (plan.padsize[1]+M+1):size(plan.workmat, 1)
+        Cuplus2di!(plan.workvecz, plan.workmat, i)
+    end
+    Cuplus1di!(plan.workmat, plan.workvecz, M+plan.padsize[1])
+
+    plan.workvecx .= zero(T)
+    for j = 1:plan.padsize[3]
+        Cuplus2dj!(plan.workvecx, plan.workmat, j)
+    end
+    Cuplus1dj!(plan.workmat, plan.workvecx, 1+plan.padsize[3])
+
+    plan.workvecx .= zero(T)
+    for j = (plan.padsize[3]+N+1):size(plan.workmat, 2)
+        Cuplus2dj!(plan.workvecx, plan.workmat, j)
+    end
+    Cuplus1dj!(plan.workmat, plan.workvecx, N+plan.padsize[3])
+
+    copyto!(output,
+        (@view plan.workmat[(plan.padsize[1]+1):(plan.padsize[1]+M),
+                            (plan.padsize[3]+1):(plan.padsize[3]+N)]),
+    )
+
+    return output
+end
+
+
+"""
+    Cufft_conv_adj(img, ker; T)
+GPU version of the adjoint of `fft_conv`.
+"""
+function Cufft_conv_adj(
+    img::AbstractMatrix{I},
+    ker::AbstractMatrix{K};
+    T::DataType = promote_type(I, K, Float32),
+) where {I <: Number, K <: Number}
+
+    ker ≈ reverse(ker, dims=:) || throw("asymmetric kernel")
+    nx, nz = size(img)
+    px, pz = size(ker)
+    plan = CuPlanFFT( ; nx, nz, px, pz, T)
+    output = similar(Matrix{T}, size(img))
+    Cufft_conv_adj!(output, img, ker, plan)
+    return output
+end
diff --git a/src/Cuhelper.jl b/src/Cuhelper.jl
new file mode 100644
index 0000000..f5f2f12
--- /dev/null
+++ b/src/Cuhelper.jl
@@ -0,0 +1,316 @@
+# Cuhelper.jl
+# A lot of helper functions
+
+export Cupadzero!, Cupadrepl!, Cupad2sizezero!
+export fftshift!, ifftshift!
+export Cuplus1di!, Cuplus1dj!
+export Cuplus2di!, Cuplus2dj!
+export Cuplus3di!, Cuplus3dj!, Cuplus3dk!
+export Cuscale3dj!, Cumul3dj!, Cucopy3dj!
+export Curotl90!, Curotr90!, Curot180!, Curot_f90!
+using LinearAlgebra
+using FFTW
+
+const RealU = Number
+Power2 = x -> 2^(ceil(Int, log2(x)))
+_padup(nx, px)    =  ceil(Int, (Power2(nx + px - 1) - nx) / 2)
+_paddown(nx, px)  = floor(Int, (Power2(nx + px - 1) - nx) / 2)
+_padleft(nz, pz)  =  ceil(Int, (Power2(nz + pz - 1) - nz) / 2)
+_padright(nz, pz) = floor(Int, (Power2(nz + pz - 1) - nz) / 2)
+
+
+"""
+    Cupadzero!(output, img, padsize)
+GPU version of padzero!
+"""
+function Cupadzero!(
+    output::AbstractMatrix{T},
+    img::AbstractMatrix,
+    padsize::NTuple{4, <:Int}, # up, down, left, right
+    ) where {T}
+
+    @boundscheck size(output) ==
+        size(img) .+ (padsize[1] + padsize[2], padsize[3] + padsize[4]) || throw("size")
+    M, N = size(img)
+    output .= zero(T)
+    (@view output[padsize[1] + 1:padsize[1] + M, padsize[3] + 1:padsize[3] + N]) .= img
+    return output
+
+end
+
+
+"""
+    Cupadrepl!(output, img, padsize)
+GPU version of padrepl!
+"""
+function Cupadrepl!(output::AbstractMatrix{T},
+                    img::AbstractMatrix,
+                    padsize::NTuple{4, <:Int}
+                    ) where {T}
+
+    @boundscheck size(output) ==
+        size(img) .+ (padsize[1] + padsize[2], padsize[3] + padsize[4]) || throw("size")
+    M, N = size(img)
+    output .= zero(T)
+    (@view output[padsize[1] + 1:padsize[1] + M, padsize[3] + 1:padsize[3] + N]) .= img
+    (@view output[1:padsize[1], padsize[3] + 1:padsize[3] + N]) .= (@view img[1:1, :])
+    (@view output[padsize[1] + M + 1 : end, padsize[3] + 1 : padsize[3] + N]) .= (@view img[end:end, :])
+    (@view output[:, 1:padsize[3]]) .= (@view output[:, padsize[3] + 1 : padsize[3] + 1])
+    (@view output[:, padsize[3] + N + 1: end]) .= (@view output[:, padsize[3] + N : padsize[3] + N])
+    return output
+end
+
+
+"""
+    Cupad2sizezero!(output, img, padsize)
+GPU version of pad2sizezero!
+"""
+function Cupad2sizezero!(output::AbstractMatrix{T},
+                         img::AbstractMatrix,
+                         padsize::Tuple{<:Int, <:Int},
+                         ) where {T}
+
+        @boundscheck size(output) == padsize || throw("size")
+        dims = size(img)
+        pad_dims = ceil.(Int, (padsize .- dims) ./ 2)
+        output .= zero(T)
+        (@view output[pad_dims[1] + 1:pad_dims[1] + dims[1], pad_dims[2] + 1:pad_dims[2] + dims[2]]) .= img
+        return output
+
+end
+
+
+"""
+    fftshift!(dst, src), ifftshift!(dst, src)
+fftshift! and ifftshift! work for both cpu and gpu
+"""
+fftshift!(dst::AbstractArray, src::AbstractArray) = circshift!(dst, src, size(src) .÷ 2)
+
+ifftshift!(dst::AbstractArray, src::AbstractArray) = circshift!(dst, src, size(src) .÷ -2)
+
+
+"""
+    Cuplus1di!(mat2d, mat1d, i)
+GPU version of `mat2d[i, :] += mat1d`
+"""
+Base.@propagate_inbounds function Cuplus1di!(
+    mat2d::AbstractMatrix,
+    mat1d::AbstractVector,
+    i::Int,
+)
+    @boundscheck (size(mat1d, 1) == size(mat2d, 2) || throw("size2"))
+    @boundscheck (1 ≤ i ≤ size(mat2d, 1) || throw("bad i"))
+    (@view mat2d[i, :]) .+= mat1d
+    return mat2d
+end
+
+
+"""
+    Cuplus1dj!(mat2d, mat1d, j)
+GPU version of `mat2d[:, j] += mat1d`
+"""
+Base.@propagate_inbounds function Cuplus1dj!(
+    mat2d::AbstractMatrix,
+    mat1d::AbstractVector,
+    j::Int,
+)
+    @boundscheck (size(mat1d, 1) == size(mat2d, 1) || throw("size1"))
+    @boundscheck (1 ≤ j ≤ size(mat2d, 2) || throw("bad j"))
+    (@view mat2d[:, j]) .+= mat1d
+    return mat2d
+end
+
+
+"""
+    Cuplus2di!(mat1d, mat2d, i)
+GPU version of `mat1d += mat2d[i,:]`
+"""
+Base.@propagate_inbounds function Cuplus2di!(
+    mat1d::AbstractVector,
+    mat2d::AbstractMatrix,
+    i::Int,
+)
+    @boundscheck (size(mat1d, 1) == size(mat2d, 2) || throw("size2"))
+    @boundscheck (1 ≤ i ≤ size(mat2d, 1) || throw("bad i"))
+    mat1d .+= (@view mat2d[i, :])
+    return mat1d
+end
+
+
+"""
+    Cuplus2dj!(mat1d, mat2d, j)
+GPU version of plus2dj!
+"""
+Base.@propagate_inbounds function Cuplus2dj!(
+    mat1d::AbstractVector,
+    mat2d::AbstractMatrix,
+    j::Int,
+)
+    @boundscheck (size(mat1d, 1) == size(mat2d, 1) || throw("size1"))
+    @boundscheck (1 ≤ j ≤ size(mat2d, 2) || throw("bad j"))
+    mat1d .+= (@view mat2d[:, j])
+    return mat1d
+end
+
+
+"""
+    Cuplus3di!(mat2d, mat3d, i)
+GPU version of plus3di!
+"""
+Base.@propagate_inbounds function Cuplus3di!(
+    mat2d::AbstractMatrix,
+    mat3d::AbstractArray{<:Any, 3},
+    i::Int,
+)
+    @boundscheck (size(mat2d, 1) == size(mat3d, 2) || throw("size2"))
+    @boundscheck (size(mat2d, 2) == size(mat3d, 3) || throw("size3"))
+    @boundscheck (1 ≤ i ≤ size(mat3d, 1) || throw("bad i"))
+    mat2d .+= (@view mat3d[i, :, :])
+    return mat2d
+end
+
+
+"""
+    Cuplus3dj!(mat2d, mat3d, j) # use sum!
+GPU version of plus3dj!
+"""
+Base.@propagate_inbounds function Cuplus3dj!(
+    mat2d::AbstractMatrix,
+    mat3d::AbstractArray{<:Any, 3},
+    j::Int,
+)
+    @boundscheck (size(mat2d, 1) == size(mat3d, 1) || throw("size1"))
+    @boundscheck (size(mat2d, 2) == size(mat3d, 3) || throw("size3"))
+    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
+    mat2d .+= (@view mat3d[:, j, :])
+    return mat2d
+end
+
+
+"""
+    Cuplus3dk!(mat2d, mat3d, k)
+GPU version of plus3dk!
+"""
+Base.@propagate_inbounds function Cuplus3dk!(
+    mat2d::AbstractMatrix,
+    mat3d::AbstractArray{<:Any, 3},
+    k::Int,
+)
+    @boundscheck (size(mat2d, 1) == size(mat3d, 1) || throw("size1"))
+    @boundscheck (size(mat2d, 2) == size(mat3d, 2) || throw("size2"))
+    @boundscheck (1 ≤ k ≤ size(mat3d, 3) || throw("bad k"))
+    mat2d .+= (@view mat3d[:, :, k])
+    return mat2d
+end
+
+
+"""
+    Cuscale3dj!(mat2d, mat3d, j, s)
+GPU version of scale3dj!
+"""
+Base.@propagate_inbounds function Cuscale3dj!(
+    mat2d::AbstractMatrix,
+    mat3d::AbstractArray{<:Any, 3},
+    j::Int,
+    s::RealU,
+)
+    @boundscheck (size(mat2d, 1) == size(mat3d, 1) || throw("size1"))
+    @boundscheck (size(mat2d, 2) == size(mat3d, 3) || throw("size3"))
+    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
+    copyto!(mat2d, (@view mat3d[:, j, :]))
+    mat2d .*= s
+    return mat2d
+end
+
+
+"""
+    Cumul3dj!(mat3d, mat2d, j)
+GPU version of mul3dj!
+"""
+Base.@propagate_inbounds function Cumul3dj!(
+    mat3d::AbstractArray{<:Any, 3},
+    mat2d::AbstractMatrix,
+    j::Int,
+)
+    @boundscheck (size(mat3d, 1) == size(mat2d, 1) || throw("size1"))
+    @boundscheck (size(mat3d, 3) == size(mat2d, 2) || throw("size3"))
+    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
+    (@view mat3d[:, j, :]) .*= mat2d
+    return mat3d
+end
+
+
+"""
+    Cucopy3dj!(mat2d, mat3d, j)
+GPU version of copy3dj!
+"""
+Base.@propagate_inbounds function Cucopy3dj!(
+    mat2d::AbstractMatrix,
+    mat3d::AbstractArray{<:Any, 3},
+    j::Int,
+)
+    @boundscheck (size(mat3d, 1) == size(mat2d, 1) || throw("size1"))
+    @boundscheck (size(mat3d, 3) == size(mat2d, 2) || throw("size3"))
+    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
+    copyto!(mat2d, (@view mat3d[:, j, :]))
+    return mat2d
+end
+
+
+"""
+    Curotl90!(B::AbstractMatrix, A::AbstractMatrix)
+In place GPU version of `rotl90`, returning rotation of `A` in `B`.
+"""
+function Curotl90!(B::AbstractMatrix, A::AbstractMatrix)
+    N = size(B, 2)
+    for i = 1:N
+        copyto!((@view B[:, i]), (@view A[i, end:-1:1]))
+    end
+    return B
+end
+
+
+"""
+    Curotr90!(B::AbstractMatrix, A::AbstractMatrix)
+In place GPU version of `rotr90`, returning rotation of `A` in `B`.
+"""
+function Curotr90!(B::AbstractMatrix, A::AbstractMatrix)
+    N = size(B, 2)
+    for i = 1:N
+        copyto!((@view B[:, N-i+1]), (@view A[i, :]))
+    end
+    return B
+end
+
+
+"""
+    Curot180!(B::AbstractMatrix, A::AbstractMatrix)
+In place GPU version of `rot180`, returning rotation of `A` in `B`.
+"""
+function Curot180!(B::AbstractMatrix, A::AbstractMatrix)
+    N = size(B, 2)
+    for i = 1:N
+        copyto!((@view B[N-i+1, end:-1:1]), (@view A[i, :]))
+    end
+    return B
+end
+
+
+"""
+    Curot_f90!(output, img, m)
+In-place GPU version of rotating an image by 90/180/270 degrees
+"""
+function Curot_f90!(output::AbstractMatrix, img::AbstractMatrix, m::Int)
+    if m == 0
+        output .= img
+    elseif m == 1
+        Curotl90!(output, img)
+    elseif m == 2
+        Curot180!(output, img)
+    elseif m == 3
+        Curotr90!(output, img)
+    else
+        throw("invalid m!")
+    end
+    return output
+end
diff --git a/src/Cuinterp.jl b/src/Cuinterp.jl
new file mode 100644
index 0000000..cd746c5
--- /dev/null
+++ b/src/Cuinterp.jl
@@ -0,0 +1,40 @@
+# Cuinterp.jl
+# Perform linear interpolation for GPU arrays
+using CUDA
+export Cuinterp
+
+struct Cuinterp{T}
+    src::CuArray{T}
+    tex::CuTexture{T}
+    query::Function
+    query!::Function
+    function Cuinterp(src::AbstractArray;
+                      interpmeth::Symbol = :linear,
+                      extrapolation::Symbol = :zero)
+        interp_dict = [:nearest, :linear, :cubic]
+        interp_func = [CUDA.NearestNeighbour(), CUDA.LinearInterpolation(), CUDA.CubicInterpolation()]
+        extra_dict = [:zero, :replicate]
+        extra_func = [CUDA.ADDRESS_MODE_BORDER, CUDA.ADDRESS_MODE_CLAMP]
+        @assert interpmeth in interp_dict || throw("invalid interp method!")
+        @assert extrapolation in extra_dict || throw("invalid extrapolation!")
+        T = eltype(src)
+        src = typeof(src) <: CuArray ? src : CuArray(src)
+        tex = CuTexture(CuTextureArray(src);
+                        interpolation = interp_func[findall(x->x==interpmeth, interp_dict)[1]],
+                        address_mode = extra_func[findall(x->x==extrapolation, extra_dict)[1]])
+        function query(queryidx::Array)
+            queryidx = CuArray(queryidx)
+            dst = CuArray{T}(undef, size(queryidx))
+            broadcast!(dst, queryidx, Ref(tex)) do idx, tex
+                tex[idx...]
+            end
+            return dst
+        end
+        function query!(dst::CuArray, queryidx::CuArray)
+            broadcast!(dst, queryidx, Ref(tex)) do idx, tex
+                tex[idx...]
+            end
+        end
+        new{T}(src, tex, query, query!)
+    end
+end
diff --git a/src/Cuplan-fft.jl b/src/Cuplan-fft.jl
new file mode 100644
index 0000000..5953ad4
--- /dev/null
+++ b/src/Cuplan-fft.jl
@@ -0,0 +1,122 @@
+# Cuplan-fft.jl
+
+export CuPlanFFT
+import AbstractFFTs
+import FFTW
+using CUDA
+
+"""
+    CuPlanFFT{T,Tf,Ti}( ; nx::Int, nz::Int, px::Int, pz::Int, T::DataType)
+Struct for storing work arrays and factors for 2D convolution for one thread.
+Each PSF is `px × pz`
+- `T` datatype of work arrays (subtype of `AbstractFloat`)
+- `nx::Int = 128` (`ny` implicitly the same)
+- `nz::Int = nx` image size is `[nx,nx,nz]`
+- `px::Int = 1`
+- `pz::Int = px` (PSF size)
+- `padsize::Tuple` : `(padup, paddown, padleft, padright)`
+- `workmat [nx+padup+paddown, nz+padleft+padright]` 2D padded image for FFT convolution
+- `workvecx [nx+padup+paddown,]`: 1D work vector
+- `workvecz [nz+padleft+padright,]`: 1D work vector
+- `img_compl [nx+padup+paddown, nz+padleft+padright]`: 2D [complex] padded image for FFT
+- `ker_compl [nx+padup+paddown, nz+padleft+padright]`: 2D [complex] padded image for FFT
+- `fft_plan::Tf` plan for doing FFT; see `plan_fft!`
+- `ifft_plan::Ti` plan for doing IFFT; see `plan_ifft!`
+"""
+struct CuPlanFFT{T, A1, A2, C2, Tf, Ti}
+    nx::Int
+    nz::Int
+    px::Int
+    pz::Int
+    padsize::NTuple{4, Int}
+    workmat::A2
+    workvecx::A1
+    workvecz::A1
+    img_compl::C2
+    ker_compl::C2
+    fft_plan::Tf # Union{AbstractFFTs.ScaledPlan, FFTW.cFFTWPlan}
+    ifft_plan::Ti # Union{AbstractFFTs.ScaledPlan, FFTW.cFFTWPlan}
+
+    function CuPlanFFT( ;
+        nx::Int = 128,
+        nz::Int = nx,
+        px::Int = 1,
+        pz::Int = px,
+        T::DataType = Float32)
+
+        T <: AbstractFloat || throw("invalid T=$T")
+        padup = _padup(nx, px)
+        paddown = _paddown(nx, px)
+        padleft = _padleft(nz, pz)
+        padright = _padright(nz, pz)
+        padsize = (padup, paddown, padleft, padright)
+
+        workmat = CuArray(Matrix{T}(undef, nx+padup+paddown, nz+padleft+padright))
+        workvecx = CuArray(Vector{T}(undef, nx+padup+paddown))
+        workvecz = CuArray(Vector{T}(undef, nz+padleft+padright))
+        vectype = typeof(workvecz)
+        mattype = typeof(workmat)
+        # complex padimg
+        img_compl = CuArray(Matrix{Complex{T}}(undef, nx+padup+paddown, nz+padleft+padright))
+        # complex kernel
+        ker_compl = CuArray(Matrix{Complex{T}}(undef, nx+padup+paddown, nz+padleft+padright))
+        compmattype = typeof(img_compl)
+        fft_plan = plan_fft!(convert(compmattype, ker_compl))
+        ifft_plan = plan_ifft!(convert(compmattype, ker_compl))
+        Tf = typeof(fft_plan)
+        Ti = typeof(ifft_plan)
+
+        new{T, vectype, mattype, compmattype, Tf, Ti}(
+            nx,
+            nz,
+            px,
+            pz,
+            padsize,
+            workmat,
+            workvecx,
+            workvecz,
+            img_compl,
+            ker_compl,
+            fft_plan,
+            ifft_plan,
+        )
+    end
+end
+
+
+"""
+    show(io::IO, ::MIME"text/plain", plan::CuPlanFFT)
+"""
+function Base.show(io::IO, ::MIME"text/plain", plan::CuPlanFFT{T, A1, A2, C2, Tf, Ti}) where {T, A1, A2, C2, Tf, Ti}
+    t = typeof(plan)
+    println(io, t)
+    for f in (:nx, :nz, :px, :pz, :padsize)
+        p = getproperty(plan, f)
+        t = typeof(p)
+        println(io, " ", f, "::", t, " ", p)
+    end
+    for f in (:workmat, :workvecx, :workvecz, :img_compl, :ker_compl, :fft_plan, :ifft_plan)
+        p = getproperty(plan, f)
+        println(io, " ", f, ":", " ", summary(p))
+    end
+    println(io, " (", sizeof(plan), " bytes)")
+end
+
+
+"""
+    show(io::IO, mime::MIME"text/plain", vp::Vector{<:CuPlanFFT})
+"""
+function Base.show(io::IO, mime::MIME"text/plain", vp::Vector{<: CuPlanFFT})
+    t = typeof(vp)
+    println(io, length(vp), "-element ", t, " with N=", vp[1].nx)
+#   show(io, mime, vp[1])
+end
+
+
+"""
+    sizeof(::CuPlanFFT)
+Show size in bytes of `CuPlanFFT` object.
+"""
+function Base.sizeof(ob::T) where {T <: CuPlanFFT}
+    sum(f -> sizeof(getfield(ob, f)), fieldnames(typeof(ob)))
+end
diff --git a/src/Cuplan-rotate.jl b/src/Cuplan-rotate.jl
new file mode 100644
index 0000000..f0f2319
--- /dev/null
+++ b/src/Cuplan-rotate.jl
@@ -0,0 +1,89 @@
+# Cuplan-rotate.jl
+
+using CUDA
+export CuPlanRotate
+
+"""
+    CuPlanRotate{T, G2}
+Struct for storing GPU work arrays and factors for 2D square image rotation
+- `T` datatype of work arrays (default `Float32`)
+- `nx::Int` image size
+- `padsize::Int` : pad each side of image by this much
+- `center::Int`: center of the (padded) image
+- `xaxis [nx + 2 * padsize, nx + 2 * padsize]` the index array for x axis
+- `yaxis [nx + 2 * padsize, nx + 2 * padsize]` the index array for y axis
+- `workmat1 [nx + 2 * padsize, nx + 2 * padsize]` padded work matrix
+- `workmat2 [nx + 2 * padsize, nx + 2 * padsize]` padded work matrix
+- `workmat3 [nx + 2 * padsize, nx + 2 * padsize]` padded work matrix
+- `padded_src`: matrix storing the padded input image
+- `query!`: function that returns the value of query indices
+"""
+struct CuPlanRotate{T, G2}
+    nx::Int
+    padsize::Int
+    center::T
+    xaxis::G2
+    yaxis::G2
+    workmat1::G2
+    workmat2::G2
+    workmat3::G2
+    padded_src::G2
+    query!::Function
+
+    function CuPlanRotate(
+        nx::Int ;
+        T::DataType = Float32,
+    )
+        # only support the case that the image is square
+        padsize = ceil(Int, 1 + nx * sqrt(2)/2 - nx / 2)
+
+        px = nx + 2 * padsize
+        center = (1 + px) / 2
+        axis = T.((1:px) .- center)
+        idxarray = CuArray([(i, j) for i in axis, j in axis])
+        xaxis = first.(idxarray)
+        yaxis = last.(idxarray)
+        padded_src = CuArray{T, 2}(undef, px, px)
+        workmat1 = CuArray{T, 2}(undef, px, px)
+        workmat2 = CuArray{T, 2}(undef, px, px)
+        workmat3 = CuArray{T, 2}(undef, px, px)
+        G2 = typeof(workmat1)
+
+        function query!(dst::CuArray, tex::CuTexture, queryidx::CuArray)
+            broadcast!(dst, queryidx, Ref(tex)) do idx, tex
+                tex[idx...]
+            end
+        end
+        new{T, G2}(nx, padsize, center, xaxis, yaxis, workmat1,
+                    workmat2, workmat3, padded_src, query!)
+    end
+end
+
+
+
+"""
+    show(io::IO, ::MIME"text/plain", plan::CuPlanRotate)
+"""
+function Base.show(io::IO, ::MIME"text/plain", plan::CuPlanRotate{T, G2}) where {T, G2}
+    t = typeof(plan)
+    println(io, t)
+    for f in (:nx, :padsize, :center)
+        p = getproperty(plan, f)
+        t = typeof(p)
+        println(io, " ", f, "::", t, " ", p)
+    end
+    for f in (:xaxis, :yaxis, :workmat1, :workmat2, :workmat3)
+        p = getproperty(plan, f)
+        println(io, " ", f, ":", " ", summary(p))
+    end
+    println(io, " (", sizeof(plan), " bytes)")
+end
+
+
+"""
+    sizeof(::CuPlanRotate)
+Show size in bytes of `CuPlanRotate` object.
+"""
+function Base.sizeof(ob::T) where {T <: CuPlanRotate}
+    sum(f -> sizeof(getfield(ob, f)), fieldnames(typeof(ob)))
+end
diff --git a/src/Cuproject.jl b/src/Cuproject.jl
new file mode 100644
index 0000000..698744e
--- /dev/null
+++ b/src/Cuproject.jl
@@ -0,0 +1,108 @@
+# Cuproject.jl
+
+export Cuproject, Cuproject!
+
+"""
+    Cuproject!(view, plan, image, viewidx)
+GPU version of SPECT projection of `image` into a single `view` with index `viewidx`.
+The `view` must be pre-allocated but need not be initialized to zero.
+"""
+function Cuproject!(
+    view::AbstractMatrix{<:RealU},
+    image::AbstractArray{<:RealU, 3},
+    plan::CuSPECTplan,
+    viewidx::Int,
+)
+
+    for z = 1:plan.imgsize[3]
+        # rotate images
+        Cuimrotate!(
+            (@view plan.imgr[:, :, z]),
+            (@view image[:, :, z]),
+            plan.viewangle[viewidx],
+            plan.planrot)
+        # rotate mumap
+        Cuimrotate!(
+            (@view plan.mumapr[:, :, z]),
+            (@view plan.mumap[:, :, z]),
+            plan.viewangle[viewidx],
+            plan.planrot)
+    end
+
+    for y = 1:plan.imgsize[2] # 1:ny
+
+        Cugen_attenuation!(plan, y)
+        # apply depth-dependent attenuation
+        Cumul3dj!(plan.imgr, plan.exp_mumapr, y)
+
+        Cufft_conv!(plan.add_view,
+                   (@view plan.imgr[:, y, :]),
+                   (@view plan.psfs[:, :, y, viewidx]),
+                    plan.planpsf)
+
+        view .+= plan.add_view
+    end
+
+    return view
+end
+
+
+"""
+    Cuproject!(views, image, plan; viewlist)
+Project `image` into multiple `views` with indexes `index` (defaults to `1:nview`).
+The 3D `views` array must be pre-allocated, but need not be initialized.
+"""
+function Cuproject!(
+    views::AbstractArray{<:RealU,3},
+    image::AbstractArray{<:RealU,3},
+    plan::CuSPECTplan;
+    viewlist::AbstractVector{<:Int} = 1:plan.nview, # all views
+)
+
+    # loop over each view index
+    for viewidx in viewlist
+        Cuproject!((@view views[:,:,viewidx]), image, plan, viewidx)
+    end
+
+    return views
+end
+
+
+"""
+    views = Cuproject(image, plan ; kwargs...)
+GPU version of a convenience method for SPECT forward projector that allocates and returns views.
+"""
+function Cuproject(
+    image::AbstractArray{<:RealU,3},
+    plan::CuSPECTplan;
+    kwargs...,
+)
+    views = CuArray{plan.T}(undef, plan.imgsize[1], plan.imgsize[3], plan.nview)
+    Cuproject!(views, image, plan; kwargs...)
+    return views
+end
+
+
+"""
+    views = Cuproject(image, mumap, psfs, dy; kwargs...)
+GPU version of a convenience method for SPECT forward projector that does all allocation
+including initializing `plan`.
+
+In
+* `image` : 3D array `(nx,ny,nz)`
+* `mumap` : `(nx,ny,nz)` 3D attenuation map, possibly zeros()
+* `psfs` : 4D PSF array
+* `dy::RealU` : pixel size
+"""
+function Cuproject(
+    image::AbstractArray{<:RealU, 3},
+    mumap::AbstractArray{<:RealU, 3}, # (nx,ny,nz) 3D attenuation map
+    psfs::AbstractArray{<:RealU, 4}, # (px,pz,ny,nview)
+    dy::RealU;
+    kwargs...,
+)
+    size(mumap) == size(image) || throw(DimensionMismatch("image/mumap size"))
+    size(image,2) == size(psfs,3) || throw(DimensionMismatch("image/psfs size"))
+    plan = CuSPECTplan(mumap, psfs, dy; kwargs...)
+    return Cuproject(image, plan; kwargs...)
+end
diff --git a/src/Cupsf-gauss.jl b/src/Cupsf-gauss.jl
new file mode 100644
index 0000000..ae64845
--- /dev/null
+++ b/src/Cupsf-gauss.jl
@@ -0,0 +1,55 @@
+# Cupsf-gauss.jl
+
+export Cupsf_gauss
+
+"""
+    Cupsf_gauss( ; ny, px, pz, fwhm_start, fwhm_end, fwhm, fwhm_x, fwhm_z, T)
+
+Create depth-dependent Gaussian PSFs
+having specified full-width half-maximum (FHWM) values.
+
+# Options
+- 'ny::Int = 128'
+- 'px::Int = 11' (should be odd)
+- 'pz::Int = px' (should be odd)
+- 'fwhm_start::Real = 1'
+- 'fwhm_end::Real = 4'
+- 'fwhm::AbstractVector{<:Real} = LinRange(fwhm_start, fwhm_end, ny)'
+- 'fwhm_x::AbstractVector{<:Real} = fwhm,
+- 'fwhm_z::AbstractVector{<:Real} = fwhm_x'
+- 'T::DataType == Float32'
+
+Returned `psf` is `[px, pz, ny]` where each PSF sums to 1.
+"""
+function Cupsf_gauss( ;
+    ny::Int = 128,
+    px::Int = 11,
+    pz::Int = px,
+    fwhm_start::Real = 1,
+    fwhm_end::Real = 4,
+    fwhm::AbstractVector{<:Real} = LinRange(fwhm_start, fwhm_end, ny),
+    fwhm_x::AbstractVector{<:Real} = fwhm,
+    fwhm_z::AbstractVector{<:Real} = fwhm_x,
+    T::DataType = Float32,
+)
+    isodd(px) || @warn("even px = $px ?")
+    isodd(pz) || @warn("even pz = $pz ?")
+    psf = CuArray(zeros(T, px, pz, ny))
+
+    rx = (-(px-1)÷2):((px-1)÷2)
+    rz = (-(pz-1)÷2):((pz-1)÷2)
+    for iy in 1:ny # depth-dependent blur
+        psf[:,:,iy] = Cupsf_gauss(rx, fwhm_x[iy]) * Cupsf_gauss(rz, fwhm_z[iy])'
+    end
+    return psf
+end
+
+function Cupsf_gauss(r::AbstractVector, fwhm::Real)
+   if fwhm == 0
+        any(==(0), r) || throw("must have some r=0 if fwhm=0")
+        return r .== 0 # Kronecker impulse
+   end
+   σ = fwhm / sqrt(log(256)) # FWHM to Gaussian σ
+   psf = @. exp(-0.5 * abs2(r / σ))
+   return psf / sum(psf)
+end
diff --git a/src/Curotate.jl b/src/Curotate.jl
new file mode 100644
index 0000000..03c3214
--- /dev/null
+++ b/src/Curotate.jl
@@ -0,0 +1,52 @@
+# Curotate.jl
+
+using CUDA
+export Cuimrotate!
+
+"""
+    Cuimrotate!(outarr, inarr, θ, plan)
+GPU version of imrotate.
+- `outarr`: output array
+- `inarr`: input array
+- `theta`: rotate angle (in rad)
+- `plan`: rotate plan, see "Cuplan-rotate.jl"
+Note: functions "Curot_f90!" and "Cupadzero!" are from "Cuhelper.jl"
+"""
+function Cuimrotate!(outarr::CuArray,
+                     inarr::CuArray,
+                     θ::Real,
+                     plan::CuPlanRotate)
+    if mod(θ, 2π) ≈ 0
+        outarr .= inarr
+        return outarr
+    end
+    m = mod(floor(Int, 0.5 + θ/(π/2)), 4)
+    if θ ≈ m * (π/2)
+        Curot_f90!(outarr, inarr, m)
+    else
+        mod_theta = θ - m * (π/2)
+        sin_θ, cos_θ = sincos(mod_theta)
+        Cupadzero!(plan.padded_src, inarr,
+                    (plan.padsize, plan.padsize, plan.padsize, plan.padsize))
+        Curot_f90!(plan.workmat1, plan.padded_src, m)
+        copyto!(plan.padded_src, plan.workmat1)
+        # multiply with the rotation matrix [cos(theta) sin(theta); -sin(theta) cos(theta)]
+        broadcast!(*, plan.workmat1, plan.xaxis, cos_θ)
+        broadcast!(*, plan.workmat2, plan.yaxis, sin_θ)
+        broadcast!(+, plan.workmat1, plan.workmat1, plan.workmat2)
+        broadcast!(*, plan.workmat2, plan.xaxis, sin_θ)
+        broadcast!(*, plan.workmat3, plan.yaxis, cos_θ)
+        broadcast!(-, plan.workmat2, plan.workmat3, plan.workmat2)
+        broadcast!(+, plan.workmat1, plan.workmat1, plan.center)
+        broadcast!(+, plan.workmat2, plan.workmat2, plan.center)
+        # line 43 - 46 surely allocates some GPU memory
+        allocated_idx = tuple.(plan.workmat1, plan.workmat2)
+        tex = CuTexture(CuTextureArray(plan.padded_src);
+                        interpolation = CUDA.LinearInterpolation(), # linear interpolation
+                        address_mode = CUDA.ADDRESS_MODE_BORDER) # filling out-of-boundary values with zeros
+        plan.query!(plan.workmat1, tex, allocated_idx) # give interpolated values
+        copyto!(outarr, (@view plan.workmat1[plan.padsize .+ (1:plan.nx),
+                                             plan.padsize .+ (1:plan.nx)]))
+    end
+    return outarr
+end

From 62fd0b9e1b357e633c9ee49224dc3a1522885718 Mon Sep 17 00:00:00 2001
From: Zongyu Li <54018575+ZongyuLi-umich@users.noreply.github.com>
Date: Thu, 19 May 2022 13:09:53 -0400
Subject: [PATCH 2/5] test code

---
 test/Cubackproject.jl |  68 ++++++
 test/Cufftconv.jl     |  54 +++++
 test/Cuhelper.jl      | 224 ++++++++++++++++++++
 test/Cuproject.jl     |  33 +++
 test/Cupsf-gauss.jl   |  28 +++
 test/Curotate.jl      |  42 ++++
 test/Curuntests.jl    |  18 ++
 test/Manifest.toml    | 467 ++++++++++++++++++++++++++++++++++++++++++
 test/Project.toml     |   9 +
 9 files changed, 943 insertions(+)
 create mode 100644 test/Cubackproject.jl
 create mode 100644 test/Cufftconv.jl
 create mode 100644 test/Cuhelper.jl
 create mode 100644 test/Cuproject.jl
 create mode 100644 test/Cupsf-gauss.jl
 create mode 100644 test/Curotate.jl
 create mode 100644 test/Curuntests.jl
 create mode 100644 test/Manifest.toml
 create mode 100644 test/Project.toml

diff --git a/test/Cubackproject.jl b/test/Cubackproject.jl
new file mode 100644
index 0000000..dc304cc
--- /dev/null
+++ b/test/Cubackproject.jl
@@ -0,0 +1,68 @@
+# Cubackproject.jl
+
+using SPECTrecon: SPECTplan, backproject!
+using LinearAlgebra: dot
+using Random: seed!
+
+@testset "Cubackproject" begin
+    T = Float32
+    nx = 8; ny = nx
+    nz = 6
+    nview = 7
+
+    mumap = rand(T, nx, ny, nz)
+
+    px = 5
+    pz = 3
+    psfs = rand(T, px, pz, ny, nview)
+    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2]) # symmetrize
+    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])
+
+    dy = T(4.7952)
+    plan = SPECTplan(mumap, psfs, dy; T, interpmeth = :two, mode = :fast)
+    image = zeros(T, nx, ny, nz) # images must be initialized to zero
+    views = rand(T, nx, nz, nview)
+    Cuimage = CuArray(zeros(T, nx, ny, nz))
+    Cuviews = CuArray(views)
+
+    Cumumap = CuArray(mumap)
+    Cupsfs = CuArray(psfs)
+    Cuplan = CuSPECTplan(Cumumap, Cupsfs, dy; T)
+
+    backproject!(image, views, plan)
+    Cubackproject!(Cuimage, Cuviews, Cuplan)
+    @test isapprox(Array(Cuimage), image; rtol = 0.6) # this is very *large* difference
+end
+
+
+@testset "adjoint-Cuproject" begin
+    seed!(0)
+    T = Float32
+    nx = 8; ny = nx
+    nz = 6
+    nview = 7
+
+    mumap = rand(T, nx, ny, nz)
+
+    px = 5
+    pz = 3
+    psfs = rand(T, px, pz, ny, nview)
+    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2]) # symmetrize
+    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])
+
+    dy = T(4.7952)
+
+    Cumumap = CuArray(mumap)
+    Cupsfs = CuArray(psfs)
+    Cuplan = CuSPECTplan(Cumumap, Cupsfs, dy; T)
+
+    image = CuArray(rand(T, nx, ny, nz))
+    backimage = CuArray(zeros(T, nx, ny, nz))
+    views = CuArray(rand(T, nx, nz, nview))
+    forviews = CuArray(zeros(T, nx, nz, nview))
+
+    Cuproject!(forviews, image, Cuplan)
+    Cubackproject!(backimage, views, Cuplan)
+    @test isapprox(dot(forviews, views), dot(backimage, image); rtol = 0.05)
+    # have to put rtol = 0.05 to pass the test
+end
diff --git a/test/Cufftconv.jl b/test/Cufftconv.jl
new file mode 100644
index 0000000..545fcaa
--- /dev/null
+++ b/test/Cufftconv.jl
@@ -0,0 +1,54 @@
+# Cufftconv.jl
+# test FFT-based convolution and
+# adjoint consistency for FFT convolution methods on very small case
+
+using SPECTrecon: fft_conv!, fft_conv_adj!, plan_psf
+using LinearAlgebra: dot
+
+
+@testset "Cufftconv!" begin
+    T = Float32
+    nx = 12
+    nz = 8
+    px = 7
+    pz = 7
+    for i = 1:4
+        img = randn(T, nx, nz)
+        ker = rand(T, px, pz)
+        ker_sym = ker .+ reverse(ker, dims=:)
+        ker_sym /= sum(ker_sym)
+        out = similar(img)
+        Cuout = CuArray(out)
+        Cuimg = CuArray(img)
+        Cuker_sym = CuArray(ker_sym)
+
+        plan = plan_psf( ; nx, nz, px, pz, T, nthread = 1)[1]
+        fft_conv!(out, img, ker_sym, plan)
+        Cuplan = CuPlanFFT( ; nx, nz, px, pz, T)
+        Cufft_conv!(Cuout, Cuimg, Cuker_sym, Cuplan)
+        @test isapprox(out, Array(Cuout))
+    end
+end
+
+
+@testset "adjoint-Cufftconv!" begin
+    nx = 20
+    nz = 14
+    px = 5
+    pz = 5
+    T = Float32
+    for i = 1:4 # test with different kernels
+        x = cu(randn(T, nx, nz))
+        out_x = similar(x)
+        y = cu(randn(T, nx, nz))
+        out_y = similar(y)
+        ker = cu(rand(T, px, pz))
+        ker = ker .+ reverse(reverse(ker, dims=1), dims=2)
+        ker /= sum(ker)
+        plan = CuPlanFFT( ; nx, nz, px, pz, T)
+        Cufft_conv!(out_x, x, ker, plan)
+        Cufft_conv_adj!(out_y, y, ker, plan)
+
+        @test dot(Array(out_x), Array(y)) ≈ dot(Array(out_y), Array(x))
+    end
+end
diff --git a/test/Cuhelper.jl b/test/Cuhelper.jl
new file mode 100644
index 0000000..897fc6b
--- /dev/null
+++ b/test/Cuhelper.jl
@@ -0,0 +1,224 @@
+# Cuhelper.jl
+
+using SPECTrecon: rot180!, rot_f90!, rotl90!, rotr90!
+using OffsetArrays
+using ImageFiltering: BorderArray, Fill, Pad
+
+
+@testset "Cupadzero!" begin
+    T = Float32
+    x = randn(T, 7, 5)
+    Cux = CuArray(x)
+    y = randn(T, 3, 3)
+    Cuy = CuArray(y)
+    Cupadzero!(Cux, Cuy, (2, 2, 1, 1))
+    z = OffsetArrays.no_offset_view(BorderArray(y, Fill(0, (2, 1), (2, 1))))
+    x_cpu = Array(Cux)
+    @test x_cpu == z
+end
+
+
+@testset "Cupadrepl!" begin
+    T = Float32
+    x = randn(T, 10, 9)
+    Cux = CuArray(x)
+    y = randn(T, 5, 4)
+    Cuy = CuArray(y)
+    Cupadrepl!(Cux, Cuy, (1, 4, 3, 2)) # up, down, left, right
+    z = OffsetArrays.no_offset_view(BorderArray(y, Pad(:replicate, (1, 3), (4, 2)))) # up, left, down, right
+    x_cpu = Array(Cux)
+    @test x_cpu == z
+end
+
+
+@testset "Cupad2sizezero!" begin
+    T = Float32
+    x = reshape(Int16(1):Int16(15), 5, 3)
+    Cux = CuArray(x)
+    padsize = (8, 6)
+    z = CuArray(randn(T, padsize))
+    Cupad2sizezero!(z, Cux, padsize)
+    tmp = OffsetArrays.no_offset_view(BorderArray(x, Fill(0, (2, 2), (1, 1))))
+    z_cpu = Array(z)
+    @test tmp == z_cpu
+end
+
+
+@testset "Cuplus1di!" begin
+    T = Float32
+    x = randn(T, 4, 9)
+    Cux = CuArray(x)
+    v = randn(T, 9)
+    Cuv = CuArray(v)
+    y = x[2, :] .+ v
+    Cuplus1di!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu[2, :] == y
+end
+
+
+@testset "Cuplus1dj!" begin
+    T = Float32
+    x = randn(T, 9, 4)
+    Cux = CuArray(x)
+    v = randn(T, 9)
+    Cuv = CuArray(v)
+    y = x[:, 2] .+ v
+    Cuplus1dj!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu[:, 2] == y
+end
+
+
+@testset "Cuplus2di!" begin
+    T = Float32
+    x = randn(9)
+    Cux = CuArray(x)
+    v = randn(4, 9)
+    Cuv = CuArray(v)
+    y = x .+ v[2, :]
+    Cuplus2di!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Cuplus2dj!" begin
+    T = Float32
+    x = randn(T, 9)
+    Cux = CuArray(x)
+    v = randn(T, 9, 4)
+    Cuv = CuArray(v)
+    y = x .+ v[:, 2]
+    Cuplus2dj!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Cuplus3di!" begin
+    T = Float32
+    x = randn(T, 9, 7)
+    Cux = CuArray(x)
+    v = randn(T, 4, 9, 7)
+    Cuv = CuArray(v)
+    y = x .+ v[2, :, :]
+    Cuplus3di!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Cuplus3dj!" begin
+    T = Float32
+    x = randn(T, 9, 7)
+    Cux = CuArray(x)
+    v = randn(T, 9, 4, 7)
+    Cuv = CuArray(v)
+    y = x .+ v[:, 2, :]
+    Cuplus3dj!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Cuplus3dk!" begin
+    T = Float32
+    x = randn(T, 9, 7)
+    Cux = CuArray(x)
+    v = randn(T, 9, 7, 4)
+    Cuv = CuArray(v)
+    y = x .+ v[:, :, 2]
+    Cuplus3dk!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Cuscale3dj!" begin
+    T = Float32
+    x = randn(T, 9, 7)
+    Cux = CuArray(x)
+    v = randn(T, 9, 4, 7)
+    Cuv = CuArray(v)
+    s = -0.5
+    y = s * v[:, 2, :]
+    Cuscale3dj!(Cux, Cuv, 2, s)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Cumul3dj!" begin
+    T = Float32
+    x = randn(T, 9, 4, 7)
+    Cux = CuArray(x)
+    v = randn(T, 9, 7)
+    Cuv = CuArray(v)
+    y = x[:,2,:] .* v
+    Cumul3dj!(Cux, Cuv, 2)
+    x_cpu = Array(Cux[:,2,:])
+    @test x_cpu == y
+end
+
+
+@testset "Cucopy3dj!" begin
+    T = Float32
+    x = randn(T, 9, 7)
+    Cux = CuArray(x)
+    v = randn(T, 9, 4, 7)
+    Cuv = CuArray(v)
+    y = v[:,2,:]
+    Cucopy3dj!(Cux, Cuv, 2)
+    x_cpu = Array(Cux)
+    @test x_cpu == y
+end
+
+
+@testset "Curotl90!" begin
+    T = Float32
+    N = 20
+    A = rand(T, N, N)
+    B = similar(A, N, N)
+    CuA = CuArray(A)
+    CuB = CuArray(B)
+    Curotl90!(CuB, CuA)
+    rotl90!(B, A)
+    @test isequal(B, Array(CuB))
+end
+
+
+@testset "Curotr90!" begin
+    T = Float32
+    N = 20
+    A = rand(T, N, N)
+    B = similar(A)
+    CuA = CuArray(A)
+    CuB = CuArray(B)
+    Curotr90!(CuB, CuA)
+    rotr90!(B, A)
+    @test isequal(B, Array(CuB))
+end
+
+
+@testset "Curot180!" begin
+    T = Float32
+    N = 20
+    A = rand(T, N, N)
+    B = similar(A)
+    CuA = CuArray(A)
+    CuB = CuArray(B)
+    Curot180!(CuB, CuA)
+    rot180!(B, A)
+    @test isequal(B, Array(CuB))
+end
+
+
+@testset "Curot_f90!" begin
+    T = Float32
+    N = 20
+    A = CuArray(rand(T, N, N))
+    B = CuArray(rand(T, N, N))
+    @test_throws String Curot_f90!(A, B, -1)
+    @test_throws String Curot_f90!(A, B, 4)
+end
diff --git a/test/Cuproject.jl b/test/Cuproject.jl
new file mode 100644
index 0000000..54f72d9
--- /dev/null
+++ b/test/Cuproject.jl
@@ -0,0 +1,33 @@
+# Cuproject.jl
+
+using SPECTrecon: SPECTplan, project!
+
+
+@testset "Cuproject" begin
+    T = Float32
+    nx = 8; ny = nx
+    nz = 6
+    nview = 7
+
+    mumap = rand(T, nx, ny, nz)
+
+    px = 5
+    pz = 3
+    psfs = rand(T, px, pz, ny, nview)
+    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2]) # symmetrize
+    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])
+
+    dy = T(4.7952)
+    plan = SPECTplan(mumap, psfs, dy; T, interpmeth = :two, mode = :fast)
+    x = randn(T, nx, ny, nz)
+    views = zeros(T, nx, nz, nview)
+
+    Cumumap = CuArray(mumap)
+    Cupsfs = CuArray(psfs)
+    Cuplan = CuSPECTplan(Cumumap, Cupsfs, dy; T)
+    Cux = CuArray(x)
+    Cuviews = CuArray(zeros(T, nx, nz, nview))
+    project!(views, x, plan)
+    Cuproject!(Cuviews, Cux, Cuplan)
+    @test isapprox(Array(Cuviews), views; rtol = 1e-2)
+end
diff --git a/test/Cupsf-gauss.jl b/test/Cupsf-gauss.jl
new file mode 100644
index 0000000..8ab03d4
--- /dev/null
+++ b/test/Cupsf-gauss.jl
@@ -0,0 +1,28 @@
+# Cupsf-gauss.jl
+
+
+@testset "Cupsf-gauss" begin
+    psf = @inferred Cupsf_gauss()
+    @test Array(psf) isa Array{Float32,3}
+
+    ny = 4
+    px = 5
+    pz = 3
+    psf = @inferred Cupsf_gauss(; ny, px, pz, fwhm = zeros(ny))
+    tmp = zeros(px,pz)
+    tmp[(end+1)÷2,(end+1)÷2] = 1 # Kronecker impulse
+    tmp = repeat(tmp, 1, 1, ny)
+    @test Array(psf) == tmp
+
+    ny = 4
+    px = 5
+    pz = 3
+    psf = @inferred Cupsf_gauss(; ny, px, pz,
+        fwhm_x = fill(Inf, ny),
+        fwhm_z = zeros(ny),
+    )
+    tmp = zeros(px,pz)
+    tmp[:,(end+1)÷2] .= 1/px # wide in x
+    tmp = repeat(tmp, 1, 1, ny)
+    @test Array(psf) ≈ tmp
+end
diff --git a/test/Curotate.jl b/test/Curotate.jl
new file mode 100644
index 0000000..a8dfbc5
--- /dev/null
+++ b/test/Curotate.jl
@@ -0,0 +1,42 @@
+# Curotatez.jl
+
+using SPECTrecon: imrotate!, imrotate_adj!, plan_rotate
+using LinearAlgebra: dot
+
+
+@testset "rotate-gpu" begin
+    nx = 20
+    θ_list = (1:23) / 12 * π
+    T = Float32
+    image2 = rand(T, nx, nx)
+    plan = plan_rotate(nx; T, method = :two)[1]
+    result = similar(image2)
+    Cuimage2 = CuArray(image2)
+    Curesult = similar(Cuimage2)
+    Cuplan = CuPlanRotate(nx; T)
+    for θ in θ_list
+        imrotate!(result, image2, θ, plan)
+        Cuimrotate!(Curesult, Cuimage2, θ, Cuplan)
+        @test isapprox(result, Array(Curesult); rtol = 2e-3)
+        imrotate_adj!(result, image2, θ, plan)
+        Cuimrotate!(Curesult, Cuimage2, -θ, Cuplan) # rotate the image by "-angle"
+        @test isapprox(result, Array(Curesult); rtol = 2e-1)
+    end
+end
+
+
+@testset "Cuadjoint-rotate" begin
+    nx = 20
+    θ_list = (1:23) / 12 * π
+    T = Float32
+    x = CuArray(rand(T, nx, nx))
+    y = CuArray(rand(T, nx, nx))
+    out_x = similar(x)
+    out_y = similar(y)
+    Cuplan = CuPlanRotate(nx; T)
+    for θ in θ_list
+        Cuimrotate!(out_x, x, θ, Cuplan)
+        Cuimrotate!(out_y, y, -θ, Cuplan)
+        @test isapprox(dot(out_x, y), dot(out_y, x); rtol = 1e-2)
+    end
+end
diff --git a/test/Curuntests.jl b/test/Curuntests.jl
new file mode 100644
index 0000000..4ec6842
--- /dev/null
+++ b/test/Curuntests.jl
@@ -0,0 +1,18 @@
+# Curuntests.jl
+
+include("../src/CuSPECTrecon.jl")
+using Main.CuSPECTrecon
+using Test: @test, @testset, @test_throws, @inferred, detect_ambiguities
+using CUDA
+CUDA.allowscalar(false)
+
+include("Cuhelper.jl")
+include("Curotate.jl")
+include("Cufftconv.jl")
+include("Cupsf-gauss.jl")
+include("Cuproject.jl")
+include("Cubackproject.jl")
+
+@testset "CuSPECTrecon" begin
+    @test isempty(detect_ambiguities(CuSPECTrecon))
+end
diff --git a/test/Manifest.toml b/test/Manifest.toml
new file mode 100644
index 0000000..348275f
--- /dev/null
+++ b/test/Manifest.toml
@@ -0,0 +1,467 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.7.1"
+manifest_format = "2.0"
+
+[[deps.AbstractFFTs]]
+deps = ["ChainRulesCore", "LinearAlgebra"]
+git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
+uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+version = "1.1.0"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.3.3"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[deps.ArrayTools]]
+git-tree-sha1 = "ac94caaccbabd2688e2d1b28ce5d12248e305818"
+uuid = "1dc0ca97-c5ce-4e77-ac6d-c576ac9d7f27"
+version = "0.2.3"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.CatIndices]]
+deps = ["CustomUnitRanges", "OffsetArrays"]
+git-tree-sha1 = "a0f80a09780eed9b1d106a1bf62041c2efc995bc"
+uuid = "aafaddc9-749c-510e-ac4f-586e18779b91"
+version = "0.2.2"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.14.0"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.3"
+
+[[deps.ColorTypes]]
+deps = ["FixedPointNumbers", "Random"]
+git-tree-sha1 = "a985dc37e357a3b22b260a5def99f3530fb415d3"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.11.2"
+
+[[deps.ColorVectorSpace]]
+deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
+git-tree-sha1 = "3f1f500312161f1ae067abe07d13b40f78f32e07"
+uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
+version = "0.9.8"
+
+[[deps.Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.12.8"
+
+[[deps.Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.43.0"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.ComputationalResources]]
+git-tree-sha1 = "52cb3ec90e8a8bea0e62e275ba577ad0f74821f7"
+uuid = "ed09eef8-17a6-5b46-8889-db040fac31e3"
+version = "0.3.2"
+
+[[deps.CustomUnitRanges]]
+git-tree-sha1 = "1a3f97f907e6dd8983b744d2642651bb162a3f7a"
+uuid = "dc8bdbbb-1ca9-579f-8c36-e416f6a65cce"
+version = "1.0.2"
+
+[[deps.DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "cc1a8e22627f33c789ab60b36a9132ac050bbf75"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.12"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[deps.Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.6"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[deps.FFTViews]]
+deps = ["CustomUnitRanges", "FFTW"]
+git-tree-sha1 = "cbdf14d1e8c7c8aacbe8b19862e0179fd08321c2"
+uuid = "4f61f5a4-77b1-5117-aa51-3ab5ef4ef0cd"
+version = "0.3.2"
+
+[[deps.FFTW]]
+deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
+git-tree-sha1 = "505876577b5481e50d089c1c68899dfb6faebc62"
+uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+version = "1.4.6"
+
+[[deps.FFTW_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
+uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
+version = "3.3.10+0"
+
+[[deps.FixedPointNumbers]]
+deps = ["Statistics"]
+git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.8.4"
+
+[[deps.Graphics]]
+deps = ["Colors", "LinearAlgebra", "NaNMath"]
+git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
+uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
+version = "1.1.1"
+
+[[deps.ImageCore]]
+deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
+git-tree-sha1 = "9a5c62f231e5bba35695a20988fc7cd6de7eeb5a"
+uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+version = "0.9.3"
+
+[[deps.ImageFiltering]]
+deps = ["CatIndices", "ComputationalResources", "DataStructures", "FFTViews", "FFTW", "ImageCore", "LinearAlgebra", "OffsetArrays", "Requires", "SparseArrays", "StaticArrays", "Statistics", "TiledIteration"]
+git-tree-sha1 = "79dac52336910325a5675813053b1eee3eb5dcc6"
+uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
+version = "0.6.22"
+
+[[deps.IntelOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2018.0.3+2"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.InterpolationKernels]]
+deps = ["InteractiveUtils", "Printf"]
+git-tree-sha1 = "4c9909e2649ceaeaa2fa2fc262af1e70515a7bee"
+uuid = "16730964-a2ec-11e9-36fa-47a4f82bfac6"
+version = "0.1.3"
+
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "336cc738f03e069ef2cac55a104eb823455dca75"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.4"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.1"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[deps.LazyAlgebra]]
+deps = ["AbstractFFTs", "ArrayTools", "FFTW", "LinearAlgebra", "MayOptimize", "Printf", "Requires", "SparseArrays", "StructuredArrays", "ZippedArrays"]
+git-tree-sha1 = "2fcaeb99c3525d40f7382431add5c98fe3166249"
+uuid = "aeaa251e-9d76-11e9-311b-83636bf44738"
+version = "0.2.4"
+
+[[deps.LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LinearInterpolators]]
+deps = ["ArrayTools", "InteractiveUtils", "InterpolationKernels", "LazyAlgebra", "LinearAlgebra", "Printf", "SparseArrays", "TwoDimensional"]
+git-tree-sha1 = "fd1c148e6571a58dde76880bf85c6be520c9cd5c"
+uuid = "2fc109c4-9d77-11e9-32ff-e5c650334c16"
+version = "0.1.6"
+
+[[deps.LinearMaps]]
+deps = ["LinearAlgebra", "SparseArrays", "Statistics"]
+git-tree-sha1 = "9e9c3fa4de0b4f146d97eed3485711928789865b"
+uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
+version = "3.6.2"
+
+[[deps.LinearMapsAA]]
+deps = ["LinearAlgebra", "LinearMaps", "Requires", "SparseArrays", "Test"]
+git-tree-sha1 = "efb70fe4d954f6f5c9e310e1940e1c4663205144"
+uuid = "599c1a8e-b958-11e9-0d14-b1e6b2ecea07"
+version = "0.9.0"
+
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.15"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.MKL_jll]]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "e595b205efd49508358f7dc670a940c790204629"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2022.0.0+0"
+
+[[deps.MappedArrays]]
+git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
+uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
+version = "0.4.1"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.MayOptimize]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "1be2f05530c5f05d6fba91057f2efa3f2338c507"
+uuid = "049513e6-5a5e-4281-9d48-334ebac8954e"
+version = "0.3.1"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[deps.MosaicViews]]
+deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
+git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
+uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
+version = "0.3.3"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[deps.NaNMath]]
+git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "0.3.7"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[deps.OffsetArrays]]
+deps = ["Adapt"]
+git-tree-sha1 = "e6c5f47ba51b734a4e264d7183b6750aec459fa0"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.11.1"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[deps.PaddedViews]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "03a7a85b76381a3d04c7a1656039197e70eda03d"
+uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
+version = "0.5.11"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.Reexport]]
+git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.2.2"
+
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.0"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.SPECTrecon]]
+deps = ["AbstractFFTs", "FFTW", "ImageFiltering", "LinearAlgebra", "LinearInterpolators", "OffsetArrays"]
+git-tree-sha1 = "ce7c269e4437fa023d41246b7fc02808af5cb809"
+uuid = "ab1be465-a7f0-4423-9048-0ee774b70ed9"
+version = "0.0.3"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "bc40f042cfcc56230f781d92db71f0e21496dffd"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.1.5"
+
+[[deps.StackViews]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
+uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
+version = "0.1.1"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "cd56bf18ed715e8b09f06ef8c6b781e6cdc49911"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.4.4"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.StructuredArrays]]
+deps = ["ArrayTools"]
+git-tree-sha1 = "003ec02cfda48e9adaebc6be6724de1e7f09fa6b"
+uuid = "2e8cd123-6858-488d-b42d-068777248635"
+version = "0.2.3"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[deps.TensorCore]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
+uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
+version = "0.1.1"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.TiledIteration]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297"
+uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
+version = "0.3.1"
+
+[[deps.TwoDimensional]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "11e927e476d669e3039c26bc629cacfbeb4bbf36"
+uuid = "1907e7ba-7586-4310-a2ba-dd01462aeb50"
+version = "0.2.1"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.ZippedArrays]]
+git-tree-sha1 = "0fe5b07a9b0bd9d56bcf0f013e603b9abd1aa033"
+uuid = "8051824b-fc5e-4603-89c4-0e6ea551cb88"
+version = "0.1.1"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..b9e8016
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,9 @@
+[deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearMapsAA = "599c1a8e-b958-11e9-0d14-b1e6b2ecea07"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SPECTrecon = "ab1be465-a7f0-4423-9048-0ee774b70ed9"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

From 20439ab57e22495d22f90e54da9c8d1b9d9961a3 Mon Sep 17 00:00:00 2001
From: Zongyu Li <54018575+ZongyuLi-umich@users.noreply.github.com>
Date: Thu, 19 May 2022 13:10:11 -0400
Subject: [PATCH 3/5] toml

---
 Manifest.toml | 340 ++++++++++++++++++++++++++++++++++++++++++++++++++
 Project.toml  |  12 +-
 2 files changed, 344 insertions(+), 8 deletions(-)
 create mode 100644 Manifest.toml

diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000..d1f1763
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,340 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.7.1"
+manifest_format = "2.0"
+
+[[deps.AbstractFFTs]]
+deps = ["ChainRulesCore", "LinearAlgebra"]
+git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
+uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+version = "1.1.0"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.3.3"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.BFloat16s]]
+deps = ["LinearAlgebra", "Printf", "Random", "Test"]
+git-tree-sha1 = "a598ecb0d717092b5539dbbe890c98bac842b072"
+uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
+version = "0.2.0"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.CEnum]]
+git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
+uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+version = "0.4.2"
+
+[[deps.CUDA]]
+deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "TimerOutputs"]
+git-tree-sha1 = "19fb33957a5f85efb3cc10e70cf4dd4e30174ac9"
+uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
+version = "3.10.0"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.14.0"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.3"
+
+[[deps.Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.43.0"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[deps.Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.6"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[deps.ExprTools]]
+git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d"
+uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
+version = "0.1.8"
+
+[[deps.FFTW]]
+deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
+git-tree-sha1 = "505876577b5481e50d089c1c68899dfb6faebc62"
+uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+version = "1.4.6"
+
+[[deps.FFTW_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
+uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
+version = "3.3.10+0"
+
+[[deps.GPUArrays]]
+deps = ["Adapt", "LLVM", "LinearAlgebra", "Printf", "Random", "Serialization", "Statistics"]
+git-tree-sha1 = "c783e8883028bf26fb05ed4022c450ef44edd875"
+uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+version = "8.3.2"
+
+[[deps.GPUCompiler]]
+deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "TimerOutputs", "UUIDs"]
+git-tree-sha1 = "d8c5999631e1dc18d767883f621639c838f8e632"
+uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
+version = "0.15.2"
+
+[[deps.IntelOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2018.0.3+2"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "336cc738f03e069ef2cac55a104eb823455dca75"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.4"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.1"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[deps.LLVM]]
+deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"]
+git-tree-sha1 = "c8d47589611803a0f3b4813d9e267cd4e3dbcefb"
+uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
+version = "4.11.1"
+
+[[deps.LLVMExtra_jll]]
+deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg", "TOML"]
+git-tree-sha1 = "771bfe376249626d3ca12bcd58ba243d3f961576"
+uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
+version = "0.0.16+0"
+
+[[deps.LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.15"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.MKL_jll]]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "e595b205efd49508358f7dc670a940c790204629"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2022.0.0+0"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.Random123]]
+deps = ["Random", "RandomNumbers"]
+git-tree-sha1 = "afeacaecf4ed1649555a19cb2cad3c141bbc9474"
+uuid = "74087812-796a-5b5d-8853-05524746bad3"
+version = "1.5.0"
+
+[[deps.RandomNumbers]]
+deps = ["Random", "Requires"]
+git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111"
+uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
+version = "1.5.3"
+
+[[deps.Reexport]]
+git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.2.2"
+
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.0"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "bc40f042cfcc56230f781d92db71f0e21496dffd"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.1.5"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.TimerOutputs]]
+deps = ["ExprTools", "Printf"]
+git-tree-sha1 = "7638550aaea1c9a1e86817a231ef0faa9aca79bd"
+uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+version = "0.5.19"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/Project.toml b/Project.toml
index f521d06..81a18a4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,9 +1,5 @@
-name = "CuSPECTrecon"
-uuid = "b1555097-7ec8-4ced-9fc9-cef29521044b"
-authors = ["Zongyu Li and Jeff Fessler"]
-version = "0.0.1"
-
 [deps]
-
-[compat]
-julia = "1.6"
+AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 700608323cac8c1c90ec078cb623851c693245c0 Mon Sep 17 00:00:00 2001
From: Zongyu Li <54018575+ZongyuLi-umich@users.noreply.github.com>
Date: Tue, 24 May 2022 11:01:43 -0400
Subject: [PATCH 4/5] correction of backproject

---
 src/Cubackproject.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Cubackproject.jl b/src/Cubackproject.jl
index 3e32bcb..bcd9939 100644
--- a/src/Cubackproject.jl
+++ b/src/Cubackproject.jl
@@ -18,7 +18,7 @@ function Cubackproject!(
         # rotate mumap by "-angle"
         Cuimrotate!((@view plan.mumapr[:, :, z]),
                     (@view plan.mumap[:, :, z]),
-                    -plan.viewangle[viewidx],
+                    plan.viewangle[viewidx],
                     plan.planrot)
     end
 

From e6cb351e8606ca21b7799ef0a26098f861cd7aac Mon Sep 17 00:00:00 2001
From: Zongyu Li <54018575+ZongyuLi-umich@users.noreply.github.com>
Date: Tue, 24 May 2022 11:02:12 -0400
Subject: [PATCH 5/5] test backproject

---
 test/Cubackproject.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/test/Cubackproject.jl b/test/Cubackproject.jl
index dc304cc..00a118a 100644
--- a/test/Cubackproject.jl
+++ b/test/Cubackproject.jl
@@ -31,7 +31,7 @@ using Random: seed!
 
     backproject!(image, views, plan)
     Cubackproject!(Cuimage, Cuviews, Cuplan)
-    @test isapprox(Array(Cuimage), image; rtol = 0.6) # this is very *large* difference
+    @test isapprox(Array(Cuimage), image; rtol = 5e-2) 
 end
 
 
@@ -63,6 +63,5 @@ end
 
     Cuproject!(forviews, image, Cuplan)
     Cubackproject!(backimage, views, Cuplan)
-    @test isapprox(dot(forviews, views), dot(backimage, image); rtol = 0.05)
-    # have to put rtol = 0.05 to pass the test
+    @test isapprox(dot(forviews, views), dot(backimage, image); rtol = 1e-2)
 end