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

• 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 has no default, but is inferred automatically as size(image, 4) when calling T2mapSEcorr(image; kwargs...).

• TE::Real: Interecho spacing (Units: time, must match T1 and T2Range). 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: time, must match TE). This argument has no default.

• T1::Real: Assumed value of T1 (Units: time, must match TE). 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: 90.0

• nRefAngles::Int64: Maximum number of candidate angles spanning [MinRefAngle, 180] during flip angle optimization. Default: 500

• nRefAnglesMin::Int64: Initial number of angles evaluated before refinement near likely optima. Setting nRefAnglesMin equal to nRefAngles evaluates every candidate angle. Default: min(15, nRefAngles)

• Reg::String: Regularization routine to use. One of "gcv", "lcurve", "reginska", "chi2", "mdp", or "none", representing the Generalized Cross-Validation method, the L-Curve method, Regińska's minimum-product criterion, Chi2Factor-based regularization, the Morozov discrepancy principle, or no regularization, respectively.

• RegNorm::String: Regularization penalty norm. One of "l2" or "l1", penalizing $\mu^2 ||x||_2^2$ or $\mu ||x||_1$, respectively. Ignored when Reg == "none". Default: l2

• 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) (Units: same as input signal) (MatrixSize 3D array)
    • "ggm": Map of general geometric mean of T2-distribution (Units: time, must match T2Range) (MatrixSize 3D array)
    • "gva": Map of general variance of the T2-distribution (Units: none) (MatrixSize 3D array)
    • "fnr": Map of fit to noise ratio = gdn / √(sum(residuals.^2) / (nTE-1)) (Units: none) (MatrixSize 3D array)
    • "snr": Map of signal to noise ratio = maximum(signal) / std(residuals) (Units: none) (MatrixSize 3D array)
    • "alpha": Map of optimized refocusing pulse flip angle (Units: degrees) (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 the 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"   => [90.0, 90.1804, 90.3607, 90.5411, 90.7214, 90.9018…
  "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_lasso
• lsqnonneg_gcv
• lsqnonneg_lcurve
• lsqnonneg_lcurve_lasso
• lsqnonneg_reginska
• lsqnonneg_reginska_lasso
• lsqnonneg_chi2
• lsqnonneg_chi2_lasso
• lsqnonneg_mdp
• lsqnonneg_mdp_lasso
• 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

• 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 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: time, must match the T2 distribution). 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: time, must match T2Range). 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: time, must match T2Range). 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 on either side of the SPWin upper limit where the sigmoid curve reaches 10% and 90% (Units: time, must match T2Range). 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

Outputs

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

  • "sfr": small pool fraction, e.g. myelin water fraction (Units: none, in [0, 1]) (MatrixSize 3D array)
  • "sgm": small pool geometric mean T2 (Units: time, matching T2Range) (MatrixSize 3D array)
  • "mfr": medium pool fraction, e.g. intra/extracellular water fraction (Units: none, in [0, 1]) (MatrixSize 3D array)
  • "mgm": medium pool geometric mean T2 (Units: time, matching T2Range) (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::Vector{String} = 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

load_image(filename; ndims::Int = 4)

Load an image from filename as an ndims-dimensional Array{Float64}. Supported file types are .mat, .nii, .nii.gz, .par, .xml, and .rec.

NIfTI values are scaled as scl_slope * x + scl_inter, unless scl_slope == 0, in which case the raw values are returned. For MAT files, an array of the requested dimensionality must be present; if several are found, the first in sorted order is used and a warning is issued. An array with fewer than ndims dimensions is reshaped with additional trailing singleton dimensions; an array with more than ndims dimensions is truncated to the first volume along its trailing dimensions.

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

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

\[x_{0} = \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_lasso(A::AbstractMatrix, b::AbstractVector, μ::Real)

Compute the $\ell^1$-regularized nonnegative least-squares (NNLS) solution $x_{\mu}$ of the problem:

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

Nonnegativity makes $||x||_1 = \mathbf{1}^T x$ linear, so this is a smooth bound-constrained convex quadratic program, solved by a Lawson-Hanson active-set method.

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_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_lcurve(A::AbstractMatrix, b::AbstractVector; max_slope = $(LCURVE_SLOPE_MAX_DEFAULT))

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 at a corner of the "L-curve"[1], a local maximum of the curvature of $\mu \mapsto (\log||Ax_\mu - b||_2^2, \log||x_\mu||_2^2)$. The L-curve may have several corners; max_slope excludes those in its near-vertical $\mu \to 0$ tail, where the fit is barely regularized. If no corner is found, $\mu = 0$ and the unregularized solution is returned. Details of L-curve theory can be found in Hansen (1992)[2].

Arguments

• A::AbstractMatrix: Decay basis matrix
• b::AbstractVector: Decay curve data
• max_slope::Real = $(LCURVE_SLOPE_MAX_DEFAULT): reject corners at which $||Ax_\mu - b||_2^2 / (\mu^2 ||x_\mu||_2^2)$ exceeds max_slope. Pass Inf to accept any corner.

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_lcurve_lasso(A::AbstractMatrix, b::AbstractVector; max_slope = $(LCURVE_SLOPE_MAX_DEFAULT))

Compute the $\ell^1$-regularized nonnegative least-squares (NNLS) solution $x_{\mu}$ of the problem:

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

where $\mu$ is chosen at a corner of the $\ell^1$ "L-curve"[1] $\mu \mapsto (\log||Ax_\mu - b||_2^2, 2\log||x_\mu||_1)$, the first positive local maximum of its turning rate $\omega = d\theta/d\log\mu$ that max_slope admits.

As for lsqnonneg_lcurve, max_slope excludes corners in the near-vertical $\mu \to 0$ tail. If no corner is found, $\mu = 0$ and the unregularized solution is returned.

Arguments

• A::AbstractMatrix: Decay basis matrix
• b::AbstractVector: Decay curve data
• max_slope::Real = $(LCURVE_SLOPE_MAX_DEFAULT): reject corners at which $2 ||Ax_\mu - b||_2^2 / (\mu ||x_\mu||_1)$ exceeds max_slope. Pass Inf to accept any corner.

Outputs

• x::AbstractVector: Regularized NNLS solution
• mu::Real: Resulting regularization parameter $\mu$
• chi2::Real: Resulting ratio $||Ax_{\mu} - b||_2^2 / ||Ax_0 - b||_2^2$ of squared residual norms

References

1. J. Nasehi Tehrani, A. McEwan, C. Jin and A. van Schaik, "L1 regularization method in electrical impedance tomography by using the L1-curve (Pareto frontier curve)". Applied Mathematical Modelling, 36(3), 1095-1105, 2012, https://doi.org/10.1016/j.apm.2011.07.055.

lsqnonneg_reginska(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 ||x||_2^2\]

where $\mu$ is chosen by Regińska's minimum-product criterion[1]:

\[\mu = \underset{\nu > 0}{\operatorname{argmin}}\; \Psi(\nu) = ||Ax_{\nu} - b||_2^2 \, ||x_{\nu}||_2^2,\]

taking the smallest local minimizer of $\Psi$. Stationarity of $\Psi$ is equivalent to a log-log L-curve tangent slope of $-1$, so the selected $\mu$ is the balance point $||Ax_{\mu} - b|| = \mu ||x_{\mu}||$. The smallest local minimizer is taken because $\Psi \to 0$ trivially as $\mu \to \infty$, where $x_{\mu} \to 0$. If $\Psi$ has no interior local minimum, the unregularized solution is returned with $\mu = 0$.

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. T. Regińska, "A Regularization Parameter in Discrete Ill-Posed Problems". SIAM Journal on Scientific Computing, 17(3), 740-749, 1996, https://doi.org/10.1137/S1064827593252672.

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

Compute the $\ell^1$-regularized nonnegative least-squares (NNLS) solution $x_{\mu}$ of the problem:

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

where $\mu$ is the smallest positive local minimizer of an $\ell^1$ analogue of Regińska's minimum-product criterion[1]:

\[\Psi(\nu) = ||Ax_{\nu} - b||_2^2 \, ||x_{\nu}||_1.\]

This is the $\ell^1$ counterpart of lsqnonneg_reginska. Stationarity of $\Psi$ is equivalent to a log-log tangent slope of $-1$ for the pair $(||Ax_{\nu} - b||_2^2, ||x_{\nu}||_1)$, so the selected $\mu$ is the balance point $||Ax_{\mu} - b||_2^2 = \mu ||x_{\mu}||_1$ at which the two terms of the objective contribute equally.

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 ratio $||Ax_{\mu} - b||_2^2 / ||Ax_0 - b||_2^2$ of squared residual norms

References

1. T. Regińska, "A Regularization Parameter in Discrete Ill-Posed Problems". SIAM Journal on Scientific Computing, 17(3), 740-749, 1996, https://doi.org/10.1137/S1064827593252672.

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_chi2_lasso(A::AbstractMatrix, b::AbstractVector, chi2_target::Real)

Compute the $\ell^1$-regularized nonnegative least-squares (NNLS) solution $x_{\mu}$ of the problem:

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

where $\mu$ is determined by solving:

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

This is the $\ell^1$ counterpart of lsqnonneg_chi2.

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 Discrepancy 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.

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

Compute the $\ell^1$-regularized nonnegative least-squares (NNLS) solution $x_{\mu}$ of the problem:

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

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

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

This is the $\ell^1$ counterpart of lsqnonneg_mdp.

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 ratio $||Ax_{\mu} - b||_2^2 / ||Ax_0 - b||_2^2$ of squared residual norms

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, t₀, P₀, sig₀; kwargs...)

Locate a corner of the L-curve, following Cultrera and Callegaro (2020)[1] with a dynamically bracketed search.

f(t) returns the L-curve point at $t = \log\mu$ with its curvature, angular velocity, and active-set digest; P₀ and sig₀ correspond to the unregularized solution at which the curve terminates. Returns $\log\mu$ at a corner, or NaN if none is found. A returned corner has positive curvature, is a local maximum, and satisfies the max_slope guard.

Keyword arguments

• init_width: width in $\log\mu$ of the initial bracket, and the scale at which the first corner is sought.
• xtol, Ptol: absolute tolerances on $\log\mu$ and on the chord between log-log curve evaluation points.
• max_expand, max_backtrack, nsweep, max_candidates: search budgets.
• max_slope: reject corners whose log-log tangent slope exceeds this.
• bounds: interval of $\log\mu$ to search within.

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: time, must match T1 and T2)
• T2::Real: transverse relaxation time (Units: time, must match TE)
• T1::Real: longitudinal relaxation time (Units: time, must match TE)
• β::Real: refocusing pulse control angle (Units: degrees)

Outputs

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

The four-argument method omits β, fixing the refocusing control angle at 180 degrees, which is the standard CPMG sequence.

EPGdecaycurve!(work::AbstractEPGWorkspace, θ::EPGParameterization)
EPGdecaycurve!(decay_curve::AbstractVector, work::AbstractEPGWorkspace, θ::EPGParameterization)

In-place form of EPGdecaycurve, reusing the workspace work allocated by EPGdecaycurve_work. The two-argument method writes into and returns the decay curve held by work; the three-argument method writes into decay_curve.

See also:

• EPGdecaycurve
• EPGdecaycurve_work

EPGdecaycurve_work(θ::EPGParameterization)
EPGdecaycurve_work(T::Type, ETL::Int)
EPGdecaycurve_work(T::Type, ::Val{ETL})

Allocate a reusable workspace for EPGdecaycurve! with element type T and echo train length ETL. Passing ETL as a Val selects the workspace specialized for a statically known echo train length.

See also:

• EPGdecaycurve!