Skip to content

Commit

Permalink
Improve error messages, formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
kamesy committed Aug 26, 2022
1 parent 0639948 commit 6fca4ce
Show file tree
Hide file tree
Showing 23 changed files with 267 additions and 572 deletions.
3 changes: 1 addition & 2 deletions src/QSM.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
14 changes: 7 additions & 7 deletions src/bgremove/ismv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

= 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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/bgremove/lbv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
21 changes: 10 additions & 11 deletions src/bgremove/pdf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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

= 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

Expand All @@ -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

= 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

Expand Down
58 changes: 29 additions & 29 deletions src/bgremove/sharp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

= 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
Expand All @@ -117,7 +117,7 @@ function _sharp!(

= 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̂)
Expand Down Expand Up @@ -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
= 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

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
= 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

Expand All @@ -391,15 +391,15 @@ 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]
end
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]
Expand Down
29 changes: 11 additions & 18 deletions src/inversion/direct.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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]))

= mul!(F̂, P, fp)
Expand Down Expand Up @@ -243,22 +236,22 @@ function _kdiv_ikernel!(
δ = convert(T, lambda)
= 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)

Expand All @@ -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
Expand Down
Loading

0 comments on commit 6fca4ce

Please sign in to comment.