diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 35c28c4..bf6a5e1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,7 +23,7 @@ jobs: strategy: fail-fast: false matrix: - version: ['1.6', '1', 'nightly'] + version: ['1', 'nightly'] os: [ubuntu-latest, windows-latest, macOS-latest] steps: - uses: actions/checkout@v3 diff --git a/README.md b/README.md index e7d33e8..dff4566 100644 --- a/README.md +++ b/README.md @@ -24,15 +24,20 @@ The examples include an illustration of how to integrate deep learning into SPECT reconstruction. -For a GPU-focused version, see +Tested with Julia ≥ 1.8. + +## Related packages + +- Pytorch/GPU version: +https://github.com/ZongyuLi-umich/SPECTrecon-pytorch + +- Julia/GPU version: [CuSPECTrecon.jl](https://github.com/JuliaImageRecon/CuSPECTrecon.jl). Designed for use with the [Michigan Image Reconstruction Toolbox (MIRT)](https://github.com/JeffFessler/MIRT.jl) or similar frameworks. -Tested with Julia ≥ 1.6. - [action-img]: https://github.com/JuliaImageRecon/SPECTrecon.jl/workflows/CI/badge.svg [action-url]: https://github.com/JuliaImageRecon/SPECTrecon.jl/actions diff --git a/src/backproject.jl b/src/backproject.jl index f767c22..c02613d 100644 --- a/src/backproject.jl +++ b/src/backproject.jl @@ -13,7 +13,7 @@ function backproject!( viewidx::Int, ) - Threads.@threads for z = 1:plan.imgsize[3] # 1:nz + Threads.@threads for z in 1:plan.imgsize[3] # 1:nz thid = Threads.threadid() # thread id # rotate mumap @@ -26,11 +26,11 @@ function backproject!( end # adjoint of convolving img with psf and applying attenuation map - Threads.@threads for y = 1:plan.imgsize[2] # 1:ny + Threads.@threads for y in 1:plan.imgsize[2] # 1:ny thid = Threads.threadid() # thread id # account for half of the final slice thickness scale3dj!(plan.exp_mumapr[thid], plan.mumapr, y, -0.5) - for j = 1:y + for j in 1:y plus3dj!(plan.exp_mumapr[thid], plan.mumapr, j) end @@ -48,7 +48,7 @@ function backproject!( end # adjoint of rotating image - Threads.@threads for z = 1:plan.imgsize[3] # 1:nz + Threads.@threads for z in 1:plan.imgsize[3] # 1:nz thid = Threads.threadid() imrotate_adj!((@view image[:, :, z]), @@ -74,7 +74,7 @@ function backproject!( viewidx::Int, ) - for z = 1:plan.imgsize[3] # 1:nz + for z in 1:plan.imgsize[3] # 1:nz # thid = Threads.threadid() # thread id # rotate mumap @@ -87,11 +87,11 @@ function backproject!( end # adjoint of convolving img with psf and applying attenuation map - for y = 1:plan.imgsize[2] # 1:ny + for y in 1:plan.imgsize[2] # 1:ny thid = Threads.threadid() # thread id # account for half of the final slice thickness scale3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], y, -0.5) - for j = 1:y + for j in 1:y plus3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], j) end @@ -109,7 +109,7 @@ function backproject!( end # adjoint of rotating image - for z = 1:plan.imgsize[3] # 1:nz + for z in 1:plan.imgsize[3] # 1:nz imrotate_adj!((@view plan.imgr[thid][:, :, z]), (@view plan.imgr[thid][:, :, z]), plan.viewangle[viewidx], @@ -138,13 +138,13 @@ function backproject!( # loop over each view index image .= zero(plan.T) # must be initialized as zero if plan.mode === :fast - [plan.add_img[i] .= zero(plan.T) for i = 1:plan.nthread] + [plan.add_img[i] .= zero(plan.T) for i in 1:plan.nthread] Threads.@threads for (i, viewidx) in collect(enumerate(index)) thid = Threads.threadid() backproject!(plan.add_img[thid], (@view views[:, :, i]), plan, thid, viewidx) end - for i = 1:plan.nthread + for i in 1:plan.nthread broadcast!(+, image, image, plan.add_img[i]) end else diff --git a/src/fft_convolve.jl b/src/fft_convolve.jl index 333ea41..1a60790 100644 --- a/src/fft_convolve.jl +++ b/src/fft_convolve.jl @@ -93,25 +93,25 @@ function fft_conv_adj!( (M, N) = size(img) # adjoint of replicate padding plan.workvecz .= zero(T) - for i = 1:plan.padsize[1] + for i in 1:plan.padsize[1] plus2di!(plan.workvecz, plan.workmat, i) end plus1di!(plan.workmat, plan.workvecz, 1+plan.padsize[1]) plan.workvecz .= zero(T) - for i = (plan.padsize[1]+M+1):size(plan.workmat, 1) + for i in (plan.padsize[1]+M+1):size(plan.workmat, 1) plus2di!(plan.workvecz, plan.workmat, i) end plus1di!(plan.workmat, plan.workvecz, M+plan.padsize[1]) plan.workvecx .= zero(T) - for j = 1:plan.padsize[3] + for j in 1:plan.padsize[3] plus2dj!(plan.workvecx, plan.workmat, j) end plus1dj!(plan.workmat, plan.workvecx, 1+plan.padsize[3]) plan.workvecx .= zero(T) - for j = (plan.padsize[3]+N+1):size(plan.workmat, 2) + for j in (plan.padsize[3]+N+1):size(plan.workmat, 2) plus2dj!(plan.workvecx, plan.workmat, j) end plus1dj!(plan.workmat, plan.workvecx, N+plan.padsize[3]) diff --git a/src/helper.jl b/src/helper.jl index cc1ee38..f39cec4 100644 --- a/src/helper.jl +++ b/src/helper.jl @@ -29,7 +29,7 @@ Base.@propagate_inbounds function padzero!( size(img) .+ (padsize[1] + padsize[2], padsize[3] + padsize[4]) || throw("size") M, N = size(img) output .= zero(T) - for j = padsize[3] + 1:padsize[3] + N, i = padsize[1] + 1:padsize[1] + M + for j in padsize[3] + 1:padsize[3] + N, i in padsize[1] + 1:padsize[1] + M @inbounds output[i, j] = img[i - padsize[1], j - padsize[3]] end return output @@ -50,7 +50,7 @@ Base.@propagate_inbounds function padrepl!( @boundscheck size(output) == size(img) .+ (padsize[1] + padsize[2], padsize[3] + padsize[4]) || throw("size") M, N = size(img) - for j = 1:size(output, 2), i = 1:size(output, 1) + for j in 1:size(output, 2), i in 1:size(output, 1) @inbounds output[i, j] = img[clamp(i - padsize[1], 1, M), clamp(j - padsize[3], 1, N)] end return output @@ -75,7 +75,7 @@ Base.@propagate_inbounds function pad2sizezero!( dims = size(img) pad_dims = ceil.(Int, (padsize .- dims) ./ 2) output .= zero(T) - for j = pad_dims[2]+1:pad_dims[2]+dims[2], i = pad_dims[1]+1:pad_dims[1]+dims[1] + for j in pad_dims[2]+1:pad_dims[2]+dims[2], i in pad_dims[1]+1:pad_dims[1]+dims[1] @inbounds output[i, j] = img[i - pad_dims[1], j - pad_dims[2]] end return output @@ -119,7 +119,7 @@ Base.@propagate_inbounds function fftshift2!( if size(src,2) == 1 @boundscheck iseven(size(src, 1)) || throw("odd $(size(src))") m = size(src,1) ÷ 2 - for i = 1:m + for i in 1:m @inbounds dst[i, 1] = src[m+i, 1] @inbounds dst[i+m, 1] = src[i, 1] end @@ -128,16 +128,16 @@ Base.@propagate_inbounds function fftshift2!( @boundscheck (iseven(size(src, 1)) && iseven(size(src, 2))) || throw("odd") m, n = div.(size(src), 2) - for j = 1:n, i = 1:m + for j in 1:n, i in 1:m @inbounds dst[i, j] = src[m+i, n+j] end - for j = n+1:2n, i = 1:m + for j in n+1:2n, i in 1:m @inbounds dst[i, j] = src[m+i, j-n] end - for j = 1:n, i = m+1:2m + for j in 1:n, i in m+1:2m @inbounds dst[i, j] = src[i-m, j+n] end - for j = n+1:2n, i = m+1:2m + for j in n+1:2n, i in m+1:2m @inbounds dst[i, j] = src[i-m, j-n] end return dst diff --git a/src/ml-os-em.jl b/src/ml-os-em.jl index 1f775a8..1591797 100644 --- a/src/ml-os-em.jl +++ b/src/ml-os-em.jl @@ -36,7 +36,7 @@ function mlem!( back = similar(x0) copyto!(out, x0) time0 = time() - for iter = 1:niter + for iter in 1:niter chat && (@show iter, extrema(out), time() - time0) mul!(ybar, A, out) @. yratio = ynoisy / (ybar + background) # coalesce broadcast! @@ -70,7 +70,7 @@ function Ablock(plan::SPECTplan, nblocks::Int) nview = plan.nview rem(nview, nblocks) == 0 || throw("nview must be divisible by nblocks!") Ab = Vector{LinearMapAO}(undef, nblocks) - for nb = 1:nblocks + for nb in 1:nblocks viewidx = nb:nblocks:nview forw!(y,x) = project!(y, x, plan; index = viewidx) back!(x,y) = backproject!(x, y, plan; index = viewidx) @@ -106,7 +106,7 @@ function osem!( nx, nz, nview = size(ynoisy) nblocks = length(Ab) asum = Vector{Array{eltype(ynoisy), 3}}(undef, nblocks) - for nb = 1:nblocks + for nb in 1:nblocks asum[nb] = Ab[nb]' * ones(eltype(ynoisy), nx, nz, nview÷nblocks) (asum[nb])[asum[nb] .== 0] .= Inf # avoid divide by zero end @@ -115,8 +115,8 @@ function osem!( back = similar(x0) copyto!(out, x0) time0 = time() - for iter = 1:niter - for nb = 1:nblocks + for iter in 1:niter + for nb in 1:nblocks chat && (@show iter, nb, extrema(out), time() - time0) mul!(ybar, Ab[nb], out) @. yratio = (@view ynoisy[:,:,nb:nblocks:nview]) / diff --git a/src/plan-psf.jl b/src/plan-psf.jl index 93da53f..6b3460f 100644 --- a/src/plan-psf.jl +++ b/src/plan-psf.jl @@ -104,7 +104,7 @@ function plan_psf( ; nthread::Int = Threads.nthreads(), T::DataType = Float32, ) - return [PlanPSF( ; nx, nz, px, pz, T) for id = 1:nthread] + return [PlanPSF( ; nx, nz, px, pz, T) for id in 1:nthread] end diff --git a/src/plan-rotate.jl b/src/plan-rotate.jl index 55162b0..f791858 100644 --- a/src/plan-rotate.jl +++ b/src/plan-rotate.jl @@ -99,7 +99,7 @@ function plan_rotate( method::Symbol = :two, nthread::Int = Threads.nthreads(), ) - return [PlanRotate(nx; T, method) for id = 1:nthread] + return [PlanRotate(nx; T, method) for id in 1:nthread] end diff --git a/src/project.jl b/src/project.jl index 35a2844..0a0bc31 100644 --- a/src/project.jl +++ b/src/project.jl @@ -15,7 +15,7 @@ function project!( ) # rotate image and mumap using multiple processors - Threads.@threads for z = 1:plan.imgsize[3] # 1:nz + Threads.@threads for z in 1:plan.imgsize[3] # 1:nz thid = Threads.threadid() # thread id # rotate image in plan.imgr imrotate!( @@ -33,12 +33,12 @@ function project!( ) end - Threads.@threads for y = 1:plan.imgsize[2] # 1:ny + Threads.@threads for y in 1:plan.imgsize[2] # 1:ny thid = Threads.threadid() # thread id # account for half of the final slice thickness scale3dj!(plan.exp_mumapr[thid], plan.mumapr, y, -0.5) - for j = 1:y + for j in 1:y plus3dj!(plan.exp_mumapr[thid], plan.mumapr, j) end @@ -82,7 +82,7 @@ function project!( ) # rotate image and mumap using multiple processors - for z = 1:plan.imgsize[3] # 1:nz + for z in 1:plan.imgsize[3] # 1:nz # rotate image in plan.imgr imrotate!( (@view plan.imgr[thid][:, :, z]), @@ -99,11 +99,11 @@ function project!( ) end - for y = 1:plan.imgsize[2] # 1:ny + for y in 1:plan.imgsize[2] # 1:ny # account for half of the final slice thickness scale3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], y, -0.5) - for j = 1:y + for j in 1:y plus3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], j) end diff --git a/src/psf-gauss.jl b/src/psf-gauss.jl index c2198d5..0dbff28 100644 --- a/src/psf-gauss.jl +++ b/src/psf-gauss.jl @@ -14,7 +14,7 @@ having specified full-width half-maximum (FHWM) values. - 'pz::Int = px' (should be odd) - 'fwhm_start::Real = 1' - 'fwhm_end::Real = 4' -- 'fwhm::AbstractVector{<:Real} = LinRange(fwhm_start, fwhm_end, ny)' +- 'fwhm::AbstractVector{<:Real} = range(fwhm_start, fwhm_end, ny)' - 'fwhm_x::AbstractVector{<:Real} = fwhm, - 'fwhm_z::AbstractVector{<:Real} = fwhm_x' - 'T::DataType == Float32' @@ -27,7 +27,7 @@ function psf_gauss( ; pz::Int = px, fwhm_start::Real = 1, fwhm_end::Real = 4, - fwhm::AbstractVector{<:Real} = LinRange(fwhm_start, fwhm_end, ny), + fwhm::AbstractVector{<:Real} = range(fwhm_start, fwhm_end, ny), fwhm_x::AbstractVector{<:Real} = fwhm, fwhm_z::AbstractVector{<:Real} = fwhm_x, T::DataType = Float32, diff --git a/src/rotatez.jl b/src/rotatez.jl index cd5b8c3..5abd9dc 100644 --- a/src/rotatez.jl +++ b/src/rotatez.jl @@ -26,7 +26,7 @@ function linearinterp!( dec = ceil(Int, s) - s ncoeff = length(A.C) ncol = len - for i = 1:ncoeff + for i in 1:ncoeff if isodd(i) A.C[i] = dec else @@ -36,12 +36,12 @@ function linearinterp!( if e <= ncol A.J[end] = ceil(Int, e) - for i = ncoeff-1:-1:1 + for i in ncoeff-1:-1:1 A.J[i] = max(1, A.J[end] - ceil(Int, (ncoeff - i) / 2)) end else A.J[1] = floor(Int, s) - for i = 2:ncoeff + for i in 2:ncoeff A.J[i] = min(ncol, A.J[1] + ceil(Int, (i - 1) / 2)) end end @@ -153,7 +153,7 @@ In place version of `rotl90`, returning rotation of `A` in `B`. function rotl90!(B::AbstractMatrix, A::AbstractMatrix) ind1, ind2 = axes(A) n = first(ind2) + last(ind2) - for i = axes(A, 1), j = ind2 + for i in axes(A, 1), j in ind2 B[n - j, i] = A[i, j] end return B @@ -167,7 +167,7 @@ In place version of `rotr90`, returning rotation of `A` in `B`. function rotr90!(B::AbstractMatrix, A::AbstractMatrix) ind1, ind2 = axes(A) m = first(ind1) + last(ind1) - for i = ind1, j = axes(A, 2) + for i in ind1, j in axes(A, 2) B[j, m - i] = A[i, j] end return B @@ -181,7 +181,7 @@ In place version of `rot180`, returning rotation of `A` in `B`. function rot180!(B::AbstractMatrix, A::AbstractMatrix) ind1, ind2 = axes(A,1), axes(A,2) m, n = first(ind1)+last(ind1), first(ind2)+last(ind2) - for j=ind2, i=ind1 + for j in ind2, i in ind1 B[m-i,n-j] = A[i, j] end return B @@ -512,7 +512,7 @@ function _imrotate!( size(image3,1) == plans[1].nx || throw("nx") length(plans) == Threads.nthreads() || throw("#threads") - Threads.@threads for z = 1:size(image3,3) # 1:nz + Threads.@threads for z in 1:size(image3,3) # 1:nz id = Threads.threadid() # thread id fun( (@view output[:, :, z]), diff --git a/src/spectplan.jl b/src/spectplan.jl index f5483fa..6905626 100644 --- a/src/spectplan.jl +++ b/src/spectplan.jl @@ -91,11 +91,11 @@ struct SPECTplan{T} if mode === :fast # imgr stores 3D image in different view angles - imgr = [Array{T, 3}(undef, nx, ny, nz) for id = 1:nthread] + imgr = [Array{T, 3}(undef, nx, ny, nz) for id in 1:nthread] # add_img stores 3d image for backprojection - add_img = [Array{T, 3}(undef, nx, ny, nz) for id = 1:nthread] + add_img = [Array{T, 3}(undef, nx, ny, nz) for id in 1:nthread] # mumapr stores 3D mumap in different view angles - mumapr = [Array{T, 3}(undef, nx, ny, nz) for id = 1:nthread] + mumapr = [Array{T, 3}(undef, nx, ny, nz) for id in 1:nthread] else # imgr stores 3D image in different view angles imgr = Array{T, 3}(undef, nx, ny, nz) @@ -105,7 +105,7 @@ struct SPECTplan{T} mumapr = Array{T, 3}(undef, nx, ny, nz) end - exp_mumapr = [Matrix{T}(undef, nx, nz) for id = 1:nthread] + exp_mumapr = [Matrix{T}(undef, nx, nz) for id in 1:nthread] planrot = plan_rotate(nx; T, method = interpmeth)