API Reference

T2mapOptions(; kwargs...)
T2mapOptions(image::Array{T,4}; kwargs...) where {T}

Options structure for T2mapSEcorr. This struct collects keyword arguments passed to T2mapSEcorr, performs checks on parameter types and values, and assigns default values to unspecified parameters.

Arguments

• legacy::Bool: Perform T2-mapping using legacy algorithms. Default: false

• Threaded::Bool: Perform T2-mapping using multiple threads. Default: Threads.nthreads() > 1

• MatrixSize::Tuple{Int64, Int64, Int64}: Size of first 3 dimensions of input 4D image. This argument has no default, but is inferred automatically as size(image)[1:3] when calling T2mapSEcorr(image; kwargs...).

• nTE::Int64: Number of echoes in input signal. This argument is has no default, but is inferred automatically as size(image, 4) when calling T2mapSEcorr(image; kwargs...).

• TE::Real: Interecho spacing (Units: seconds). This argument has no default.

• nT2::Int64: Number of T2 times to estimate in the multi-exponential analysis. This argument has no default.

• T2Range::Tuple{T, T} where T<:Real: Tuple of min and max T2 values (Units: seconds). This argument has no default.

• T1::Real: Assumed value of T1 (Units: seconds). Default: 1.0

• Threshold::Real: First echo intensity cutoff for empty voxels. Default: 0.0

• MinRefAngle::Real: Minimum refocusing angle for flip angle optimization (Units: degrees). Default: 50.0

• nRefAngles::Int64: During flip angle optimization, goodness of fit is checked for up to nRefAngles angles in the range [MinRefAngle, 180]. The optimal angle is then determined through interpolation from these samples. Default: if !legacy 64 else 8 end

• nRefAnglesMin::Int64: Initial number of angles to check during flip angle optimization before refinement near likely optima. Setting nRefAnglesMin equal to nRefAngles forces all angles to be checked. Default: if !legacy min(5, nRefAngles) else nRefAngles end

• Reg::String: Regularization routine to use. One of "none", "lcurve", "gcv", "chi2", or "mdp", representing no regularization, the L-Curve method, the Generalized Cross-Validation method, Chi2Factor-based Tikhonov regularization, or the Morozov discrepancy principle, respectively.

• Chi2Factor::Union{Nothing, T} where T<:Real: Constraint on $\chi^2$ used for regularization when Reg == "chi2". Default: nothing

• NoiseLevel::Union{Nothing, T} where T<:Real: Estimate of the homoscedastic noise level $|b_i - \hat{b}_i|$, where $b$ is the unknown true signal and $\hat{b}$ is the measured corrupted signal. For Gaussian noise, this is the standard deviation. Default: nothing

• RefConAngle::Real: Refocusing pulse control angle (Units: degrees). Default: 180.0

• SetFlipAngle::Union{Nothing, T} where T<:Real: Instead of optimizing flip angle, use SetFlipAngle for all voxels (Units: degrees). Default: nothing

• SaveResidualNorm::Bool: Boolean flag to include a 3D array of the $\ell^2$-norms of the residuals from the NNLS fits in the output maps dictionary. Default: false

• SaveDecayCurve::Bool: Boolean flag to include a 4D array of the time domain decay curves resulting from the NNLS fits in the output maps dictionary. Default: false

• SaveRegParam::Bool: Boolean flag to include 3D arrays of the regularization parameters $\mu$ and resulting $\chi^2$-factors in the output maps dictionary. Default: false

• SaveNNLSBasis::Bool: Boolean flag to include a 5D (or 2D if SetFlipAngle is used) array of NNLS basis matrices in the output maps dictionary. Default: false

• Silent::Bool: Suppress printing to the console. Default: false

See also:

• T2mapSEcorr

T2mapSEcorr(image::Array{T,4}; <keyword arguments>)
T2mapSEcorr(image::Array{T,4}, opts::T2mapOptions{T})

Uses nonnegative least squares (NNLS) to compute T2 distributions in the presence of stimulated echos by optimizing the refocusing pulse flip angle. Records parameter maps and T2 distributions for further partitioning.

Arguments

• image: 4D array with intensity data as (row, column, slice, echo)
• A series of optional keyword argument settings which will be used to construct a T2mapOptions struct internally, or a T2mapOptions struct directly

Outputs

• maps: dictionary containing parameter maps with the following fields:

  • Default Fields

    • "echotimes" Echo times of time signal (length nTE 1D array)
    • "t2times" T2 times corresponding to T2-distributions (length nT2 1D array)
    • "refangleset" Refocusing angles used during flip angle optimization (length nRefAngles 1D array by default; scalar if SetFlipAngle is used)
    • "decaybasisset" Decay basis sets corresponding to "refangleset" (nTE x nT2 x nRefAngles 3D array by default; nTE x nT2 2D array if SetFlipAngle is used)
    • "gdn": Map of general density = sum(T2distribution) (MatrixSize 3D array)
    • "ggm": Map of general geometric mean of T2-distribution (MatrixSize 3D array)
    • "gva": Map of general variance (MatrixSize 3D array)
    • "fnr": Map of fit to noise ratio = gdn / √(sum(residuals.^2) / (nTE-1)) (MatrixSize 3D array)
    • "snr": Map of signal to noise ratio = maximum(signal) / std(residuals) (MatrixSize 3D array)
    • "alpha": Map of optimized refocusing pulse flip angle (MatrixSize 3D array)

  • Optional Fields

    • "resnorm": $\ell^2$-norm of NNLS fit residuals; see SaveResidualNorm option (MatrixSize 3D array)
    • "decaycurve": Signal decay curve resulting from NNLS fit; see SaveDecayCurve option (MatrixSize x nTE 4D array)
    • "mu": Regularization parameter used during from NNLS fit; see SaveRegParam option (MatrixSize 3D array)
    • "chi2factor": $\chi^2$ increase factor relative to unregularized NNLS fit; see SaveRegParam option (MatrixSize 3D array)
    • "decaybasis": Decay bases resulting from flip angle optimization; see SaveNNLSBasis option (MatrixSize x nTE x nT2 5D array, or nTE x nT2 2D array if SetFlipAngle is used)

• distributions: T2-distribution array with data as (row, column, slice, T2 amplitude) (MatrixSize x nT2 4D array)

Examples

julia> image = DECAES.mock_image(; MatrixSize = (100, 100, 1), nTE = 48); # mock image with size 100x100x1x48

julia> maps, dist = T2mapSEcorr(image; TE = 10e-3, nT2 = 40, T2Range = (10e-3, 2.0), Reg = "lcurve", Silent = true); # compute the T2-maps and T2-distribution

julia> maps
Dict{String, Any} with 10 entries:
  "echotimes"     => [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08,…
  "t2times"       => [0.01, 0.0114551, 0.013122, 0.0150315, 0.0172188…
  "refangleset"   => [50.0, 54.1935, 58.3871, 62.5806, 66.7742, 70.96…
  "gdn"           => [1.26381 1.27882 … 1.2463 1.25091; 1.29848 1.243…
  "fnr"           => [379.9 437.541 … 446.88 386.396; 485.27 360.591 …
  "alpha"         => [165.461 166.286 … 164.614 164.389; 163.735 164.…
  "gva"           => [0.691794 0.440231 … 0.0490302 0.1253; 0.849798 …
  "ggm"           => [0.0663333 0.0705959 … 0.056455 0.0576729; 0.053…
  "snr"           => [312.773 364.031 … 363.463 313.372; 372.631 313.…
  "decaybasisset" => [0.0277684 0.0315296 … 0.0750511 0.0751058; 0.04…

See also:

• T2partSEcorr
• lsqnonneg
• lsqnonneg_tikh
• lsqnonneg_lcurve
• lsqnonneg_gcv
• lsqnonneg_chi2
• lsqnonneg_mdp
• EPGdecaycurve

T2partOptions(; kwargs...)
T2partOptions(t2dist::Array{T,4}; kwargs...) where {T}

Options structure for T2partSEcorr. This struct collects keyword arguments passed to T2partSEcorr, performs checks on parameter types and values, and assigns default values to unspecified parameters.

Arguments

• legacy::Bool: Calculate T2-parts using legacy algorithms. Default: false

• Threaded::Bool: Perform T2-parts using multiple threads. Default: Threads.nthreads() > 1

• MatrixSize::Tuple{Int64, Int64, Int64}: Size of first 3 dimensions of input 4D T2 distribution. This argument is has no default, but is inferred automatically as size(t2dist)[1:3] when calling T2partSEcorr(t2dist; kwargs...).

• nT2::Int64: Number of T2 times to use. This argument has no default.

• T2Range::Tuple{T, T} where T<:Real: Tuple of min and max T2 values (Units: seconds). This argument has no default.

• SPWin::Tuple{T, T} where T<:Real: Tuple of min and max T2 values of the short peak window (Units: seconds). This argument has no default.

• MPWin::Tuple{T, T} where T<:Real: Tuple of min and max T2 values of the middle peak window (Units: seconds). This argument has no default.

• Sigmoid::Union{Nothing, T} where T<:Real: Apply sigmoidal weighting to the upper limit of the short peak window in order to smooth the hard small peak window cutoff time. Sigmoid is the delta-T2 parameter, which is the distance in seconds on either side of the SPWin upper limit where the sigmoid curve reaches 10% and 90% (Units: seconds). Default: nothing

• Silent::Bool: Suppress printing to the console. Default: false

See also:

• T2partSEcorr

T2partSEcorr(T2distributions::Array{T,4}; <keyword arguments>)
T2partSEcorr(T2distributions::Array{T,4}, opts::T2partOptions{T})

Analyzes T2 distributions produced by T2mapSEcorr to produce data maps of a series of parameters.

Arguments

• T2distributions: 4D array with data as (row, column, slice, T2 amplitude)
• A series of optional keyword argument settings which will be used to construct a T2partOptions struct internally, or a T2partOptions struct directly

Ouputs

• maps: a dictionary containing the following 3D data maps as fields:

  • "sfr": small pool fraction, e.g. myelin water fraction (MatrixSize 3D array)
  • "sgm": small pool geometric mean T2 (MatrixSize 3D array)
  • "mfr": medium pool fraction, e.g. intra/extracellular water fraction (MatrixSize 3D array)
  • "mgm": medium pool geometric mean T2 (MatrixSize 3D array)

Examples

julia> dist = DECAES.mock_t2dist(; MatrixSize = (100, 100, 1), nT2 = 40); # mock distribution with size 100x100x1x40

julia> maps = T2partSEcorr(dist; T2Range = (10e-3, 2.0), SPWin = (10e-3, 25e-3), MPWin = (25e-3, 200e-3), Silent = true); # compute T2-parts maps

julia> maps
Dict{String, Any} with 4 entries:
  "sgm" => [0.014202 0.0106354 … 0.0125409 0.0114035; 0.0119888 0.0110439 …
  "mfr" => [0.86938 0.886926 … 0.901487 0.835647; 0.840086 0.890914 … 0.88…
  "sfr" => [0.13062 0.112288 … 0.0985133 0.163075; 0.159914 0.109086 … 0.1…
  "mgm" => [0.0871951 0.0481156 … 0.0612596 0.0475037; 0.0629991 0.0738904…

See also:

• T2mapSEcorr

main(command_line_args = ARGS)

Entry point function for command line interface, parsing the command line arguments ARGS and subsequently calling one or both of T2mapSEcorr and T2partSEcorr with the parsed settings. See the Arguments section for available options.

See also:

• T2mapSEcorr
• T2partSEcorr

lsqnonneg(A::AbstractMatrix, b::AbstractVector)

Compute the nonnegative least-squares (NNLS) solution $X$ of the problem:

\[X = \underset{x \ge 0}{\operatorname{argmin}}\; ||Ax - b||_2^2.\]

Arguments

• A::AbstractMatrix: Left hand side matrix acting on x
• b::AbstractVector: Right hand side vector

Outputs

• X::AbstractVector: NNLS solution

lsqnonneg_tikh(A::AbstractMatrix, b::AbstractVector, μ::Real)

Compute the Tikhonov-regularized nonnegative least-squares (NNLS) solution $X_{\mu}$ of the problem:

\[X_{\mu} = \underset{x \ge 0}{\operatorname{argmin}}\; ||Ax - b||_2^2 + \mu^2 ||x||_2^2.\]

Arguments

• A::AbstractMatrix: Left hand side matrix acting on x
• b::AbstractVector: Right hand side vector
• μ::Real: Regularization parameter

Outputs

• X::AbstractVector: NNLS solution

lsqnonneg_lcurve(A::AbstractMatrix, b::AbstractVector)

Compute the Tikhonov-regularized nonnegative least-squares (NNLS) solution $X_{\mu}$ of the problem:

\[X_{\mu} = \underset{x \ge 0}{\operatorname{argmin}}\; ||Ax - b||_2^2 + \mu^2 ||L x||_2^2\]

where $L$ is the identity matrix, and $\mu$ is chosen by locating the corner of the "L-curve"[1]. Details of L-curve theory can be found in Hansen (1992)[2].

Arguments

• A::AbstractMatrix: Decay basis matrix
• b::AbstractVector: Decay curve data

Outputs

• X::AbstractVector: Regularized NNLS solution
• mu::Real: Resulting regularization parameter $\mu$
• chi2::Real: Resulting increase in residual norm relative to the unregularized $\mu = 0$ solution

References

1. A. Cultrera and L. Callegaro, "A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems". IOPSciNotes, vol. 1, no. 2, p. 025004, Aug. 2020, https://doi.org/10.1088/2633-1357/abad0d.
2. Hansen, P.C., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580, https://doi.org/10.1137/1034115.

lsqnonneg_gcv(A::AbstractMatrix, b::AbstractVector)

Compute the Tikhonov-regularized nonnegative least-squares (NNLS) solution $X_{\mu}$ of the problem:

\[X_{\mu} = \underset{x \ge 0}{\operatorname{argmin}}\; ||Ax - b||_2^2 + \mu^2 ||L x||_2^2\]

where $L$ is the identity matrix, and $\mu$ is chosen via the Generalized Cross-Validation (GCV) method:

\[\mu = \underset{\nu \ge 0}{\operatorname{argmin}}\; \frac{||AX_{\nu} - b||_2^2}{\mathcal{T}(\nu)^2}\]

where $\mathcal{T}(\mu)$ is the "degrees of freedom" of the regularized system

\[\mathcal{T}(\mu) = \operatorname{tr}(I - A (A^T A + \mu^2 L^T L) A^T).\]

Details of the GCV method can be found in Hansen (1992)[1].

Arguments

• A::AbstractMatrix: Decay basis matrix
• b::AbstractVector: Decay curve data

Outputs

• X::AbstractVector: Regularized NNLS solution
• mu::Real: Resulting regularization parameter $\mu$
• chi2::Real: Resulting increase in residual norm relative to the unregularized $\mu = 0$ solution

References

1. Hansen, P.C., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580, https://doi.org/10.1137/1034115.

lsqnonneg_chi2(A::AbstractMatrix, b::AbstractVector, chi2_target::Real)

Compute the Tikhonov-regularized nonnegative least-squares (NNLS) solution $X_{\mu}$ of the problem:

\[X_{\mu} = \underset{x \ge 0}{\operatorname{argmin}}\; ||Ax - b||_2^2 + \mu^2 ||x||_2^2\]

where $\mu$ is determined by solving:

\[\chi^2(\mu) = \frac{||AX_{\mu} - b||_2^2}{||AX_{0} - b||_2^2} = \chi^2_{\mathrm{target}}.\]

That is, $\mu$ is chosen such that the squared residual norm of the regularized problem is chi2_target times larger than the squared residual norm of the unregularized problem.

Arguments

• A::AbstractMatrix: Decay basis matrix
• b::AbstractVector: Decay curve data
• chi2_target::Real: Target $\chi^2(\mu)$; typically a small value, e.g. 1.02 representing a 2% increase

Outputs

• X::AbstractVector: Regularized NNLS solution
• mu::Real: Resulting regularization parameter $\mu$
• chi2::Real: Resulting $\chi^2(\mu)$, which should be approximately equal to chi2_target

lsqnonneg_mdp(A::AbstractMatrix, b::AbstractVector, δ::Real)

Compute the Tikhonov-regularized nonnegative least-squares (NNLS) solution $X_{\mu}$ of the problem:

\[X_{\mu} = \underset{x \ge 0}{\operatorname{argmin}}\; ||Ax - b||_2^2 + \mu^2 ||x||_2^2\]

where $\mu$ is chosen using Morozov's Discrepency Principle (MDP)[1,2]:

\[\mu = \operatorname{sup}\; \left\{ \nu \ge 0 : ||AX_{\nu} - b|| \le \delta \right\}.\]

That is, $\mu$ is maximized subject to the constraint that the residual norm of the regularized problem is at most $\delta$[1].

Arguments

• A::AbstractMatrix: Decay basis matrix
• b::AbstractVector: Decay curve data
• δ::Real: Upper bound on regularized residual norm

Outputs

• X::AbstractVector: Regularized NNLS solution
• mu::Real: Resulting regularization parameter $\mu$
• chi2::Real: Resulting increase in residual norm relative to the unregularized $\mu = 0$ solution

References

1. Morozov VA. Methods for Solving Incorrectly Posed Problems. Springer Science & Business Media, 2012.
2. Clason C, Kaltenbacher B, Resmerita E. Regularization of Ill-Posed Problems with Non-negative Solutions. In: Bauschke HH, Burachik RS, Luke DR (eds) Splitting Algorithms, Modern Operator Theory, and Applications. Cham: Springer International Publishing, pp. 113–135.

lcurve_corner(f, xlow, xhigh)

Find the corner of the L-curve via curvature maximization using a modified version of Algorithm 1 from Cultrera and Callegaro (2020)[1].

References

1. A. Cultrera and L. Callegaro, "A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems". IOPSciNotes, vol. 1, no. 2, p. 025004, Aug. 2020, https://doi.org/10.1088/2633-1357/abad0d.

EPGdecaycurve(ETL::Int, α::Real, TE::Real, T2::Real, T1::Real, β::Real)

Computes the normalized echo decay curve for a multi spin echo sequence using the extended phase graph algorithm using the given input parameters.

The sequence of flip angles used is slight generalization of the standard 90 degree excitation pulse followed by 180 degree pulse train. Here, the sequence used is A*90, A*180, A*β, A*β, ... where A = α/180 accounts for B1 inhomogeneities. Equivalently, the pulse sequence can be written as α/2, α, α * (β/180), α * (β/180), .... Note that if α = β = 180, we recover the standard 90, 180, 180, ... pulse sequence.

Arguments

• ETL::Int: echo train length, i.e. number of echos
• α::Real: angle of refocusing pulses (Units: degrees)
• TE::Real: inter-echo time (Units: seconds)
• T2::Real: transverse relaxation time (Units: seconds)
• T1::Real: longitudinal relaxation time (Units: seconds)
• β::Real: value of Refocusing Pulse Control Angle (Units: degrees)

Outputs

• decay_curve::AbstractVector: normalized echo decay curve with length ETL