Use in, abandon 1.6, add gpu urls (#44)

JeffFessler · web-flow · commit 25d9fcf5174f · 2023-02-13T15:45:41.000-05:00
* Use `in`

* Abandon 1.6

* urls

* fix label
diff --git 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
@@ -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.
-
 <!-- URLs -->
 [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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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

Original file line number	Diff line number	Diff line change
`@@ -104,7 +104,7 @@ function plan_psf( ;`
`104`	`104`	`nthread::Int = Threads.nthreads(),`
`105`	`105`	`T::DataType = Float32,`
`106`	`106`	`)`
`107`		`- return [PlanPSF( ; nx, nz, px, pz, T) for id = 1:nthread]`
	`107`	`+ return [PlanPSF( ; nx, nz, px, pz, T) for id in 1:nthread]`
`108`	`108`	`end`
`109`	`109`
`110`	`110`
Original file line number	Diff line number	Diff line change
`@@ -99,7 +99,7 @@ function plan_rotate(`
`99`	`99`	`method::Symbol = :two,`
`100`	`100`	`nthread::Int = Threads.nthreads(),`
`101`	`101`	`)`
`102`		`- return [PlanRotate(nx; T, method) for id = 1:nthread]`
	`102`	`+ return [PlanRotate(nx; T, method) for id in 1:nthread]`
`103`	`103`	`end`
`104`	`104`
`105`	`105`