From 6fca4ce7a5a85888be98a136fd74a6b481a645db Mon Sep 17 00:00:00 2001 From: kamesy Date: Fri, 26 Aug 2022 15:53:25 -0700 Subject: [PATCH] Improve error messages, formatting --- src/QSM.jl | 3 +- src/bgremove/ismv.jl | 14 +- src/bgremove/lbv.jl | 2 +- src/bgremove/pdf.jl | 21 +- src/bgremove/sharp.jl | 58 ++--- src/inversion/direct.jl | 29 +-- src/inversion/nltv.jl | 31 ++- src/inversion/rts.jl | 32 ++- src/inversion/tv.jl | 31 ++- src/unwrap/laplacian.jl | 83 ++++---- src/utils/fd.jl | 48 ++--- src/utils/generate.jl | 84 -------- src/utils/kernels.jl | 46 ++-- src/utils/lsmr.jl | 6 +- src/utils/models.jl | 199 ------------------ src/utils/multi_echo.jl | 37 ++-- src/utils/poisson_solver/mgpcg.jl | 16 +- src/utils/poisson_solver/multigrid/poisson.jl | 22 +- .../poisson_solver/multigrid/smoothers.jl | 4 +- .../poisson_solver/multigrid/transfer.jl | 8 +- src/utils/poisson_solver/poisson_solver.jl | 15 +- src/utils/r2star.jl | 12 +- src/utils/utils.jl | 38 ++-- 23 files changed, 267 insertions(+), 572 deletions(-) delete mode 100644 src/utils/generate.jl delete mode 100644 src/utils/models.jl diff --git a/src/QSM.jl b/src/QSM.jl index 7695178..1dde493 100644 --- a/src/QSM.jl +++ b/src/QSM.jl @@ -80,8 +80,7 @@ const FFTW_NTHREADS = Ref{Int}(known(num_cores())) end function fftw_set_threading(lib::Symbol = :FFTW) - lib ∈ (:FFTW, :Polyester, :Threads) || - throw(ArgumentError("lib must be one of :FFTW, :Polyester or :Threads, got :$(lib)")) + checkopts(lib, (:FFTW, :Polyester, :Threads), :lib) if lib ∈ (:Polyester, :Threads) && nthreads() < 2 @warn "Cannot use $lib with FFTW. Defaulting to FFTW multi-threading" Threads.nthreads() diff --git a/src/bgremove/ismv.jl b/src/bgremove/ismv.jl index 48a8ec5..7bc4f8d 100644 --- a/src/bgremove/ismv.jl +++ b/src/bgremove/ismv.jl @@ -88,24 +88,24 @@ function _ismv!( S = _smv_kernel!(S, F̂, s, vsz, r, P) # constants - _δ = one(eltype(s)) - sqrt(eps(eltype(s))) + δ = one(eltype(s)) - sqrt(eps(eltype(s))) # erode mask s = tcopyto!(s, m) # in-place type conversion, reuse smv var m0 = tcopyto!(m0, m) F̂ = mul!(F̂, P, s) - @inbounds @batch for I in eachindex(F̂) + @batch for I in eachindex(F̂) F̂[I] *= S[I] end s = mul!(s, iP, F̂) - @inbounds @batch for I in eachindex(m) - m[I] = s[I] > _δ + @batch for I in eachindex(m) + m[I] = s[I] > δ end # iSMV - @inbounds for t in axes(f, 4) + for t in axes(f, 4) if verbose && size(f, 4) > 1 @printf("Echo: %d/%d\n", t, size(f, 4)) end @@ -135,7 +135,7 @@ end function __ismv!(f::AbstractArray{T}, f0, S, bc, m, iP, F̂, P, tol, maxit, verbose) where {T} - @inbounds @batch for I in eachindex(f0) + @batch for I in eachindex(f0) f0[I] = m[I]*f[I] end @@ -146,7 +146,7 @@ function __ismv!(f::AbstractArray{T}, f0, S, bc, m, iP, F̂, P, tol, maxit, verb @printf("iter%6s\tresnorm\n", "") end - @inbounds for i in 1:maxit + for i in 1:maxit if nr ≤ ϵ break end diff --git a/src/bgremove/lbv.jl b/src/bgremove/lbv.jl index bbde24e..554163f 100644 --- a/src/bgremove/lbv.jl +++ b/src/bgremove/lbv.jl @@ -118,7 +118,7 @@ function _lbv!( end # set boundaries - @inbounds for t in axes(f, 4) + for t in axes(f, 4) fct = @view(f[Rc,t]) @batch for I in eachindex(fct, mc) fct[I] *= mc[I] diff --git a/src/bgremove/pdf.jl b/src/bgremove/pdf.jl index cf2a03d..f97e8db 100644 --- a/src/bgremove/pdf.jl +++ b/src/bgremove/pdf.jl @@ -85,8 +85,7 @@ function _pdf!( checkshape(W, f, (:W, :f)) end - Dkernel ∈ (:k, :kspace, :i, :ispace) || - throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)")) + checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel) # pad to fast fft size fp = padfastfft(@view(f[:,:,:,1]), pad, rfft=true) @@ -126,7 +125,7 @@ function _pdf!( D = _dipole_kernel!(D, F̂, b, sz0, vsz, bdir, P, Dkernel, :rfft) # background mask - @inbounds @batch for I in eachindex(m) + @batch for I in eachindex(m) M̃[I] = !m[I] end @@ -138,7 +137,7 @@ function _pdf!( elseif ndims(W) == 3 # same weights for all echoes b = padarray!(b, W) - @inbounds @batch for I in eachindex(MW) + @batch for I in eachindex(MW) MW[I] = m[I] * b[I] end @@ -147,7 +146,7 @@ function _pdf!( # VOID end - @inbounds for t in axes(f, 4) + for t in axes(f, 4) if verbose && size(f, 4) > 1 @printf("Echo: %d/%d\n", t, size(f, 4)) end @@ -205,17 +204,17 @@ function _A_pdf!(Av, v, W, iP, D, F̂, P, M̃) v = reshape(v, size(W)) x = reshape(Av, size(W)) - @inbounds @batch for I in eachindex(x) + @batch for I in eachindex(x) x[I] = M̃[I] * v[I] end F̂ = mul!(F̂, P, x) - @inbounds @batch for I in eachindex(F̂) + @batch for I in eachindex(F̂) F̂[I] *= D[I] end x = mul!(x, iP, F̂) - @inbounds @batch for I in eachindex(x) + @batch for I in eachindex(x) x[I] *= W[I] end @@ -227,17 +226,17 @@ function _At_pdf!(Atu, u, M̃, iP, D, F̂, P, W) u = reshape(u, size(W)) x = reshape(Atu, size(W)) - @inbounds @batch for I in eachindex(x) + @batch for I in eachindex(x) x[I] = W[I] * u[I] end F̂ = mul!(F̂, P, x) - @inbounds @batch for I in eachindex(F̂) + @batch for I in eachindex(F̂) F̂[I] *= D[I] # conj(D), D is real end x = mul!(x, iP, F̂) - @inbounds @batch for I in eachindex(x) + @batch for I in eachindex(x) x[I] *= M̃[I] end diff --git a/src/bgremove/sharp.jl b/src/bgremove/sharp.jl index 95e7ff9..be96890 100644 --- a/src/bgremove/sharp.jl +++ b/src/bgremove/sharp.jl @@ -80,27 +80,27 @@ function _sharp!( S = _smv_kernel!(S, F̂, s, vsz, r, P) # constants - _thr = convert(eltype(S), thr) - _one = one(eltype(S)) - _zero = zero(eltype(F̂)) - _δ = one(eltype(s)) - sqrt(eps(eltype(s))) + thrT = convert(eltype(S), thr) + oneT = one(eltype(S)) + zeroT = zero(eltype(F̂)) + δ = one(eltype(s)) - sqrt(eps(eltype(s))) # erode mask s = tcopyto!(s, m) # in-place type conversion, reuse smv var F̂ = mul!(F̂, P, s) - @inbounds @batch for I in eachindex(F̂) + @batch for I in eachindex(F̂) F̂[I] *= S[I] - S[I] = _one - S[I] + S[I] = oneT - S[I] end s = mul!(s, iP, F̂) - @inbounds @batch for I in eachindex(m) - m[I] = s[I] > _δ + @batch for I in eachindex(m) + m[I] = s[I] > δ end # SHARP - @inbounds for t in axes(f, 4) + for t in axes(f, 4) if t > 1 fp = padarray!(fp, @view(f[Rc,t])) end @@ -117,7 +117,7 @@ function _sharp!( F̂ = mul!(F̂, P, fp) @batch for I in eachindex(F̂) - F̂[I] = ifelse(abs(S[I]) < _thr, _zero, F̂[I]*inv(S[I])) + F̂[I] = ifelse(abs(S[I]) < thrT, zeroT, F̂[I]*inv(S[I])) end fp = mul!(fp, iP, F̂) @@ -226,34 +226,34 @@ function _sharp!( fp = tfill!(fp, 0) # constants - _thr = convert(eltype(S), thr) - _one = one(eltype(S)) - _zero = zero(eltype(F̂)) - _δ = one(eltype(s)) - sqrt(eps(eltype(s))) + thrT = convert(eltype(S), thr) + oneT = one(eltype(S)) + zeroT = zero(eltype(F̂)) + δ = one(eltype(s)) - sqrt(eps(eltype(s))) # fft of original mask s = tcopyto!(s, mr) # in-place type conversion M̂ = mul!(M̂, P, s) - @inbounds for (i, r) in enumerate(rs) + for (i, r) in enumerate(rs) # get smv kernel S = _smv_kernel!(S, F̂, s, vsz, r, P) # erode mask @batch for I in eachindex(F̂) F̂[I] = S[I] * M̂[I] - S[I] = _one - S[I] + S[I] = oneT - S[I] end s = mul!(s, iP, F̂) @batch for I in eachindex(mr) - mr[I] = s[I] > _δ + mr[I] = s[I] > δ end # high-pass filter first (largest) kernel for deconvolution if i == 1 @batch for I in eachindex(iS) - iS[I] = ifelse(abs(S[I]) < _thr, _zero, inv(S[I])) + iS[I] = ifelse(abs(S[I]) < thrT, zeroT, inv(S[I])) end end @@ -301,7 +301,7 @@ function _sharp!( checkshape(fl, f, (:fl, :f)) checkshape(smask, mask, (:smask, :mask)) checkshape(axes(mask), axes(f)[1:3], (:mask, :f)) - length(r) > 0 || throw(ArgumentError("r must not be empty")) + !isempty(r) || throw(ArgumentError("r must not be empty")) rs = sort!(unique(collect(r)), rev=true) @@ -338,34 +338,34 @@ function _sharp!( fp = tfill!(fp, 0) # constants - _thr = convert(eltype(S), thr) - _one = one(eltype(S)) - _zero = zero(eltype(F̂)) - _δ = one(eltype(s)) - sqrt(eps(eltype(s))) + thrT = convert(eltype(S), thr) + oneT = one(eltype(S)) + zeroT = zero(eltype(F̂)) + δ = one(eltype(s)) - sqrt(eps(eltype(s))) # fft of original mask s = tcopyto!(s, mr) # in-place type conversion M̂ = mul!(M̂, P, s) - @inbounds for (i, r) in enumerate(rs) + for (i, r) in enumerate(rs) # get smv kernel S = _smv_kernel!(S, F̂, s, vsz, r, P) # erode mask @batch for I in eachindex(F̂) F̂[I] = S[I] * M̂[I] - S[I] = _one - S[I] + S[I] = oneT - S[I] end s = mul!(s, iP, F̂) @batch for I in eachindex(mr) - mr[I] = s[I] > _δ + mr[I] = s[I] > δ end # high-pass filter first (largest) kernel for deconvolution if i == 1 @batch for I in eachindex(iS) - iS[I] = ifelse(abs(S[I]) < _thr, _zero, inv(S[I])) + iS[I] = ifelse(abs(S[I]) < thrT, zeroT, inv(S[I])) end end @@ -391,7 +391,7 @@ function _sharp!( # deconvolution + high-pass filter F̂p = mul!(F̂p, P4, fp) - @inbounds for t in axes(F̂p, 4) + for t in axes(F̂p, 4) F̂t = @view(F̂p[:,:,:,t]) @batch for I in eachindex(F̂t) F̂t[I] *= iS[I] @@ -399,7 +399,7 @@ function _sharp!( end fp = mul!(fp, iP4, F̂p) - @inbounds for t in axes(fp, 4) + for t in axes(fp, 4) ft = @view(fp[:,:,:,t]) @batch for I in eachindex(ft) ft[I] *= m[I] diff --git a/src/inversion/direct.jl b/src/inversion/direct.jl index d008bf1..aac5ccc 100644 --- a/src/inversion/direct.jl +++ b/src/inversion/direct.jl @@ -175,16 +175,9 @@ function _kdiv!( checkshape(x, f, (:x, :f)) checkshape(axes(mask), axes(f)[1:3]) - Dkernel ∈ (:k, :kspace, :i, :ispace) || - throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)")) - - method ∈ (:tkd, :tsvd, :tikh) || - throw(ArgumentError("method must be one of :tkd, :tsvd, got :$(method)")) - - if method == :tikh - reg ∈ (:identity, :gradient, :laplacian) || - throw(ArgumentError("reg must be one of :identity, :gradient, :laplacian, got :$(reg)")) - end + checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel) + checkopts(method, (:tkd, :tsvd, :tikh), :method) + method == :tikh && checkopts(reg, (:identity, :gradient, :laplacian), :reg) # pad to fast fft size fp = padfastfft(@view(f[:,:,:,1]), pad, rfft=true) @@ -208,7 +201,7 @@ function _kdiv!( # inverse k-space kernel iD = _kdiv_ikernel!(D, F̂, fp, vsz, P, lambda, reg, method) - @inbounds for t in axes(f, 4) + for t in axes(f, 4) fp = padarray!(fp, @view(f[:,:,:,t])) F̂ = mul!(F̂, P, fp) @@ -243,22 +236,22 @@ function _kdiv_ikernel!( δ = convert(T, lambda) iδ = inv(δ) - @inbounds @batch for I in eachindex(D) + @batch for I in eachindex(D) D[I] = ifelse(abs(D[I]) ≤ δ, copysign(iδ, D[I]), inv(D[I])) end elseif method == :tsvd δ = convert(T, lambda) - _zero = zero(T) + zeroT = zero(T) - @inbounds @batch for I in eachindex(D) - D[I] = ifelse(abs(D[I]) ≤ δ, _zero, inv(D[I])) + @batch for I in eachindex(D) + D[I] = ifelse(abs(D[I]) ≤ δ, zeroT, inv(D[I])) end elseif method == :tikh λ = convert(T, lambda) - _zero = zero(T) + zeroT = zero(T) Γ = similar(D) @@ -273,8 +266,8 @@ function _kdiv_ikernel!( Γ = tmap!(abs2, Γ) end - @inbounds @batch for I in eachindex(D) - D[I] = ifelse(iszero(D[I]) && iszero(Γ[I]), _zero, inv(D[I]*D[I] + λ*Γ[I])*D[I]) + @batch for I in eachindex(D) + D[I] = ifelse(iszero(D[I]) && iszero(Γ[I]), zeroT, inv(D[I]*D[I] + λ*Γ[I])*D[I]) end end diff --git a/src/inversion/nltv.jl b/src/inversion/nltv.jl index d192d47..0120fba 100644 --- a/src/inversion/nltv.jl +++ b/src/inversion/nltv.jl @@ -121,12 +121,11 @@ function _nltv!( end end - Dkernel ∈ (:k, :kspace, :i, :ispace) || - throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)")) + checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel) # convert scalars - _zero = zero(T) - _one = one(T) + zeroT = zero(T) + oneT = one(T) ρ = convert(T, rho) μ = convert(T, mu) @@ -201,7 +200,7 @@ function _nltv!( D = _dipole_kernel!(D, F̂, xp, sz0, vsz, bdir, P, Dkernel, :rfft) L = _laplace_kernel!(L, F̂, xp, vsz, P, negative=true) - @inbounds for t in axes(f, 4) + for t in axes(f, 4) if verbose && size(f, 4) > 1 @printf("Echo: %d/%d\n", t, size(f, 4)) end @@ -211,7 +210,7 @@ function _nltv!( # iA = (μ*D^H*D - ρ*Δ)^-1 @batch for I in eachindex(iA) a = μ*conj(D[I])*D[I] + ρ*L[I] - iA[I] = ifelse(iszero(a), _zero, inv(a)) + iA[I] = ifelse(iszero(a), zeroT, inv(a)) end F̂ = tfill!(F̂, 0) @@ -234,9 +233,9 @@ function _nltv!( end end - ux = tfill!(ux, _zero) - uy = tfill!(uy, _zero) - uz = tfill!(uz, _zero) + ux = tfill!(ux, zeroT) + uy = tfill!(uy, zeroT) + uz = tfill!(uz, zeroT) dx, dy, dz = _gradfp!(dx, dy, dz, fp, vsz) if verbose @@ -297,7 +296,7 @@ function _nltv!( @batch for I in eachindex(ux) λiρ = λiρWx[I] u = ux[I] + dx[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) ux[I] = u - z dx[I] = z - u + z end @@ -305,7 +304,7 @@ function _nltv!( @batch for I in eachindex(uy) λiρ = λiρWy[I] u = uy[I] + dy[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uy[I] = u - z dy[I] = z - u + z end @@ -313,7 +312,7 @@ function _nltv!( @batch for I in eachindex(uz) λiρ = λiρWz[I] u = uz[I] + dz[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uz[I] = u - z dz[I] = z - u + z end @@ -321,21 +320,21 @@ function _nltv!( else @batch for I in eachindex(ux) u = ux[I] + dx[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) ux[I] = u - z dx[I] = z - u + z end @batch for I in eachindex(uy) u = uy[I] + dy[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uy[I] = u - z dy[I] = z - u + z end @batch for I in eachindex(uz) u = uz[I] + dz[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uz[I] = u - z dz[I] = z - u + z end @@ -368,7 +367,7 @@ function _nltv!( φ = fp[I] s, c = sincos_fast(z - φ) - δ = (muladd(w, s, z) - u) / muladd(w, c, _one) + δ = (muladd(w, s, z) - u) / muladd(w, c, oneT) ze[I] = z - δ threadlocal[1] = muladd(δ, δ, threadlocal[1]) diff --git a/src/inversion/rts.jl b/src/inversion/rts.jl index ffc39a7..7bb8eb8 100644 --- a/src/inversion/rts.jl +++ b/src/inversion/rts.jl @@ -91,12 +91,10 @@ function _rts!( checkshape(x, f, (:x, :f)) checkshape(axes(mask), axes(f)[1:3]) - - Dkernel ∈ (:k, :kspace, :i, :ispace) || - throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)")) + checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel) # convert scalars - _zero = zero(T) + zeroT = zero(T) ρ = convert(T, rho) δ = convert(T, delta) @@ -159,11 +157,11 @@ function _rts!( L = _laplace_kernel!(L, F̂, xp, vsz, P, negative=true) # mask of well-conditioned frequencies - @inbounds @batch for I in eachindex(M) - M[I] = ifelse(abs(D[I]) > δ, μ, _zero) + @batch for I in eachindex(M) + M[I] = ifelse(abs(D[I]) > δ, μ, zeroT) end - @inbounds for t in axes(f, 4) + for t in axes(f, 4) if verbose && size(f, 4) > 1 @printf("Echo: %d/%d\n", t, size(f, 4)) end @@ -174,7 +172,7 @@ function _rts!( # Step 1: Well-conditioned ###################################################################### B̂ = mul!(B̂, P, xp) - lsmr!(WS; lambda=_zero, atol=_zero, btol=_zero, maxit=lstol) + lsmr!(WS; lambda=zeroT, atol=zeroT, btol=zeroT, maxit=lstol) xp = mul!(xp, iP, X̂) @@ -190,8 +188,8 @@ function _rts!( @batch for I in eachindex(F̂) a = muladd(ρ, L[I], M[I]) if iszero(a) - F̂[I] = _zero - iA[I] = _zero + F̂[I] = zeroT + iA[I] = zeroT else ia = inv(a) F̂[I] = ia * M[I] * X̂[I] @@ -200,9 +198,9 @@ function _rts!( end nr = typemax(T) - px = tfill!(px, _zero) - py = tfill!(py, _zero) - pz = tfill!(pz, _zero) + px = tfill!(px, zeroT) + py = tfill!(py, zeroT) + pz = tfill!(pz, zeroT) vx, vy, vz = _gradfp!(vx, vy, vz, xp, vsz) if verbose @@ -258,7 +256,7 @@ function _rts!( @batch for I in eachindex(px) d = dx[I] p = px[I] + d - y = ifelse(abs(p) ≤ iρ, _zero, copysign(abs(p)-iρ, p)) + y = ifelse(abs(p) ≤ iρ, zeroT, copysign(abs(p)-iρ, p)) px[I] = p - y vx[I] = y - p + y dx[I] = d - y @@ -267,7 +265,7 @@ function _rts!( @batch for I in eachindex(py) d = dy[I] p = py[I] + d - y = ifelse(abs(p) ≤ iρ, _zero, copysign(abs(p)-iρ, p)) + y = ifelse(abs(p) ≤ iρ, zeroT, copysign(abs(p)-iρ, p)) py[I] = p - y vy[I] = y - p + y dy[I] = d - y @@ -276,7 +274,7 @@ function _rts!( @batch for I in eachindex(pz) d = dz[I] p = pz[I] + d - y = ifelse(abs(p) ≤ iρ, _zero, copysign(abs(p)-iρ, p)) + y = ifelse(abs(p) ≤ iρ, zeroT, copysign(abs(p)-iρ, p)) pz[I] = p - y vz[I] = y - p + y dz[I] = d - y @@ -336,7 +334,7 @@ end function _A_rts!(Dv, v, D) - @inbounds @batch for I in eachindex(Dv) + @batch for I in eachindex(Dv) Dv[I] = D[I]*v[I] end return Dv diff --git a/src/inversion/tv.jl b/src/inversion/tv.jl index 057e91a..95003ad 100644 --- a/src/inversion/tv.jl +++ b/src/inversion/tv.jl @@ -114,11 +114,10 @@ function _tv!( end end - Dkernel ∈ (:k, :kspace, :i, :ispace) || - throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)")) + checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel) # convert scalars - _zero = zero(T) + zeroT = zero(T) ρ = convert(T, rho) μ = convert(T, mu) @@ -189,7 +188,7 @@ function _tv!( D = _dipole_kernel!(D, F̂, xp, sz0, vsz, bdir, P, Dkernel, :rfft) L = _laplace_kernel!(L, F̂, xp, vsz, P, negative=true) - @inbounds for t in axes(f, 4) + for t in axes(f, 4) if verbose && size(f, 4) > 1 @printf("Echo: %d/%d\n", t, size(f, 4)) end @@ -205,7 +204,7 @@ function _tv!( @batch for I in eachindex(iA) a = μ*conj(D[I])*D[I] + ρ*L[I] - iA[I] = ifelse(iszero(a), _zero, inv(a)) + iA[I] = ifelse(iszero(a), zeroT, inv(a)) end fw = padarray!(fw, @view(W[:,:,:,min(t, end)])) @@ -225,8 +224,8 @@ function _tv!( @batch for I in eachindex(F̂) a = conj(D[I])*D[I] + ρ*L[I] if iszero(a) - F̂[I] = _zero - iA[I] = _zero + F̂[I] = zeroT + iA[I] = zeroT else ia = inv(a) F̂[I] = ia * conj(D[I]) * X̂[I] @@ -246,9 +245,9 @@ function _tv!( end end - ux = tfill!(ux, _zero) - uy = tfill!(uy, _zero) - uz = tfill!(uz, _zero) + ux = tfill!(ux, zeroT) + uy = tfill!(uy, zeroT) + uz = tfill!(uz, zeroT) dx, dy, dz = _gradfp!(dx, dy, dz, xp, vsz) if verbose @@ -309,7 +308,7 @@ function _tv!( @batch for I in eachindex(ux) λiρ = λiρWx[I] u = ux[I] + dx[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) ux[I] = u - z dx[I] = z - u + z end @@ -317,7 +316,7 @@ function _tv!( @batch for I in eachindex(uy) λiρ = λiρWy[I] u = uy[I] + dy[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uy[I] = u - z dy[I] = z - u + z end @@ -325,7 +324,7 @@ function _tv!( @batch for I in eachindex(uz) λiρ = λiρWz[I] u = uz[I] + dz[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uz[I] = u - z dz[I] = z - u + z end @@ -333,21 +332,21 @@ function _tv!( else @batch for I in eachindex(ux) u = ux[I] + dx[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) ux[I] = u - z dx[I] = z - u + z end @batch for I in eachindex(uy) u = uy[I] + dy[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uy[I] = u - z dy[I] = z - u + z end @batch for I in eachindex(uz) u = uz[I] + dz[I] - z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u)) + z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u)) uz[I] = u - z dz[I] = z - u + z end diff --git a/src/unwrap/laplacian.jl b/src/unwrap/laplacian.jl index 123056d..865bf15 100644 --- a/src/unwrap/laplacian.jl +++ b/src/unwrap/laplacian.jl @@ -48,10 +48,9 @@ function unwrap_laplacian( solver::Symbol = :mgpcg ) where {T<:AbstractFloat, N} N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d")) - checkshape(axes(mask), axes(phas)[1:3], (:mask, :phas)) - solver ∈ (:dct, :fft, :mgpcg) || - throw(ArgumentError("solver must be one of :dct, :fft, :mgpcg")) + checkshape(axes(mask), axes(phas)[1:3], (:mask, :phas)) + checkopts(solver, (:dct, :fft, :mgpcg), :solver) d2uphas = wrapped_laplacian(phas, vsz) @@ -78,7 +77,7 @@ function unwrap_laplacian( ) # set boundaries - @inbounds for t in axes(d2uphas, 4) + for t in axes(d2uphas, 4) d2ut = @view(d2uphas[:,:,:,t]) @batch for I in eachindex(d2ut, mask) d2ut[I] *= mask[I] @@ -113,7 +112,7 @@ function wrapped_laplacian!( tsz = padded_tilesize(T, (2, 2, 2), 1) R = vec(collect(TileIterator((2:nx-1, 2:ny-1, 2:nz-1), tsz))) - @inbounds @batch per=thread for (I, J, K) in R + @batch per=thread for (I, J, K) in R for k in K for j in J for i in I @@ -145,20 +144,20 @@ function wrapped_laplacian!( tsz = padded_tilesize(T, (2, 2, 2), 1) R = vec(collect(TileIterator((2:nx-1, 2:ny-1, 2:nz-1), tsz))) - @inbounds for t in axes(u, 4) - _u = @view(u[:,:,:,t]) - _d2u = @view(d2u[:,:,:,t]) + for t in axes(u, 4) + ut = @view(u[:,:,:,t]) + d2ut = @view(d2u[:,:,:,t]) @batch per=thread for (I, J, K) in R for k in K for j in J for i in I - Δ = -dz2 *_wrap(_u[i,j,k] - _u[i,j,k-1]) - Δ = muladd(-dy2, _wrap(_u[i,j,k] - _u[i,j-1,k]), Δ) - Δ = muladd(-dx2, _wrap(_u[i,j,k] - _u[i-1,j,k]), Δ) - Δ = muladd( dx2, _wrap(_u[i+1,j,k] - _u[i,j,k]), Δ) - Δ = muladd( dy2, _wrap(_u[i,j+1,k] - _u[i,j,k]), Δ) - Δ = muladd( dz2, _wrap(_u[i,j,k+1] - _u[i,j,k]), Δ) - _d2u[i,j,k] = Δ + Δ = -dz2 *_wrap(ut[i,j,k] - ut[i,j,k-1]) + Δ = muladd(-dy2, _wrap(ut[i,j,k] - ut[i,j-1,k]), Δ) + Δ = muladd(-dx2, _wrap(ut[i,j,k] - ut[i-1,j,k]), Δ) + Δ = muladd( dx2, _wrap(ut[i+1,j,k] - ut[i,j,k]), Δ) + Δ = muladd( dy2, _wrap(ut[i,j+1,k] - ut[i,j,k]), Δ) + Δ = muladd( dz2, _wrap(ut[i,j,k+1] - ut[i,j,k]), Δ) + d2ut[i,j,k] = Δ end end end @@ -179,6 +178,7 @@ function wrapped_laplacian_boundary_neumann!( dx::NTuple{3, Real} ) where {T<:AbstractFloat, N} N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d")) + checkshape(d2u, u, (:d2u, :u)) zeroT = zero(T) @@ -187,37 +187,37 @@ function wrapped_laplacian_boundary_neumann!( R = edge_indices(axes(u)[1:3]) - @inbounds for t in axes(u, 4) - _u = @view(u[:,:,:,t]) - _d2u = @view(d2u[:,:,:,t]) + for t in axes(u, 4) + ut = @view(u[:,:,:,t]) + d2ut = @view(d2u[:,:,:,t]) @batch per=thread for (i, j, k) in R Δ = zeroT if k > 1 - Δ = muladd(-dz2, _wrap(_u[i,j,k] - _u[i,j,k-1]), Δ) + Δ = muladd(-dz2, _wrap(ut[i,j,k] - ut[i,j,k-1]), Δ) end if j > 1 - Δ = muladd(-dy2, _wrap(_u[i,j,k] - _u[i,j-1,k]), Δ) + Δ = muladd(-dy2, _wrap(ut[i,j,k] - ut[i,j-1,k]), Δ) end if i > 1 - Δ = muladd(-dx2, _wrap(_u[i,j,k] - _u[i-1,j,k]), Δ) + Δ = muladd(-dx2, _wrap(ut[i,j,k] - ut[i-1,j,k]), Δ) end if i < nx - Δ = muladd( dx2, _wrap(_u[i+1,j,k] - _u[i,j,k]), Δ) + Δ = muladd( dx2, _wrap(ut[i+1,j,k] - ut[i,j,k]), Δ) end if j < ny - Δ = muladd( dy2, _wrap(_u[i,j+1,k] - _u[i,j,k]), Δ) + Δ = muladd( dy2, _wrap(ut[i,j+1,k] - ut[i,j,k]), Δ) end if k < nz - Δ = muladd( dz2, _wrap(_u[i,j,k+1] - _u[i,j,k]), Δ) + Δ = muladd( dz2, _wrap(ut[i,j,k+1] - ut[i,j,k]), Δ) end - _d2u[i,j,k] = Δ + d2ut[i,j,k] = Δ end end @@ -231,6 +231,7 @@ function wrapped_laplacian_boundary_periodic!( dx::NTuple{3, Real} ) where {T<:AbstractFloat, N} N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d")) + checkshape(d2u, u, (:d2u, :u)) nx, ny, nz = size(u)[1:3] @@ -238,29 +239,29 @@ function wrapped_laplacian_boundary_periodic!( R = edge_indices(axes(u)[1:3]) - @inbounds for t in axes(u, 4) - _u = @view(u[:,:,:,t]) - _d2u = @view(d2u[:,:,:,t]) + for t in axes(u, 4) + ut = @view(u[:,:,:,t]) + d2ut = @view(d2u[:,:,:,t]) @batch per=thread for (i, j, k) in R - v = k == 1 ? _u[i,j,end] : _u[i,j,k-1] - Δ = -dz2 * _wrap(_u[i,j,k] - v) + v = k == 1 ? ut[i,j,end] : ut[i,j,k-1] + Δ = -dz2 * _wrap(ut[i,j,k] - v) - v = j == 1 ? _u[i,end,k] : _u[i,j-1,k] - Δ = muladd(-dy2, _wrap(_u[i,j,k] - v), Δ) + v = j == 1 ? ut[i,end,k] : ut[i,j-1,k] + Δ = muladd(-dy2, _wrap(ut[i,j,k] - v), Δ) - v = i == 1 ? _u[end,j,k] : _u[i-1,j,k] - Δ = muladd(-dx2, _wrap(_u[i,j,k] - v), Δ) + v = i == 1 ? ut[end,j,k] : ut[i-1,j,k] + Δ = muladd(-dx2, _wrap(ut[i,j,k] - v), Δ) - v = i == nx ? _u[1,j,k] : _u[i+1,j,k] - Δ = muladd(dx2, _wrap(v - _u[i,j,k]), Δ) + v = i == nx ? ut[1,j,k] : ut[i+1,j,k] + Δ = muladd(dx2, _wrap(v - ut[i,j,k]), Δ) - v = j == ny ? _u[i,1,k] : _u[i,j+1,k] - Δ = muladd(dy2, _wrap(v - _u[i,j,k]), Δ) + v = j == ny ? ut[i,1,k] : ut[i,j+1,k] + Δ = muladd(dy2, _wrap(v - ut[i,j,k]), Δ) - v = k == nz ? _u[i,j,1] : _u[i,j,k+1] - Δ = muladd(dz2, _wrap(v - _u[i,j,k]), Δ) + v = k == nz ? ut[i,j,1] : ut[i,j,k+1] + Δ = muladd(dz2, _wrap(v - ut[i,j,k]), Δ) - _d2u[i,j,k] = Δ + d2ut[i,j,k] = Δ end end diff --git a/src/utils/fd.jl b/src/utils/fd.jl index 4ef2649..81aceec 100644 --- a/src/utils/fd.jl +++ b/src/utils/fd.jl @@ -64,7 +64,7 @@ function _gradfp!( hy = convert(Tdy, inv(h[2])) hz = convert(Tdz, inv(h[3])) - @inbounds @batch for k in 1:nz-1 + @batch for k in 1:nz-1 for j in 1:ny-1 for i in 1:nx-1 dx[i,j,k] = hx*(u[i+1,j,k] - u[i,j,k]) @@ -88,7 +88,7 @@ function _gradfp!( dz[nx,ny,k] = hz*(u[nx,ny,k+1] - u[nx,ny,k]) end - @inbounds @batch for j in 1:ny-1 + @batch for j in 1:ny-1 for i in 1:nx-1 dz[i,j,nz] = hz*(u[i,j,1] - u[i,j,nz]) dx[i,j,nz] = hx*(u[i+1,j,nz] - u[i,j,nz]) @@ -100,17 +100,15 @@ function _gradfp!( dy[nx,j,nz] = hy*(u[nx,j+1,nz] - u[nx,j,nz]) end - @inbounds @batch for i in 1:nx-1 + @batch for i in 1:nx-1 dz[i,ny,nz] = hz*(u[i,ny,1] - u[i,ny,nz]) dy[i,ny,nz] = hy*(u[i,1,nz] - u[i,ny,nz]) dx[i,ny,nz] = hx*(u[i+1,ny,nz] - u[i,ny,nz]) end - @inbounds begin - dz[nx,ny,nz] = hz*(u[nx,ny,1] - u[nx,ny,nz]) - dy[nx,ny,nz] = hy*(u[nx,1,nz] - u[nx,ny,nz]) - dx[nx,ny,nz] = hx*(u[1,ny,nz] - u[nx,ny,nz]) - end + dz[nx,ny,nz] = hz*(u[nx,ny,1] - u[nx,ny,nz]) + dy[nx,ny,nz] = hy*(u[nx,1,nz] - u[nx,ny,nz]) + dx[nx,ny,nz] = hx*(u[1,ny,nz] - u[nx,ny,nz]) return dx, dy, dz end @@ -185,14 +183,12 @@ function _gradfp_adj!( hy = convert(Tdy, -inv(h[2])) hz = convert(Tdz, -inv(h[3])) - @inbounds begin - d2u[1,1,1] = - hx*(dx[1,1,1] - dx[nx,1,1]) + - hy*(dy[1,1,1] - dy[1,ny,1]) + - hz*(dz[1,1,1] - dz[1,1,nz]) - end + d2u[1,1,1] = + hx*(dx[1,1,1] - dx[nx,1,1]) + + hy*(dy[1,1,1] - dy[1,ny,1]) + + hz*(dz[1,1,1] - dz[1,1,nz]) - @inbounds @batch for i in 2:nx + @batch for i in 2:nx ux = dx[i,1,1] - dx[i-1,1,1] uy = dy[i,1,1] - dy[i,ny,1] uz = dz[i,1,1] - dz[i,1,nz] @@ -203,7 +199,7 @@ function _gradfp_adj!( d2u[i,1,1] = u end - @inbounds @batch for j in 2:ny + @batch for j in 2:ny uy = dy[1,j,1] - dy[1,j-1,1] ux = dx[1,j,1] - dx[nx,j,1] uz = dz[1,j,1] - dz[1,j,nz] @@ -225,7 +221,7 @@ function _gradfp_adj!( end end - @inbounds @batch for k in 2:nz + @batch for k in 2:nz uz = dz[1,1,k] - dz[1,1,k-1] ux = dx[1,1,k] - dx[nx,1,k] uy = dy[1,1,k] - dy[1,ny,k] @@ -322,17 +318,14 @@ function _lap!( u::AbstractArray{T, 3}, h::NTuple{3, Real}, ) where {T<:AbstractFloat} - idx2 = convert(T, inv(h[1]*h[1])) - idy2 = convert(T, inv(h[2]*h[2])) - idz2 = convert(T, inv(h[3]*h[3])) + idx2, idy2, idz2 = convert.(T, inv.(h.*h)) D = -2*(idx2 + idy2 + idz2) - ax = map(a -> first(a)+1:last(a)-1, axes(u)) - tsz = padded_tilesize(T, (2, 2, 2), 1) + ax = map(a -> first(a)+1:last(a)-1, axes(u)) R = vec(collect(TileIterator(ax, tsz))) - @inbounds @batch for (I, J, K) in R + @batch for (I, J, K) in R for k in K for j in J for i in I @@ -367,17 +360,14 @@ function _lap!( u::AbstractArray{T, 4}, h::NTuple{3, Real}, ) where {T<:AbstractFloat} - idx2 = convert(T, inv(h[1]*h[1])) - idy2 = convert(T, inv(h[2]*h[2])) - idz2 = convert(T, inv(h[3]*h[3])) + idx2, idy2, idz2 = convert.(T, inv.(h.*h)) D = -2*(idx2 + idy2 + idz2) - ax = map(a -> first(a)+1:last(a)-1, axes(u)[1:3]) - tsz = padded_tilesize(T, (2, 2, 2), 1) + ax = map(a -> first(a)+1:last(a)-1, axes(u)[1:3]) R = vec(collect(TileIterator(ax, tsz))) - @inbounds for t in axes(u, 4) + for t in axes(u, 4) ut = @view(u[:,:,:,t]) d2ut = @view(d2u[:,:,:,t]) @batch for (I, J, K) in R diff --git a/src/utils/generate.jl b/src/utils/generate.jl deleted file mode 100644 index aa933fe..0000000 --- a/src/utils/generate.jl +++ /dev/null @@ -1,84 +0,0 @@ -using FFTW -using MAT -using Random -using Printf -using StaticArrays - -include("models.jl") - - -""" -Generate magnetic susceptibility maps and the corresponding fields due to -randomly placed spheres and randomly oriented infinite cylinders subject to a -static B-field inside a 3d rectangular volume. - -Method ------- - 1. Pick array size, voxel size for grid (N/n grid) - 2. for n iterations - i. pick number of tries, ntry, to check for valid config - ii. randomly place sources inside grid until no valid config can be - found ntry times in a row - iii. valid config means no overlapping structures. concentric allowed. -""" -function generate( - sz::NTuple{N, Integer}, - vsz::NTuple{3, Real}; - outpath::Union{Nothing, AbstractString} = nothing, - nsources::Integer = typemax(Int)-1, - density::Real = 0.51, - bdir::NTuple{3, Real} = (0, 0, 1), - χμ::Float64 = 1e-6, - χσ::Float64 = 0.07, - seed::Int = 12345, -) where {N} - @assert nsources > 0 - @assert density > 0 && density < 1 - Random.seed!(seed) - - if !(norm(bdir) ≈ 1) - @warn "B direction vector is not normalized. Normalizing..." - bdir = bdir ./ norm(bdir) - end - - x = zeros(sz) - f = zeros(sz) - m = Array{Bool}(falses(sz)) - - rmax = cld(minimum(sz), 5) - rmax = 7 - Cs = ntuple(i -> rmax+1:sz[i]-(rmax+1)-1, Val(3)) - Cs = ntuple(i -> 10*rmax+1:sz[i]-(10*rmax+1)-1, Val(3)) - - ns = 0 - while true - model = rand((:sphere, :cylinder)) - model = :sphere - - r = rand(3:rmax) - c = rand.(Cs) - χ = χμ + χσ*randn() - - if model == :sphere - x, f, m = sphere!(x, f, m, c, r, χ, bdir, vsz) - ns += 1 - - elseif model == :cylinder - θ = π*rand() - ϕ = 2π*rand() - - x, f, m = cylinder!(x, f, m, c, r, θ, ϕ, χ, bdir, vsz) - ns += 1 - end - - if sum(m)/length(m) > density || ns >= nsources - break - end - end - - if outpath !== nothing - matwrite(outpath, Dict(:x => x, :f => f, :m => m), compress=true) - end - - return x, f, m -end diff --git a/src/utils/kernels.jl b/src/utils/kernels.jl index 6f189f9..cdbe0ea 100644 --- a/src/utils/kernels.jl +++ b/src/utils/kernels.jl @@ -59,16 +59,13 @@ function dipole_kernel( transform::Union{Nothing, Symbol} = nothing, shift::Bool = false ) where {T<:AbstractFloat} - !any(<=(0), sz) || throw(ArgumentError("invalid array size")) - !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0")) - !iszero(norm(bdir)) || throw(DomainError(bdir, "norm of bdir must not be zero")) - !any(<(0), dsz) || throw(DomainError(dsz, "dsz must be ≥ 0")) + all(>(0), sz) || throw(ArgumentError("invalid array size")) + all(>(0), vsz) || throw(DomainError(vsz, "vsz must be > 0")) + all(>(0), dsz) || throw(DomainError(dsz, "dsz must be > 0")) + norm(bdir) > 0 || throw(DomainError(bdir, "bdir must not be zero")) - method ∈ (:k, :kspace, :i, :ispace) || - throw(ArgumentError("method must be one of :k, :kspace, :i, :ispace, got :$(method)")) - - transform === nothing || transform ∈ (:fft, :rfft) || - throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)")) + checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel) + transform !== nothing && checkopts(transform, (:fft, :rfft), :transform) if method == :k || method == :kspace D = Array{T, 3}(undef, transform == :rfft ? (sz[1]>>1 + 1, sz[2], sz[3]) : sz) @@ -178,7 +175,7 @@ function _dipole_kernel!( X, Y, Z = _grid(sz, dx, fft=true, shift=shift) a = 1 / 3 - @inbounds @batch for k in 1:nz + @batch for k in 1:nz for j in 1:ny for i in 1:nx k̂ = SVector{3, Float64}(X[i], Y[j], Z[k]) @@ -201,10 +198,10 @@ function _dipole_kernel!( X, Y, Z = _grid(size(D), dx, shift=shift) a = prod(dx) * inv(4*pi) - @inbounds @batch for k in 1:nz + @batch for k in 1:nz z = Z[k] zout = z < Z1 || z > Z2 - @fastpow for j in 1:ny + for j in 1:ny y = Y[j] zyout = zout || y < Y1 || y > Y2 for i in 1:nx @@ -217,7 +214,7 @@ function _dipole_kernel!( r = SVector{3, Float64}(x, y, z) rz = r⋅B r2 = r⋅r - D[i,j,k] = a * (3*rz^2-r2)/sqrt(r2^5) + @fastpow D[i,j,k] = a * (3*rz^2-r2)/sqrt(r2^5) end end end @@ -277,12 +274,9 @@ function smv_kernel( transform::Union{Nothing, Symbol} = nothing, shift::Bool = false ) where {T<:AbstractFloat} - !any(<=(0), sz) || throw(ArgumentError("invalid array size")) - !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0")) - r >= 0 || throw(DomainError(r, "r must be ≥ 0")) - - transform === nothing || transform ∈ (:fft, :rfft) || - throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)")) + all(>(0), sz) || throw(ArgumentError("invalid array size")) + all(>(0), vsz) || throw(DomainError(vsz, "vsz must be > 0")) + transform !== nothing && checkopts(transform, (:fft, :rfft), :transform) if transform === nothing s = Array{T, 3}(undef, sz) @@ -379,14 +373,14 @@ function _sphere!( X, Y, Z = _grid(sz, dx, shift=shift) - _zero = zero(T) - _one = one(T) + zeroT = zero(T) + oneT = one(T) @inbounds for k in 1:nz for j in 1:ny for i in 1:nx R = SVector{3, Float64}(X[i], Y[j], Z[k]) - s[i,j,k] = ifelse(R⋅R <= r2, _one, _zero) + s[i,j,k] = ifelse(R⋅R <= r2, oneT, zeroT) end end end @@ -454,11 +448,9 @@ function laplace_kernel( transform::Union{Nothing, Symbol} = nothing, shift::Bool = false ) where {T<:AbstractFloat} - !any(<=(0), sz) || throw(ArgumentError("invalid array size")) - !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0")) - - transform === nothing || transform ∈ (:fft, :rfft) || - throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)")) + all(>(0), sz) || throw(ArgumentError("invalid array size")) + all(>(0), vsz) || throw(DomainError(vsz, "vsz must be > 0")) + transform !== nothing && checkopts(transform, (:fft, :rfft), :transform) if transform === nothing Δ = Array{T, 3}(undef, sz) diff --git a/src/utils/lsmr.jl b/src/utils/lsmr.jl index ed60498..d6384b1 100644 --- a/src/utils/lsmr.jl +++ b/src/utils/lsmr.jl @@ -365,7 +365,7 @@ end @inline function _xpby_norm!(X::Array, b::Real, Y::Array{T}) where {T} Tr = real(T) bT = convert(Tr, b) - @inbounds @batch threadlocal=zero(T)::T for I in eachindex(Y) + @batch threadlocal=zero(T)::T for I in eachindex(Y) Y[I] = muladd(bT, Y[I], X[I]) threadlocal = muladd(conj(Y[I]), Y[I], threadlocal) end @@ -381,7 +381,7 @@ end @inline function _scal!(a, X::Array{T}) where {T} aT = convert(real(T), a) - @inbounds @batch for I in eachindex(X) + @batch for I in eachindex(X) X[I] *= aT end return X @@ -417,7 +417,7 @@ end σ = ζ / (ρ * ρbar) ν = -θnew / ρ - @inbounds @batch threadlocal=zero(T)::T for I in eachindex(x) + @batch threadlocal=zero(T)::T for I in eachindex(x) h̄[I] = muladd(δ, h̄[I], h[I]) x[I] = muladd(σ, h̄[I], x[I]) h[I] = muladd(ν, h[I], v[I]) diff --git a/src/utils/models.jl b/src/utils/models.jl deleted file mode 100644 index bf5bdbd..0000000 --- a/src/utils/models.jl +++ /dev/null @@ -1,199 +0,0 @@ -using FastPow -using LinearAlgebra -using Polyester -using StaticArrays - - -""" - sphere!(x, f, m, c, r, χ, bdir, vsz) - -Add a sphere to grid. - -### Arguments -- `x::AbstractArray{T1, 3}`: magnetic susceptibilities -- `f::AbstractArray{T2, 3}`: local field generated by magnetic susceptibilities -- `m::AbstractArray{Bool, 3}`: binary mask of magnetic susceptibilities -- `c::Union{NTuple{3}, SVector{3}}`: center of source to add; c[i] ∈ (1:size(x,i)) -- `r::Real`: radius of source to add -- `χ::Real`: magnetic susceptibility of new source in ppm -- `bdir::NTuple{3, Real}`: unit vector of B-field direction (x, y, z) -- `vsz::NTuple{3, Real}`: voxel size - -### Returns -- `x` -- `f` -- `m` - -### References - [1] Cheng, Y. C. N., Neelavalli, J., & Haacke, E. M. (2009). Limitations of - calculating field distributions and magnetic susceptibilities in MRI using - a Fourier based method. Physics in Medicine & Biology, 54(5), 1169. -""" -function sphere!( - x::AbstractArray{T1, 3}, - f::AbstractArray{T2, 3}, - m::AbstractArray{Bool, 3}, - c::Union{NTuple{3}, SVector{3}}, - r::Real, - χ::Real, - bdir::NTuple{3, Real}, - vsz::NTuple{3, Real}, -) where {T1, T2} - size(x) == size(f) || throw(DimensionMismatch()) - size(f) == size(m) || throw(DimensionMismatch()) - checkbounds(Bool, x, CartesianIndex(c...)) || throw(DimensionMismatch()) - - sz = size(x) - nx, ny, nz = sz - - C = SVector{3, Float64}(c) - B = SVector{3, Float64}(bdir) - ns = SVector{3, Float64}(sz) - dx = SVector{3, Float64}(vsz) - - if !(norm(B) ≈ 1) - @warn "B direction vector is not normalized. Normalizing..." - B = B ./ norm(B) - end - - X, Y, Z = _grid(C, ns, dx, sz) - C = C .+ first.((X, Y, Z)) - - chi = Float64(χ) - r2 = Float64(r^2) - a = Float64(1//3*χ*r^3) - - @batch for k in 1:nz - @inbounds for j in 1:ny - @fastpow for i in 1:nx - R = SVector{3, Float64}(X[i], Y[j], Z[k]) - C - R2 = R⋅R - - if R2 <= r2 # inside sphere, field = 0 - x[i,j,k] += chi - m[i,j,k] = true - else # outside sphere, compute field - Rz = R⋅B - _f = (3*Rz^2-R2)/sqrt(R2^5) - f[i,j,k] = muladd(a, _f, f[i,j,k]) - end - end - end - end - - return x, f, m -end - - -""" - cylinder!(x, f, m, c, r, θ, ϕ, χ, bdir, vsz) - -Add an infinite cylinder to grid. - -### Arguments -- `x::AbstractArray{T1, 3}`: magnetic susceptibilities -- `f::AbstractArray{T2, 3}`: local field generated by magnetic susceptibilities -- `m::AbstractArray{Bool, 3}`: binary mask of magnetic susceptibilities -- `c::Union{NTuple{3}, SVector{3}}`: center of source to add; c[i] ∈ (1:size(x,i)) -- `r::Real`: radius of source to add -- `θ::Real`: polar angle wrt B-field of cylinder -- `ϕ::Real`: azimuthal angle wrt B-field of cylinder -- `χ::Real`: magnetic susceptibility of new source in ppm -- `bdir::NTuple{3, Real}`: unit vector of B-field direction (x, y, z) -- `vsz::NTuple{3, Real}`: voxel size - -### Returns -- `x` -- `f` -- `m` - -### References - [1] Cheng, Y. C. N., Neelavalli, J., & Haacke, E. M. (2009). Limitations of - calculating field distributions and magnetic susceptibilities in MRI using - a Fourier based method. Physics in Medicine & Biology, 54(5), 1169. -""" -function cylinder!( - x::AbstractArray{T1, 3}, - f::AbstractArray{T2, 3}, - m::AbstractArray{Bool, 3}, - c::Union{NTuple{3}, SVector{3}}, - r::Real, - θ::Real, - ϕ::Real, - χ::Real, - bdir::NTuple{3, Real}, - vsz::NTuple{3, Real}, -) where {T1, T2} - size(x) == size(f) || throw(DimensionMismatch()) - size(f) == size(m) || throw(DimensionMismatch()) - checkbounds(Bool, x, CartesianIndex(c...)) || throw(DimensionMismatch()) - - sz = size(x) - nx, ny, nz = sz - - C = SVector{3, Float64}(c) - ns = SVector{3, Float64}(sz) - dx = SVector{3, Float64}(vsz) - - if !(norm(bdir) ≈ 1) - @warn "B direction vector is not normalized. Normalizing..." - bdir = bdir ./ norm(bdir) - end - - X, Y, Z = _grid(C, ns, dx, sz) - - bz = Float64(bdir[3]) - theta = Float64(θ) - phi = Float64(ϕ) - chi = Float64(χ) - ρ = Float64(r) - - B = cos(acos(bz) - theta) - F = 1/6 * chi * (3*B^2 - 1) # constant field inside cylinder - a = 0.5*chi * ρ^2 * sin(acos(B))^2 - - e = SVector{3, Float64}(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)) - - @batch for k in 1:nz - @inbounds for j in 1:ny - @fastpow for i in 1:nx - R = SVector{3, Float64}(X[i], Y[j], Z[k]) - - # mathworld.wolfram.com/Point-LineDistance3-Dimensional.html - if norm(R×(R-e)) <= ρ # inside cylinder, constant field - f[I] += F - x[I] += chi - m[I] = true - else # outside cylinder, compute field - if iszero(theta) - j = SVector{3, Float64}(0, 1, 0) - i = SVector{3, Float64}(1, 0, 0) - else - # j = k×z, z = (0, 0, 1), k = e - j = SVector{3, Float64}(e[2], -e[1], 0) - j = j .* inv(norm(j)) - i = j×e - end - Rx2 = (R⋅i)^2 - Ry2 = (R⋅j)^2 - _f = (Rx2 - Ry2)/(Rx2 + Ry2)^2 - - f[i,j,k] = muladd(a, _f, f[i,j,k]) - end - end - end - end - - return x, f, m -end - - -function _grid( - c::Union{NTuple{N}, SVector{N}}, - sz::Union{NTuple{N}, SVector{N}}, - dx::Union{NTuple{N}, SVector{N}}, - ns::Union{NTuple{N}, SVector{N}} -) where {N} - x = dx .* fld.(sz, 2) - return ntuple(i -> range(-x[i], step=dx[i], length=ns[i]), Val(N)) -end diff --git a/src/utils/multi_echo.jl b/src/utils/multi_echo.jl index da8774e..b587cc5 100644 --- a/src/utils/multi_echo.jl +++ b/src/utils/multi_echo.jl @@ -35,10 +35,11 @@ function fit_echo_linear( NT > 1 || throw(ArgumentError("data must be multi-echo")) p = similar(phas, size(phas)[1:N-1]) - if !phase_offset - return fit_echo_linear!(p, phas, W, TEs) + + return if !phase_offset + fit_echo_linear!(p, phas, W, TEs) else - return fit_echo_linear!(p, similar(p), phas, W, TEs) + fit_echo_linear!(p, similar(p), phas, W, TEs) end end @@ -70,8 +71,8 @@ function fit_echo_linear!( M > 1 || throw(ArgumentError("array must contain echoes in last dimension")) NT > 1 || throw(ArgumentError("data must be multi-echo")) - size(phas, M) == NT || throw(DimensionMismatch()) - length(p) == length(phas) ÷ NT || throw(DimensionMismatch()) + checkshape(axes(p), axes(phas)[1:M-1], (:p, :phas)) + checkshape((NT,), (size(phas, M),), (:TEs, :phas)) checkshape(W, phas, (:W, :phas)) vphas = reshape(phas, :, NT) @@ -81,9 +82,9 @@ function fit_echo_linear!( T = promote_type(Tp, Tphas, TW) tes = convert.(T, TEs) - _zeroTp = zero(Tp) + zeroTp = zero(Tp) - @inbounds @batch for I in eachindex(vp) + @batch for I in eachindex(vp) w = vW[I,1] * vW[I,1] * tes[1] den = w * tes[1] num = w * vphas[I,1] @@ -92,7 +93,7 @@ function fit_echo_linear!( den = muladd(w, tes[t], den) num = muladd(w, vphas[I,t], num) end - vp[I] = iszero(den) ? _zeroTp : num * inv(den) + vp[I] = iszero(den) ? zeroTp : num * inv(den) end return p @@ -130,9 +131,9 @@ function fit_echo_linear!( M > 1 || throw(ArgumentError("array must contain echoes in last dimension")) NT > 1 || throw(ArgumentError("data must be multi-echo")) - size(phas, M) == NT || throw(DimensionMismatch()) - length(p) == length(phas) ÷ NT || throw(DimensionMismatch()) - length(p0) == length(phas) ÷ NT || throw(DimensionMismatch()) + checkshape(axes(p), axes(phas)[1:M-1], (:p, :phas)) + checkshape(axes(p0), axes(phas)[1:M-1], (:p0, :phas)) + checkshape((NT,), (size(phas, M),), (:TEs, :phas)) checkshape(W, phas, (:W, :phas)) vphas = reshape(phas, :, NT) @@ -143,10 +144,10 @@ function fit_echo_linear!( T = promote_type(Tp, Tp0, Tphas, TW) tes = convert.(T, TEs) - _zeroTp = zero(Tp) - _zeroTp0 = zero(Tp0) + zeroTp = zero(Tp) + zeroTp0 = zero(Tp0) - @inbounds @batch for I in eachindex(vp) + @batch for I in eachindex(vp) w = vW[I,1] x = w * tes[1] y = w * vphas[I,1] @@ -157,8 +158,8 @@ function fit_echo_linear!( end if iszero(w) - vp[I] = _zeroTp - vp0[I] = _zeroTp0 + vp[I] = zeroTp + vp0[I] = zeroTp0 continue end @@ -183,8 +184,8 @@ function fit_echo_linear!( vp[I] = num * inv(den) vp0[I] = y - vp[I] * x else - vp[I] = _zeroTp - vp0[I] = _zeroTp0 + vp[I] = zeroTp + vp0[I] = zeroTp0 end end diff --git a/src/utils/poisson_solver/mgpcg.jl b/src/utils/poisson_solver/mgpcg.jl index ad6e59f..1c0b3bb 100644 --- a/src/utils/poisson_solver/mgpcg.jl +++ b/src/utils/poisson_solver/mgpcg.jl @@ -99,13 +99,13 @@ function _mgpcg!( maxit::Integer, verbose::Bool ) where {T<:AbstractFloat} - _zero = zero(T) + zeroT = zero(T) ρ = one(T) A = M.grids[1].A q = M.workspace.x[1] r = M.workspace.b[1] - p = tfill!(p, _zero) + p = tfill!(p, zeroT) # q = A * x # r -= q @@ -127,7 +127,7 @@ function _mgpcg!( for i in 1:maxit # q = M\r - tfill!(q, _zero, A) + tfill!(q, zeroT, A) for _ in 1:ncycles cycle!(M, cycle) end @@ -165,7 +165,7 @@ end function _compute_ρ(q::AbstractArray{T}, r, A::Poisson2) where {T} # ρ = q⋅r - @inbounds @batch threadlocal=zero(T)::T for I in A.RI + @batch threadlocal=zero(T)::T for I in A.RI threadlocal = muladd(q[I], r[I], threadlocal) end return sum(threadlocal::Vector{T}) @@ -174,7 +174,7 @@ end function _compute_p!(p, q, β, A::Poisson2) # p = q + β * p - @inbounds @batch for I in A.RI + @batch for I in A.RI p[I] *= β p[I] += q[I] end @@ -190,7 +190,7 @@ function _compute_q!(q::AbstractArray{T}, A::Poisson2, p) where {T} idy2 = convert(T, A.idx2[2]) idz2 = convert(T, A.idx2[3]) - @inbounds @batch per=thread threadlocal=zero(T)::T for (I, J, K) in A.R8 + @batch per=thread threadlocal=zero(T)::T for (I, J, K) in A.R8 for k in K, j in J, i in I z1 = p[i,j,k-1] y1 = p[i,j-1,k] @@ -224,7 +224,7 @@ function _compute_x_r!(x::AbstractArray{T}, r, p, q, α, A::Poisson2) where {T} # x = x + α * p # r = r - α * q # res = norm(r) - @inbounds @batch threadlocal=zero(T)::T for I in A.RI + @batch threadlocal=zero(T)::T for I in A.RI x[I] = muladd( α, p[I], x[I]) r[I] = muladd(-α, q[I], r[I]) threadlocal = muladd(r[I], r[I], threadlocal) @@ -238,7 +238,7 @@ function _init_r!(r::AbstractArray{T}, q, A, x) where {T} # r -= q # res = norm(r) q = mul!(q, A, x) - @inbounds @batch threadlocal=zero(T)::T for I in eachindex(r) + @batch threadlocal=zero(T)::T for I in eachindex(r) r[I] -= q[I] threadlocal = muladd(r[I], r[I], threadlocal) end diff --git a/src/utils/poisson_solver/multigrid/poisson.jl b/src/utils/poisson_solver/multigrid/poisson.jl index 35113c6..3bb3648 100644 --- a/src/utils/poisson_solver/multigrid/poisson.jl +++ b/src/utils/poisson_solver/multigrid/poisson.jl @@ -68,14 +68,12 @@ function LinearAlgebra.mul!( x::AbstractArray{T, 3}, ) where {T<:AbstractFloat} D = convert(T, A.D) - idx2 = convert(T, A.idx2[1]) - idy2 = convert(T, A.idx2[2]) - idz2 = convert(T, A.idx2[3]) + idx2, idy2, idz2 = convert.(T, A.idx2) # @threads ~10% faster here # Polyester needs A.interior inside or the loop slows down by 100%??? # ie. slow down when m = A.interior ... if m[i,j,k] ... - @inbounds @batch for (I, J, K) in A.R8 + @batch for (I, J, K) in A.R8 for k in K for j in J for i in I @@ -113,14 +111,12 @@ function LinearAlgebra.mul!( x::AbstractArray{T, 4}, ) where {T<:AbstractFloat} D = convert(T, A.D) - idx2 = convert(T, A.idx2[1]) - idy2 = convert(T, A.idx2[2]) - idz2 = convert(T, A.idx2[3]) + idx2, idy2, idz2 = convert.(T, A.idx2) # @threads ~10% faster here # Polyester needs A.interior inside or the loop slows down by 100%??? # ie. slow down when m = A.interior ... if m[i,j,k] ... - @inbounds for t in axes(x, 4) + for t in axes(x, 4) xt = @view(x[:,:,:,t]) d2xt = @view(d2x[:,:,:,t]) @batch for (I, J, K) in A.R8 @@ -171,7 +167,7 @@ function residual!( # @threads ~10% faster here # Polyester needs A.interior inside or the loop slows down by 100%??? # ie. slow down when m = A.interior ... if m[i,j,k] ... - @inbounds @batch for (I, J, K) in A.R8 + @batch for (I, J, K) in A.R8 for k in K for j in J for i in I @@ -216,7 +212,7 @@ function __norm_residual( idy2 = convert(T, -A.idx2[2]) idz2 = convert(T, -A.idx2[3]) - @inbounds @batch threadlocal=zero(T)::T for (I, J, K) in A.R8 + @batch threadlocal=zero(T)::T for (I, J, K) in A.R8 for k in K for j in J for i in I @@ -268,7 +264,7 @@ function boundary_mask!( end # interior - @inbounds @batch minbatch=8 for kc in 2:nzc-1 + @batch minbatch=8 for kc in 2:nzc-1 k = (kc << 1) - 1 for jc in 2:nyc-1 j = (jc << 1) - 1 @@ -297,7 +293,7 @@ function boundary_mask!( end end - @inbounds @batch minbatch=8 for k in 2:nz-1 + @batch minbatch=8 for k in 2:nz-1 for j in 2:ny-1 for i in 2:nx-1 if m[i,j,k] @@ -320,7 +316,7 @@ end function tfill!(x::AbstractArray{T, N}, v, A::Poisson2{<:AbstractFloat, N}) where {T, N} vT = convert(T, v) - @inbounds @batch minbatch=1024 for I in A.RI + @batch minbatch=1024 for I in A.RI x[I] = vT end return x diff --git a/src/utils/poisson_solver/multigrid/smoothers.jl b/src/utils/poisson_solver/multigrid/smoothers.jl index 1479085..6470fb6 100644 --- a/src/utils/poisson_solver/multigrid/smoothers.jl +++ b/src/utils/poisson_solver/multigrid/smoothers.jl @@ -29,7 +29,7 @@ struct Smoothers{N, T} <: AbstractSmoother end function (s::Smoothers{N})(x, A, b) where {N} - @inbounds for i in 1:N + for i in 1:N x = s.s[i](x, A, b) end return x @@ -203,7 +203,7 @@ function rbgs!( idz2 = -iD * A.idx2[3] for s in C - @inbounds @batch minbatch=8 for (I, J, K) in R + @batch minbatch=8 for (I, J, K) in R for k in K ks = Bool(k & 1) for j in J diff --git a/src/utils/poisson_solver/multigrid/transfer.jl b/src/utils/poisson_solver/multigrid/transfer.jl index bf801d0..b3911ac 100644 --- a/src/utils/poisson_solver/multigrid/transfer.jl +++ b/src/utils/poisson_solver/multigrid/transfer.jl @@ -17,7 +17,7 @@ end function restrict!(mc::AbstractArray{Bool, 3}, m::AbstractArray{Bool, 3}) nxc, nyc, nzc = size(mc) - @inbounds @batch minbatch=8 for kc in 2:nzc-1 + @batch minbatch=8 for kc in 2:nzc-1 k = (kc << 1) - 1 for jc in 2:nyc-1 j = (jc << 1) - 1 @@ -51,7 +51,7 @@ function restrict!( a2 = convert(T, 2//64) a3 = convert(T, 1//64) - @inbounds @batch minbatch=8 for (Ic, Jc, Kc) in Ac.R27 + @batch minbatch=8 for (Ic, Jc, Kc) in Ac.R27 for kc in Kc k = (kc << 1) - 1 for jc in Jc @@ -134,7 +134,7 @@ function prolong!( nxc, nyc, nzc = size(xc) m = A.interior - @inbounds @batch for kc in 1:nzc-1 + @batch for kc in 1:nzc-1 k = (kc << 1) - 1 for jc in 1:nyc-1 j = (jc << 1) - 1 @@ -209,7 +209,7 @@ function correct_prolong!( nxc, nyc, nzc = size(xc) m = A.interior - @inbounds @batch for kc in 1:nzc-1 + @batch for kc in 1:nzc-1 k = (kc << 1) - 1 for jc in 1:nyc-1 j = (jc << 1) - 1 diff --git a/src/utils/poisson_solver/poisson_solver.jl b/src/utils/poisson_solver/poisson_solver.jl index 25538e9..f4a36fe 100644 --- a/src/utils/poisson_solver/poisson_solver.jl +++ b/src/utils/poisson_solver/poisson_solver.jl @@ -24,6 +24,7 @@ function solve_poisson_mgpcg!( kwargs... ) where {T, N} N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d")) + checkshape(u, d2u, (:u, :d2u)) checkshape(axes(mask), axes(u)[1:3], (:mask, :u)) @@ -90,12 +91,11 @@ function solve_poisson_dct!( dx::NTuple{3, Real} ) where {N} N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d")) + checkshape(u, d2u, (:u, :d2u)) nx, ny, nz = size(d2u)[1:3] - idx2 = inv(dx[1].*dx[1]) - idy2 = inv(dx[2].*dx[2]) - idz2 = inv(dx[3].*dx[3]) + idx2, idy2, idz2 = inv.(dx.*dx) # extreme slowdown for certain sizes with lots of threads # even worse for in-place, ie dct! @@ -109,7 +109,7 @@ function solve_poisson_dct!( Y = [2*(cospi(j)-1)*idy2 for j in range(0, step=1/ny, length=ny)] Z = [2*(cospi(k)-1)*idz2 for k in range(0, step=1/nz, length=nz)] - @inbounds for t in axes(d2û, 4) + for t in axes(d2û, 4) d2ût = @view(d2û[:,:,:,t]) @batch for k in 1:nz for j in 1:ny @@ -143,12 +143,11 @@ function solve_poisson_fft!( dx::NTuple{3, Real} ) where {N} N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d")) + checkshape(u, d2u, (:u, :d2u)) nx, ny, nz = size(d2u)[1:3] - idx2 = inv(dx[1].*dx[1]) - idy2 = inv(dx[2].*dx[2]) - idz2 = inv(dx[3].*dx[3]) + idx2, idy2, idz2 = inv.(dx.*dx) _rfft = iseven(nx) FFTW.set_num_threads(FFTW_NTHREADS[]) @@ -170,7 +169,7 @@ function solve_poisson_fft!( Y = [2*(cospi(j)-1)*idy2 for j in fftfreq(size(u, 2), 2)] Z = [2*(cospi(k)-1)*idz2 for k in fftfreq(size(u, 3), 2)] - @inbounds for t in axes(d2û, 4) + for t in axes(d2û, 4) d2ût = @view(d2û[:,:,:,t]) @batch for k in 1:nz for j in 1:ny diff --git a/src/utils/r2star.jl b/src/utils/r2star.jl index a524bc6..47fd439 100644 --- a/src/utils/r2star.jl +++ b/src/utils/r2star.jl @@ -38,7 +38,7 @@ function r2star_ll( b = tmap(log, transpose(vmag)) x̂ = ldiv!(A, b) @batch for I in eachindex(vr2s) - @inbounds vr2s[I] = x̂[1,I] + vr2s[I] = x̂[1,I] end else @@ -55,14 +55,14 @@ function r2star_ll( bt = @view(b[t,:]) mt = @view(vmag[:,t]) @batch for I in eachindex(R) - @inbounds bt[I] = log(mt[R[I]]) + bt[I] = log(mt[R[I]]) end end x̂ = ldiv!(A, b) @batch for I in eachindex(R) - @inbounds vr2s[R[I]] = x̂[1,I] + vr2s[R[I]] = x̂[1,I] end end @@ -120,7 +120,7 @@ function r2star_arlo( α = convert(T, 3 / (TEs[2]-TEs[1])) @batch for I in eachindex(vr2s) - @inbounds if mask === nothing || mask[I] + if mask === nothing || mask[I] m0 = vmag[I,1] m1 = vmag[I,2] m2 = vmag[I,3] @@ -227,7 +227,7 @@ function r2star_crsi( @batch for I in eachindex(vr2s) if mask === nothing || mask[I] den = β - @inbounds for t in 1:NT-1 + for t in 1:NT-1 p0 = vP[I,t] p1 = vP[I,t+1] p = pow(p0, γ0[1]) * pow(p1, γ1[1]) @@ -321,7 +321,7 @@ function r2star_numart2s( α = convert(T, 2*(NT - 1) / (TEs[end] - TEs[1])) @batch for I in eachindex(vr2s) - @inbounds if mask === nothing || mask[I] + if mask === nothing || mask[I] den = vmag[I,1] for t in 2:NT-1 den += vmag[I,t] + vmag[I,t] diff --git a/src/utils/utils.jl b/src/utils/utils.jl index b8fa9bc..3a71e65 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -163,15 +163,7 @@ function padarray!( pad == :replicate ? getindex_replicate : pad == :symmetric ? getindex_symmetric : pad == :reflect ? getindex_reflect : - throw(ArgumentError( - "pad must be one of " * - ":fill, " * - ":circular, " * - ":replicate, " * - ":symmetric, " * - ":reflect, " * - "got :$(pad)" - )) + checkopts(pad, (:fill, :circular, :replicate, :symmetric, :reflect), :pad) return _padarray_kernel!(xp, x, getindex_pad) end @@ -186,9 +178,9 @@ function _padarray_kernel!(xp::AbstractArray, x::AbstractArray, getindex_pad) @batch for Ip in CartesianIndices(xp) I = Ip - ΔI if any(map(∉, I.I, ax)) - @inbounds xp[Ip] = getindex_pad(x, I, lo, hi) + xp[Ip] = getindex_pad(x, I, lo, hi) else - @inbounds xp[Ip] = x[I] + xp[Ip] = x[I] end end @@ -399,7 +391,7 @@ function erode_mask!( @batch for k in 1+t:nz-t for j in 1+t:ny-t for i in 1+t:nx-t - @inbounds m1[i,j,k] = __erode_kernel(m0, i, j, k) + m1[i,j,k] = __erode_kernel(m0, i, j, k) end end end @@ -469,7 +461,7 @@ function psf2otf( else _kp = tfill!(similar(k, sz), zero(T)) @batch minbatch=1024 for I in CartesianIndices(k) - @inbounds _kp[I] = k[I] + _kp[I] = k[I] end end @@ -728,3 +720,23 @@ function checkshape( end return nothing end + + +checkopts(::Type{Bool}, o::T, opts::NTuple{N, T}) where {N, T} = o ∈ opts + +function checkopts( + o::T, + opts::NTuple{N, T}, + var::Union{Symbol, AbstractString} = :x +) where {N, T} + if !checkopts(Bool, o, opts) + pp = map((o, opts...)) do x + x isa Symbol ? ":$x" : + x isa AbstractString ? "\"$x\"" : "$x" + end + s1 = first(pp) + s2 = join(Base.tail(pp), ", ", " ") + throw(ArgumentError("$var must be one of $s2, got $s1")) + end + return nothing +end