estimateIRsolid: Estimate parameters in Inversion Recovery MRI experiments...

In qMRI: Methods for Quantitative Magnetic Resonance Imaging ('qMRI')

Estimate parameters in Inversion Recovery MRI experiments mixture model for non-fluid voxel

Description

The Inversion Recovery MRI signal in non-fluid voxel follows is modeled as a mixture of a fluid and a solid compartment.

Usage

estimateIRsolid(IRfluidobj, TEScale = 100, dataScale = 1000,
verbose = TRUE, lower = c(0, 0, 0), upper = c(0.95, 2, 2))

Arguments

IRfluidobj	Object of class "IRfluid" as generated by function estimateIRfluid.
TEScale	Internal scale factor for Echo Times. This influences parameter scales in numerical calculations.
dataScale	Internal scale factor for MR signals. This influences parameter scales in numerical calculations.
verbose	Logical. Provide some runtime diagnostics.
lower	Lower bounds for parameter values.
upper	Upper bounds for parameter values.

Details

The Inversion Recovery MRI signal in non-fluid voxel follows is modeled as a mixture of a fluid and a solid compartment. The function calculates estimates of the maximum signal and recovery rate for the solid compartment and a mixture coefficient (proportion of fluid) for all voxel with segment%in%2:3using results from function estimateIRfluid.

Value

List of class IRmixed with components

IRdata	4D array containing the IRMRI data, first dimension refers to inversion times
InvTimes	vector of inversion times
segm	segmentation codes, 1 for CSF, 2 for GM, 3 for WM, 0 for out of brain
sigma	noise standard deviation, if not specified estimated fron CSF areas in image with largest inversion time
L	effective number of coils
fx	Array of fluid proportions
Sx	Array of maximal signals
Rx	Array of relaxation rates
Sf	Global estimate of maximal fluid signal
Rf	Global estimate of fluid relaxation rate
ICovx	Covariance matrix of estimates fx, Sx and Rx.
sigma	Array of provided or estimated noise standard deviations
Convx	Array of convergence indicators
rsdx	Residual standard deviations
method	"NLS" for nonlinear regression or "QL" for quasi likelihood.
varest	Method used for variance estimation

The arrays contain entries for all voxel with segments%in%1:3.

Author(s)

Karsten Tabelow [e-mail]
J\"org Polzehl [e-mail]

See Also

estimateIRfluid, estimateIR, estimateIRsolidfixed,smoothIRSolid

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (IRdata, InvTimes, segments, Sfluid, Rfluid, TEScale = 100, 
    dataScale = 1000, method = c("NLR", "QL"), sigma = NULL, 
    L = 1, maxR2star = 50, varest = c("RSS", "data"), verbose = TRUE, 
    lower = c(0, 0, 0), upper = c(0.95, 2, 2)) 
{
    mask <- segments > 1
    nvoxel <- sum(mask)
    ntimes <- length(InvTimes)
    InvTimes[InvTimes == Inf] <- 50 * max(InvTimes[InvTimes != 
        Inf])
    dimdata <- dim(IRdata)
    if (dimdata[1] != ntimes) 
        stop("estimateIRsolid: incompatible length of InvTimes")
    if (any(dimdata[-1] != dim(mask))) 
        stop("estimateIRsolid: incompatible dimension of segments")
    InvTimesScaled <- InvTimes/TEScale
    npar <- 3
    fx <- Rx <- Sx <- rsdx <- array(0, dim(mask))
    ICovx <- array(0, c(3, 3, prod(dim(mask))))
    Convx <- array(0, dim(mask))
    fx[segments == 1] <- 1
    Rx[segments == 1] <- Rfluid
    Sx[segments == 1] <- Sfluid
    Convx[segments == 1] <- 1
    ICovx[1, 1, segments == 1] <- 1e+20
    ICovx[2, 2, segments == 1] <- 1e+20
    ICovx[3, 3, segments == 1] <- 1e+20
    isConv <- array(FALSE, nvoxel)
    isThresh <- array(FALSE, nvoxel)
    modelCoeff <- array(0, c(npar, nvoxel))
    invCov <- array(0, c(npar, npar, nvoxel))
    rsigma <- array(0, nvoxel)
    if (method[1] == "QL") {
        if (is.null(sigma)) {
            method <- "NLR"
            warning("estimateIRsolid: method QL needs sigma estimated from fluid or supplied")
        }
        sig <- sigma/dataScale
        CL <- sig * sqrt(pi/2) * gamma(L + 0.5)/gamma(L)/gamma(1.5)
    }
    dim(IRdata) <- c(dimdata[1], prod(dim(segments)))
    IRdataSolid <- IRdata[, mask]
    thetas <- matrix(0, 3, nvoxel)
    thetas[3, ] <- IRdataSolid[(1:ntimes)[InvTimes == max(InvTimes)][1], 
        ]/dataScale
    thetas[2, ] <- 1/median(InvTimesScaled)
    thetas[1, ] <- 0.3
    if (verbose) {
        cat("Start estimation in", nvoxel, "voxel at", format(Sys.time()), 
            "\n")
        pb <- txtProgressBar(0, nvoxel, style = 3)
    }
    for (xyz in 1:nvoxel) {
        ivec <- IRdataSolid[, xyz]/dataScale
        th <- thetas[, xyz]
        res <- if (method[1] == "NLR") 
            try(nls(ivec ~ IRmix2(par, ITS, Sfluid, Rfluid), 
                data = list(ITS = InvTimesScaled, Sfluid = Sfluid, 
                  Rfluid = Rfluid), start = list(par = th), control = list(maxiter = 200, 
                  warnOnly = TRUE)), silent = TRUE)
        else try(nls(ivec ~ IRmix2QL(par, ITS, Sfluid, Rfluid, 
            CL, sig, L), data = list(ITS = InvTimesScaled, Sfluid = Sfluid, 
            Rfluid = Rfluid, CL = CL, sig = sig, L = L), start = list(par = th), 
            control = list(maxiter = 200, warnOnly = TRUE)), 
            silent = TRUE)
        if (class(res) != "try-error") {
            thhat <- coef(res)
            outofrange <- any(thhat != pmin(upper, pmax(lower, 
                thhat)))
        }
        if (class(res) == "try-error" || outofrange) {
            th <- pmin(upper, pmax(lower, th))
            res <- if (method[1] == "NLR") 
                try(nls(ivec ~ IRmix2(par, ITS, Sfluid, Rfluid), 
                  data = list(ITS = InvTimesScaled, Sfluid = Sfluid, 
                    Rfluid = Rfluid), start = list(par = th), 
                  algorithm = "port", control = list(maxiter = 200, 
                    warnOnly = TRUE), lower = lower, upper = upper), 
                  silent = TRUE)
            else try(nls(ivec ~ IRmix2QL(par, ITS, Sfluid, Rfluid, 
                CL, sig, L), data = list(ITS = InvTimesScaled, 
                Sfluid = Sfluid, Rfluid = Rfluid, CL = CL, sig = sig, 
                L = L), start = list(par = th), algorithm = "port", 
                control = list(maxiter = 200, warnOnly = TRUE), 
                lower = lower, upper = upper), silent = TRUE)
        }
        if (class(res) != "try-error") {
            sres <- if (varest[1] == "RSS") 
                getnlspars(res)
            else getnlspars2(res, shat[, xyz], sind)
            isConv[xyz] <- as.integer(res$convInfo$isConv)
            modelCoeff[, xyz] <- sres$coefficients
            if (sres$sigma != 0) {
                invCov[, , xyz] <- sres$invCov
                rsigma[xyz] <- sres$sigma
            }
        }
        if (verbose) 
            if (xyz%/%1000 * 1000 == xyz) 
                setTxtProgressBar(pb, xyz)
    }
    if (verbose) {
        close(pb)
        cat("Finished estimation", format(Sys.time()), "\n")
    }
    fx[mask] <- modelCoeff[1, ]
    Rx[mask] <- modelCoeff[2, ]
    Sx[mask] <- modelCoeff[3, ]
    ICovx[, , mask] <- invCov
    dim(ICovx) <- c(3, 3, dim(mask))
    Convx[mask] <- isConv
    rsdx[mask] <- rsigma
    list(fx = fx, Rx = Rx, Sx = Sx, Sf = Sfluid, Rf = Rfluid, 
        ICovx = ICovx, Convx = Convx, sigma = sigma, rsdx = rsdx)
  }

qMRI documentation built on Sept. 18, 2023, 9:08 a.m.