From ac8e978b47f24216aef23fd4bf6012881228b791 Mon Sep 17 00:00:00 2001 From: kamesy Date: Wed, 13 Apr 2022 16:11:00 -0700 Subject: [PATCH] Add R2* mapping --- docs/src/api.md | 8 + src/QSM.jl | 3 +- src/utils/r2star.jl | 413 ++++++++++++++++++++++++++++++++++++++++++++ src/utils/utils.jl | 3 +- 4 files changed, 425 insertions(+), 2 deletions(-) create mode 100644 src/utils/r2star.jl diff --git a/docs/src/api.md b/docs/src/api.md index ce87c75..4f22491 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -59,6 +59,14 @@ fit_echo_linear! ``` ## Other +### R2* mapping +```@docs +r2star_arlo +r2star_crsi +r2star_ll +r2star_numart2s +``` + ### Finite differences ```@docs gradfp diff --git a/src/QSM.jl b/src/QSM.jl index 754a41f..582c4ec 100644 --- a/src/QSM.jl +++ b/src/QSM.jl @@ -8,7 +8,7 @@ using NIfTI: NIVolume, niread, niwrite using Polyester: @batch, num_cores using PolyesterWeave: reset_workers! using Printf: @printf -using SLEEFPirates: sincos_fast +using SLEEFPirates: pow, sincos_fast using StaticArrays: SVector using ThreadingUtilities: initialize_task using TiledIteration: EdgeIterator, TileIterator, padded_tilesize @@ -21,6 +21,7 @@ export bet export gradfp, gradfp!, gradfp_adj, gradfp_adj!, lap, lap! export dipole_kernel, laplace_kernel, smv_kernel export fit_echo_linear, fit_echo_linear! +export r2star_arlo, r2star_crsi, r2star_ll, r2star_numart2s export crop_mask, crop_indices, erode_mask, erode_mask! export fastfftsize, padfastfft, padarray!, unpadarray, unpadarray!, psf2otf include("utils/utils.jl") diff --git a/src/utils/r2star.jl b/src/utils/r2star.jl new file mode 100644 index 0000000..8c92980 --- /dev/null +++ b/src/utils/r2star.jl @@ -0,0 +1,413 @@ +""" + function r2star_ll( + mag::AbstractArray{<:AbstractFloat, N > 1}, + TEs::NTuple{NT > 1, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing + ) -> typeof(similar(mag, size(mag)[1:N-1])) + +Log-linear fit. + +### Arguments +- `mag::AbstractArray{<:AbstractFloat, N > 1}`: multi-echo magnitude +- `TEs::NTuple{NT > 1, Real}`: echo times +- `mask::Union{Nothing, AbstractArray{Bool}} = nothing`: binary mask of region of interest + +### Returns +- `typeof(similar(mag, size(mag)[1:N-1]))`: R2* map (1 / units of TEs) +""" +function r2star_ll( + mag::AbstractArray{T, N}, + TEs::NTuple{NT, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing +) where {T<:AbstractFloat, N, NT} + N > 1 || throw(ArgumentError("array must contain echoes in last dimension")) + NT > 1 || throw(ArgumentError("data must be multi-echo")) + + size(mag, N) == NT || throw(DimensionMismatch()) + mask === nothing || length(mask) == length(mag) ÷ NT || throw(DimensionMismatch()) + + r2s = tzero(mag, size(mag)[1:N-1]) + + vmag = reshape(mag, :, NT) + vr2s = vec(r2s) + + A = qr(Matrix{T}([-[TEs...] ones(NT)])) + + if mask === nothing + b = tcopy(log, transpose(vmag)) + + x̂ = ldiv!(A, b) + + @inbounds @batch for I in eachindex(vr2s) + vr2s[I] = x̂[1,I] + end + + else + vmask = vec(mask) + + i = 0 + R = Vector{Int}(undef, sum(vmask)) + @inbounds for I in eachindex(vmask) + if vmask[I] + R[i += 1] = I + end + end + + b = similar(mag, (NT, length(R))) + @inbounds for t in Base.OneTo(NT) + bt = @view(b[t,:]) + mt = @view(vmag[:,t]) + @batch for I in eachindex(R) + bt[I] = log(mt[R[I]]) + end + end + + x̂ = ldiv!(A, b) + + @inbounds @batch for I in eachindex(R) + vr2s[R[I]] = x̂[1,I] + end + end + + return r2s +end + + +""" + function r2star_arlo( + mag::AbstractArray{<:AbstractFloat, N > 1}, + TEs::NTuple{NT > 1, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing + ) -> typeof(similar(mag, size(mag)[1:N-1])) + +Auto-Regression on Linear Operations (ARLO) [1]. + +### Arguments +- `mag::AbstractArray{<:AbstractFloat, N > 1}`: multi-echo magnitude +- `TEs::NTuple{NT > 1, Real}`: echo times +- `mask::Union{Nothing, AbstractArray{Bool}} = nothing`: binary mask of region of interest + +### Returns +- `typeof(similar(mag, size(mag)[1:N-1]))`: R2* map (1 / units of TEs) + +### References +[1] Pei M, Nguyen TD, Thimmappa ND, Salustri C, Dong F, Cooper MA, Li J, + Prince MR, Wang Y. Algorithm for fast monoexponential fitting based on + auto‐regression on linear operations (ARLO) of data. + Magnetic resonance in medicine. 2015 Feb;73(2):843-50. +""" +function r2star_arlo( + mag::AbstractArray{T, N}, + TEs::NTuple{NT, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing +) where {T<:AbstractFloat, N, NT} + N > 1 || throw(ArgumentError("array must contain echoes in last dimension")) + NT > 2 || throw(ArgumentError("ARLO requires at least 3 echoes")) + + size(mag, N) == NT || throw(DimensionMismatch()) + mask === nothing || length(mask) == length(mag) ÷ NT || throw(DimensionMismatch()) + + all((≈)(TEs[2]-TEs[1]), TEs[2:end].-TEs[1:end-1]) || + throw(DomainError("ARLO requires equidistant echoes")) + + r2s = tzero(mag, size(mag)[1:N-1]) + + vmag = reshape(mag, :, NT) + vr2s = vec(r2s) + + _zeroT = zero(T) + _twoT = convert(T, 2) + _fourT = convert(T, 4) + + α = convert(T, 3 / (TEs[2]-TEs[1])) + + if mask === nothing + @inbounds @batch for I in eachindex(vr2s) + m0 = vmag[I,1] + m1 = vmag[I,2] + m2 = vmag[I,3] + + δ = m0 - m2 + s = m0 + muladd(_fourT, m1, m2) + a = muladd(_twoT, m1, m0) + + num = δ * a + den = s * a + + for t in 2:NT-2 + m0 = m1 + m1 = m2 + m2 = vmag[I,t+2] + + δ = m0 - m2 + s = m0 + muladd(_fourT, m1, m2) + a = muladd(_twoT, m1, m0) + + num = muladd(δ, a, num) + den = muladd(s, a, den) + end + + vr2s[I] = iszero(den) ? _zeroT : α * num * inv(den) + end + + else + vmask = vec(mask) + + @inbounds @batch for I in eachindex(vr2s) + if vmask[I] + m0 = vmag[I,1] + m1 = vmag[I,2] + m2 = vmag[I,3] + + δ = m0 - m2 + s = m0 + muladd(_fourT, m1, m2) + a = muladd(_twoT, m1, m0) + + num = δ * a + den = s * a + + for t in 2:NT-2 + m0 = m1 + m1 = m2 + m2 = vmag[I,t+2] + + δ = m0 - m2 + s = m0 + muladd(_fourT, m1, m2) + a = muladd(_twoT, m1, m0) + + num = muladd(δ, a, num) + den = muladd(s, a, den) + end + + vr2s[I] = iszero(den) ? _zeroT : α * num * inv(den) + end + end + end + + return r2s +end + + +""" + function r2star_crsi( + mag::AbstractArray{<:AbstractFloat, N > 1}, + TEs::NTuple{NT > 1, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing; + M::Integer = 3, + sigma::Union{Nothing, Real} = nothing, + Rsz::NTuple{N-1, Integer} = size(mag)[1:N-1] .÷ 20, + ) -> typeof(similar(mag, size(mag)[1:N-1])) + +Calculation of Relaxivities by Signal Integration (CRSI) [1]. + +### Arguments +- `mag::AbstractArray{<:AbstractFloat, N > 1}`: multi-echo magnitude +- `TEs::NTuple{NT > 1, Real}`: echo times +- `mask::Union{Nothing, AbstractArray{Bool}} = nothing`: binary mask of region of interest + +### Keywords +- `M::Integer = 3`: interpolation factor +- `sigma::Union{Nothing, Real} = nothing`: noise +- `Rsz::NTuple{N-1, Integer} = size(mag)[1:N-1] .÷ 20`: + - `sigma isa Real`: unused + - `sigma isa Nothing`: size of kernels used to calculate the noise from the + background signal of the magnitude. + +### Returns +- `typeof(similar(mag, size(mag)[1:N-1]))`: R2* map (1 / units of TEs) + +### References +[1] Song R, Loeffler RB, Holtrop JL, McCarville MB, Hankins JS, Hillenbrand CM. + Fast quantitative parameter maps without fitting: Integration yields + accurate mono‐exponential signal decay rates. + Magnetic resonance in medicine. 2018 Jun;79(6):2978-85. +""" +function r2star_crsi( + mag::AbstractArray{T, N}, + TEs::NTuple{NT, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing; + M::Integer = 3, + sigma::Union{Nothing, Real} = nothing, + Rsz::NTuple{NR, Integer} = size(mag)[1:N-1] .÷ 20, +) where {T<:AbstractFloat, N, NT, NR} + N > 1 || throw(ArgumentError("array must contain echoes in last dimension")) + NT > 1 || throw(ArgumentError("data must be multi-echo")) + + size(mag, N) == NT || throw(DimensionMismatch()) + mask === nothing || length(mask) == length(mag) ÷ NT || throw(DimensionMismatch()) + sigma === nothing || NR == N-1 || throw(DimensionMismatch()) + M > 0 || throw(ArgumentError("interpolation factor M must be greater than 0")) + + P = tcopy(x -> x*x, mag) + r2s = tzero(mag, size(mag)[1:N-1]) + + vP = reshape(P, :, NT) + vr2s = vec(r2s) + + if sigma !== nothing + σ2 = convert(T, 2*sigma*sigma) + else + σ2 = _noise_crsi(P, Rsz, mask) + end + + _zeroT = zero(T) + α = convert(T, 1//2) + β = convert(T, -α * σ2 * (TEs[end] - TEs[1])) + + τ = SVector{NT-1, T}((TEs[2:end] .- TEs[1:end-1]) ./ (M+1)) + γ0 = SVector{M+1, T}([(2*M - 2*m + 1) / (2*M + 2) for m in 0:M]...) + γ1 = SVector{M+1, T}([(2*m + 1) / (2*M + 2) for m in 0:M]...) + + if mask === nothing + @inbounds @batch per=thread for I in eachindex(vr2s) + den = β + for t in Base.OneTo(NT-1) + p0 = vP[I,t] + p1 = vP[I,t+1] + p = pow(p0, γ0[1]) * pow(p1, γ1[1]) + for m in 2:M+1 + p = muladd(pow(p0, γ0[m]), pow(p1, γ1[m]), p) + end + den = muladd(τ[t], p, den) + end + + vr2s[I] = iszero(den) ? _zeroT : α * (vP[I,1] - vP[I,end]) * inv(den) + end + + else + vmask = vec(mask) + + @inbounds @batch per=thread for I in eachindex(vr2s) + if vmask[I] + den = β + for t in Base.OneTo(NT-1) + p0 = vP[I,t] + p1 = vP[I,t+1] + p = pow(p0, γ0[1]) * pow(p1, γ1[1]) + for m in 2:M+1 + p = muladd(pow(p0, γ0[m]), pow(p1, γ1[m]), p) + end + den = muladd(τ[t], p, den) + end + + vr2s[I] = iszero(den) ? _zeroT : α * (vP[I,1] - vP[I,end]) * inv(den) + end + end + end + + return r2s +end + +function _noise_crsi( + P::AbstractArray{T, N}, + Rsz::NTuple{M, Integer}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing +) where {T<:AbstractFloat, N, M} + M == N-1 || throw(DimensionMismatch()) + + sz = size(P) + rsz = ntuple(n -> min(Rsz[n], (sz[n]-2)÷2), Val(N-1)) + + lb = rsz .+ 2 + ub = sz[1:M] .- 1 .- rsz + + outer = CartesianIndices(ntuple(n -> 2:sz[n]-1, Val(N-1))) + inner = CartesianIndices(ntuple(n -> 2+rsz[n]:sz[n]-1-rsz[n], Val(N-1))) + + n = 0 + σ2 = zero(T) + + if mask === nothing + @inbounds for t in axes(P, N) + for I in EdgeIterator(outer, inner) + if all(n -> I[n] < lb[n] || I[n] > ub[n], 1:M) + σ2 += P[I,t] + n += 1 + end + end + end + + else + @inbounds for t in axes(P, N) + for I in EdgeIterator(outer, inner) + if all(n -> I[n] < lb[n] || I[n] > ub[n], 1:M) && !mask[I] + σ2 += P[I,t] + n += 1 + end + end + end + end + + return n == 0 ? zero(T) : σ2 / n +end + + +""" + function r2star_numart2s( + mag::AbstractArray{<:AbstractFloat, N > 1}, + TEs::NTuple{NT > 1, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing + ) -> typeof(similar(mag, size(mag)[1:N-1])) + +Numerical Algorithm for Real-time T2* mapping (NumART2*) [1]. + +### Arguments +- `mag::AbstractArray{<:AbstractFloat, N > 1}`: multi-echo magnitude +- `TEs::NTuple{NT > 1, Real}`: echo times +- `mask::Union{Nothing, AbstractArray{Bool}} = nothing`: binary mask of region of interest + +### Returns +- `typeof(similar(mag, size(mag)[1:N-1]))`: R2* map (1 / units of TEs) + +### References +[1] Hagberg GE, Indovina I, Sanes JN, Posse S. Real‐time quantification of T2* + changes using multiecho planar imaging and numerical methods. + Magnetic Resonance in Medicine: An Official Journal of the International + Society for Magnetic Resonance in Medicine. 2002 Nov;48(5):877-82. +""" +function r2star_numart2s( + mag::AbstractArray{T, N}, + TEs::NTuple{NT, Real}, + mask::Union{Nothing, AbstractArray{Bool}} = nothing +) where {T<:AbstractFloat, N, NT} + N > 1 || throw(ArgumentError("array must contain echoes in last dimension")) + NT > 1 || throw(ArgumentError("data must be multi-echo")) + + size(mag, N) == NT || throw(DimensionMismatch()) + mask === nothing || length(mask) == length(mag) ÷ NT || throw(DimensionMismatch()) + + r2s = tzero(mag, size(mag)[1:N-1]) + + vmag = reshape(mag, :, NT) + vr2s = vec(r2s) + + _zeroT = zero(T) + α = convert(T, 2*(NT - 1) / (TEs[end] - TEs[1])) + + if mask === nothing + @inbounds @batch for I in eachindex(vr2s) + den = vmag[I,1] + for t in 2:NT-1 + den += vmag[I,t] + vmag[I,t] + end + den += vmag[I,NT] + vr2s[I] = iszero(den) ? _zeroT : α * (vmag[I,1] - vmag[I,NT]) * inv(den) + end + + else + vmask = vec(mask) + + @inbounds @batch for I in eachindex(vr2s) + if vmask[I] + den = vmag[I,1] + for t in 2:NT-1 + den += vmag[I,t] + vmag[I,t] + end + den += vmag[I,NT] + vr2s[I] = iszero(den) ? _zeroT : α * (vmag[I,1] - vmag[I,NT]) * inv(den) + end + end + end + + return r2s +end diff --git a/src/utils/utils.jl b/src/utils/utils.jl index e3c7812..e26dd84 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -4,6 +4,7 @@ include("kernels.jl") include("lsmr.jl") include("multi_echo.jl") include("poisson_solver/poisson_solver.jl") +include("r2star.jl") ##### @@ -343,7 +344,7 @@ erode_mask(mask::AbstractArray{Bool, 3}, iter::Integer = 1) = erode_mask!(tzero(mask), mask, iter) """ - erode_mask( + erode_mask!( emask::AbstractArray{Bool, 3}, mask::AbstractArray{Bool, 3}, iter::Integer = 1