R/dcemri_spline.R

In dcemriS4: A Package for Medical Image Analysis (S4 implementation)

##
##
## Copyright (c) 2009,2010 Brandon Whitcher and Volker Schmid
## All rights reserved.
## 
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are
## met:
## 
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer. 
##     * Redistributions in binary form must reproduce the above
##       copyright notice, this list of conditions and the following
##       disclaimer in the documentation and/or other materials provided
##       with the distribution.
##     * The names of the authors may not be used to endorse or promote
##       products derived from this software without specific prior
##       written permission.
## 
## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
## "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
## LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
## A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
## HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
## DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
## THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
## 
## $Id: dcemri_spline.R 332 2010-01-29 16:54:07Z bjw34032 $
##

#############################################################################
## setGeneric("dcemri.spline")
#############################################################################

setGeneric("dcemri.spline",
           function(conc, ...) standardGeneric("dcemri.spline"))
setMethod("dcemri.spline", signature(conc="array"),
          function(conc, time, img.mask, time.input=time,
                   model="weinmann", aif="tofts.kermode",
                   user=NULL, aif.observed=NULL, nriters=500,
                   thin=5, burnin=100, ab.hyper=c(1e-5,1e-5),
                   ab.tauepsilon=c(1,1/1000), k=4, p=25, rw=2,
                   knots=NULL, nlr=FALSE, t0.compute=FALSE,
                   samples=FALSE, multicore=FALSE, verbose=FALSE, 
		   response=FALSE, fitted=FALSE, ...)
          .dcemriWrapper("dcemri.spline", conc, time, img.mask, time.input,
                         model, aif, user, aif.observed, nriters, thin,
                         burnin, ab.hyper, ab.tauepsilon, k, p, rw, knots,
                         nlr, t0.compute, samples, multicore, verbose,
                         response, fitted, ...))

.dcemri.spline.single <- function(conc, time, D, time.input, p, rw, knots,
                                  k, A, t0.compute=FALSE, nlr=FALSE,
                                  nriters=500, thin=5, burnin=100,
                                  ab.hyper=c(1e-5,1e-5),
                                  ab.tauepsilon=c(1,1/1000), silent=0,
                                  multicore=FALSE, model=NULL,
                                  model.func=NULL, model.guess=NULL,
                                  samples=FALSE, B=NULL) {

  require("minpack.lm")
  ## Sanity check: conc must not contain any missing values
  if (any(is.na(conc))) {
    stop("Concentration time curves must not contain missing values.")
  }

  samplesize <- floor(nriters/thin)

  T <- length(time)
  TT <- T*T
  tau <- matrix(1000, p-rw, samplesize) # array(1000, c(p-rw,samplesize))
  beta <- matrix(0, p, samplesize) # array(0, c(p,samplesize))
  ##MAX <- array(0, c(samplesize))
  tauepsilon <- rep(1000, samplesize) # array(1000, c(samplesize))
  burnin <- min(burnin, samplesize)
  
  result <- .C("dce_spline_run",
               as.integer(1),
               as.integer(burnin),
               as.integer(c(1,1,1,T)),
               as.double(conc),
               as.double(tau),
               as.double(tauepsilon),
               as.double(D),
               as.integer(rw),
               as.double(beta),
               as.double(c(ab.hyper[1], ab.hyper[2], ab.tauepsilon[1],
                           ab.tauepsilon[2])),
               as.integer(p),
               as.double(1:T),
               as.double(1:T),
               as.double(1:T),
               as.double(1:TT),
               as.double(1:TT),
               as.double(1:TT),
               as.double(1:TT),
               as.double(t(D)),
               as.integer(silent),
               PACKAGE="dcemriS4")

  tau <- matrix(result[[5]], p-rw, samplesize)
  tauepsilon <- result[[6]]
  beta <- matrix(result[[9]], p, samplesize)

  ##
  ## Why is this run twice?!?
  ##
  
  result <- .C("dce_spline_run",
               as.integer(samplesize),
               as.integer(thin),
               as.integer(c(1,1,1,T)),
               as.double(conc),
               as.double(tau),
               as.double(tauepsilon),
               as.double(D),
               as.integer(rw),
               as.double(beta),
               as.double(c(ab.hyper[1], ab.hyper[2], ab.tauepsilon[1],
                           ab.tauepsilon[2])),
               as.integer(p),
               as.double(1:T),
               as.double(1:T),
               as.double(1:T),
               as.double(1:TT),
               as.double(1:TT),
               as.double(1:TT),
               as.double(1:TT),
               as.double(t(D)),
               as.integer(silent),
               PACKAGE="dcemriS4")

  tau <- matrix(result[[5]], p-rw, samplesize)
  tauepsilon <- result[[6]]
  beta <- matrix(result[[9]], p, samplesize)

  t0 <- 0
  if (t0.compute) {
    d <- matrix(NA, T, samplesize)
    for (j in 1:samplesize) {
      d[,j] <- D %*% beta[,j]
    }
    d1 <- apply(d, 1, quantile, probs=c(0.005), na.rm=TRUE)
    d2 <- apply(d, 1, median, na.rm=TRUE)
    du <- min(which(d1 > 0))
    ## beta.abl <- beta.abl2 <- rep(0, p)
    B2 <- splineDesign(knots, time.input, k-2)
    B2 <- B2[,(1:p)+1]

    for (j in 1:samplesize) {
      beta.abl <- 1:p
      for (q in 1:(p - 1)) {
	beta.abl[q] <- (beta[q+1, j] - beta[q, j]) *
          k / (knots[q+k+1] - knots[q+1])
      }
      beta.abl[p] <- 0
      ABL2 <- A %*% B2 %*% beta.abl
      du2 <- time[du] - d2[du] / ABL2[du]  
      t0[j] <- du2
    }
    if (sum(!is.na(t0)) == 0) {
      t0 <- 0
    } else {
      t0 <- median(t0, na.rm=TRUE)
    }
    if (t0 < min(time)) {
      t0 <- 0
    }
    if (t0 > max(time)) {
      t0 <- 0
    }
  }

  fitted <- vector("list", samplesize) # list()
  for (i in 1:samplesize) {
    fitted[[i]] <- B %*% beta[,i]
  }

  ##if (multicore && require("parallel")) {
  ##  MAX <- unlist(mclapply(fitted, max))
  ##} else {
    MAX <- unlist(lapply(fitted, max))
  ##}
  ##MAX <- rep(NA, samplesize)
  ##for (i in 1:samplesize) {
  ##  MAX[i] <- MAX0[[i]]
  ##}
  parameters <- list()
  if (nlr) {
    if (model=="AATH") {
      model.guess[2] <- median(MAX, na.rm=TRUE)
    }
    fcn <- function(p, time, x, N.Err, fcall, jcall)
      (x - do.call("fcall", c(list(time=time), as.list(p))))

    nls.lm.single <- function(fitted, par, fn, fcall, model, time) {
      fcall2 <- fcall
      if (length(fcall) > 1) {
	fcall2 <- fcall[[1]]
      }
      fit <- nls.lm(par=par, fn=fn, fcall=fcall2, time=time, x=fitted,
                    N.Err=sqrt(300), control=list(nprint=0, ftol=10^-20))
      if (model=="AATH" && fit$par$TC < 1e-6) {
	fit <- nls.lm(par=par[-3], fn=fn, fcall=fcall[[2]], time=time,
                      x=fitted, N.Err=sqrt(300),
                      control=list(nprint=0, ftol=10^-20))
	fit$par$TC <- 0
      }
      return(fit)
    }
    
    if (multicore && require("parallel")) {
      out <- mclapply(fitted, nls.lm.single, par=model.guess,
                           fn=fcn, fcall=model.func, model=model,
                           time=time-t0)
    } else {
      out <- lapply(fitted, nls.lm.single, par=model.guess,
                         fn=fcn, fcall=model.func, model=model,
                         time=time-t0)
    }

    if (model=="AATH") {
      E <- F <- TC <- ve <- rep(NA, samplesize) # NULL
      for (i in 1:samplesize) {
	E[i] <- exp(out[[i]]$par$logE)
	F[i] <- exp(out[[i]]$par$logF)
	TC[i] <- out[[i]]$par$TC
	ve[i] <- exp(out[[i]]$par$logve)
      }
      parameters <- list("E"=median(E, na.rm=TRUE),
                         "F"=median(F, na.rm=TRUE),
                         "TC"=median(TC, na.rm=TRUE),
                         "ve"=median(ve, na.rm=TRUE))
      if (samples) {
	parameters$E.samples <- E
        parameters$F.samples <- F
        parameters$TC.samples <- TC
        parameters$ve.samples <- ve
      }
    }
    if (model == "weinmann") {
      ktrans <- kep <- rep(NA, samplesize) # NULL
      for (i in 1:samplesize) {
	ktrans[i] <- exp(out[[i]]$par$logktrans)
	kep[i] <- exp(out[[i]]$par$logkep)
      }
      parameters <- list(ktrans=median(ktrans, na.rm=TRUE),
                         kep=median(kep, na.rm=TRUE))
      if (samples) {
	parameters$ktrans.sample <- ktrans
        parameters$kep.sample <- kep
      }
    }
  }

  list(beta=beta, tau=tau, tauepsilon=tauepsilon, t0=t0,
       Fp=median(MAX, na.rm=TRUE), Fp.samples=MAX, fitted=fitted,
       par=parameters)
}

.dcemri.spline <- function(conc, time, img.mask, time.input=time,
                          model="weinmann", aif="tofts.kermode",
                          user=NULL, aif.observed=NULL, nriters=500,
                          thin=5, burnin=100, ab.hyper=c(1e-5,1e-5),
                          ab.tauepsilon=c(1,1/1000), k=4, p=25, rw=2,
                          knots=NULL, nlr=FALSE, t0.compute=FALSE,
                          samples=FALSE, multicore=FALSE, verbose=FALSE,
                          response=FALSE, fitted=FALSE, ...) {

  ## dcemri.spline - a function for fitting Bayesian Penalty Splines to 
  ## DCE-MRI images and computing kinetic parameters
  ##
  ## authors: Volker Schmid, Brandon Whitcher
  ##
  ## input:
  ##        conc: array of Gd concentration,
  ##        time: timepoints of aquisition,
  ##        img.mask: array of voxels to fit,
  ##        D(=0.1): Gd dose in mmol/kg,
  ##        model: AIF... "weinman" or "parker",
  ##
  ## output: list with ktrans, kep, ve, std.error of ktrans and kep
  ##         (ktranserror and keperror)
  ##

  require("splines")
  if (nlr) {
    require("minpack.lm")
  }

  ## function to make precision matrix for random walk
  R <- function(taux, rw) {
    RR <- matrix(0, length(taux)+rw, length(taux)+rw)
    if (rw == 0) {
      for (i in 1:length(taux)) {
	RR[i,i] <- taux[i]
      }
    }
    if (rw == 1) {
      for (i in 1:length(taux)) {
	RR[i,i] <- RR[i,i] + taux[i]
	RR[i+1,i+1] <- RR[i+1,i+1] + taux[i]
	RR[i+1,i] <- RR[i+1,i] - taux[i]
	RR[i,i+1] <- RR[i,i+1] - taux[i]
      }
    }
    if (rw == 2) {
      for (i in 1:length(taux)) {
	RR[i,i] <- RR[i,i] + taux[i]
	RR[i+1,i+1] <- RR[i+1,i+1] + 4*taux[i]
	RR[i+2,i+2] <- RR[i+2,i+2] + taux[i]
	RR[i+1,i] <- RR[i+1,i] - 2*taux[i]
	RR[i,i+1] <- RR[i,i+1] - 2*taux[i]
	RR[i+2,i+1] <- RR[i+2,i+1] - 2*taux[i]
	RR[i+1,i+2] <- RR[i+1,i+2] - 2*taux[i]
	RR[i+2,i] <- RR[i+2,i] + taux[i]
	RR[i,i+2] <- RR [i,i+2] + taux[i]
      }
    }
    return(RR)
  }

  ## Main function

  knotpoints <- p
  if (is.null(knots)) {
    a <- min(time) - 1e-3 - (k-1) * (max(time)-min(time)) / (knotpoints-k+1)
    b <- 1e-3 + max(time) + k * (max(time)-min(time)) / (knotpoints-k+1)
    knots <- seq(a, b, (2e-3 + max(time) - min(time)) / (knotpoints - k + 1))
  }
  mod <- model
  nvoxels <- sum(img.mask)
  I <- nrow(conc)
  J <- ncol(conc)
  K <- nsli(conc)
  T <- length(time)
  Ti <- length(time.input)
  
  if (!is.numeric(dim(conc))) {
    I <- J <- K <- 1
  } else {
    if (length(dim(conc)) == 2) {
      J <- K <- 1
    } else {
      if (length(dim(conc)) == 3) {
        K <- 1
      }
    }
  }
  
  if (verbose) {
    cat("  Deconstructing data...", fill=TRUE)
  }
  img.mask <- ifelse(img.mask > 0, TRUE, FALSE)
  conc.mat <- matrix(conc[img.mask], nvoxels)
  conc.mat[is.na(conc.mat)] <- 0
  conc.list <- vector("list", nvoxels) # list()
  for (i in 1:nvoxels) {
    conc.list[[i]] <- conc.mat[i,]
  }

  switch(aif,
    tofts.kermode = {
      D <- 0.1; a1 <- 3.99; a2 <- 4.78; m1 <- 0.144; m2 <- 0.0111
      input <- D*(a1*exp(-m1*time)+a2*exp(-m2*time))
    },
    fritz.hansen = {
      D <- 1; a1 <- 2.4; a2 <- 0.62; m1 <- 3.0; m2 <- 0.016
      input <- D*(a1*exp(-m1*time)+a2*exp(-m2*time))	 
    },
    observed = {
      input <- aif.observed
    },
    print("WARNING: AIF parameters must be specified!"))

  model.func <- model.guess <- NULL

  if (nlr && model != "weinmann" && model != "AATH") {
    stop("Only model=\"weinmann\" or model=\"AATH\" acceptable AIFs for nlr=TRUE")
  }

  if (model == "weinmann") {
    ktrans <- kep <- list(par=rep(NA, nvoxels), error=rep(NA, nvoxels))
    sigma2 <- rep(NA, nvoxels)
    model.func <- function(time, logktrans, logkep) {
      ktrans <- exp(logktrans)
      kep <- exp(logkep)
      erg <- ktrans * exp(-kep*(time))
      eval(erg)
    }
    model.guess <- list(logktrans=-1, logkep=0)
  }

  if (model=="AATH") {
    E <- F <- TC <- ve <- list(par=rep(NA, nvoxels), error=rep(NA, nvoxels))
    sigma2 <- rep(NA, nvoxels)

    model.func <- list()
    model.func[[1]] <- function(time, logF, logE, TC, logve) {
      TC2 <- 2*TC
      if (TC < 0) {
        TC2 <- 0
      }
      E <- exp(logE)
      F <- exp(logF)
      ve <- exp(logve)      
      kep <- E * F/ve
      erg <- E * exp(-kep*(time-TC))
      erg[time<TC] <- 1 # - time[time<TC2]*(1-E) / TC
      erg <- erg*F
      if (TC < 0) {
	erg <- rep(-10^16, T)
      }
      eval(erg)
    }
    model.func[[2]] <- function(time, logF, logE, logve) {
      E <- exp(logE)
      F <- exp(logF)
      ve <- exp(logve)      
      kep <- E * F / ve
      erg <- E * exp(-kep * time)
      erg <- erg * F
      eval(erg)
    }
    model.guess <- list(logE=log(0.6), logF=log(2), TC=0, logve=log(0.05))
  }

  ## define A and B

  p <- length(knots) - k
  B <- splineDesign(knots, time.input, k, outer.ok=TRUE)
  if (sum(B[, dim(B)[2]] == 0) == dim(B)[1]) {
    B <- B[,-dim(B)[2]]
  }
  if (sum(B[,1] == 0) == dim(B)[1]) {
    B <- B[,-1]
  }
  p <- dim(B)[2]
  A <- matrix(0, T, Ti)
  ni <- time
  for (i in 1:T) {
    for (j in 1:Ti) {
      if (time.input[j] <= time[i]) {
	ni[i] <- j
      }
    }
  }
  for (i in 1:T) {
    for (j in 1:ni[i]) {
      A[i,j] <- input[1+ni[i]-j]
    }
  }
  A <- A * mean(diff(time.input), na.rm=TRUE)
  A[is.na(A)] <- 0
  D <- A %*% B
  ## T <- length(time)
  
  if (verbose) {
    cat("  Estimating the parameters...", fill=TRUE)
  }

  if (multicore && require("parallel")) {
    fit <- mclapply(conc.list, FUN=.dcemri.spline.single, time=time,
                    D=D, time.input=time.input, p=p, rw=rw, knots=knots,
                    k=k, A=A, nriters=nriters, thin=thin, burnin=burnin,
                    ab.hyper=ab.hyper, ab.tauepsilon=ab.tauepsilon,
                    t0.compute=t0.compute, nlr=nlr, multicore=TRUE,
                    model=model, model.func=model.func,
                    model.guess=model.guess, samples=samples, B=B)
  } else {
    fit <- lapply(conc.list, FUN=.dcemri.spline.single, time=time, D=D,
                  time.input=time.input, p=p, rw=rw, knots=knots, k=k,
                  A=A, nriters=nriters, thin=thin, burnin=burnin,
                  ab.hyper=ab.hyper, ab.tauepsilon=ab.tauepsilon,
                  t0.compute=t0.compute, nlr=nlr, multicore=FALSE,
                  model=model, model.func=model.func,
                  model.guess=model.guess, samples=samples, B=B)
  }
  
  if (verbose) {
    cat("  Reconstructing results...", fill=TRUE)
  }

  samplesize <- floor(nriters/thin)

  t0 <- numeric(nvoxels)
  for (k in 1:nvoxels) {
    t0 <- fit[[k]]$t0
  }
  t0.img <- array(NA, c(I,J,K))
  t0.img[img.mask] <- t0
  t0 <- t0.img
  rm(t0.img)

  Fp <- numeric(nvoxels)
  for (k in 1:nvoxels) {
    Fp[k] <- fit[[k]]$Fp
  }
  Fp.img <- array(NA, c(I,J,K))
  Fp.img[img.mask] <- Fp
  Fp.samples <- matrix(NA, nvoxels, samplesize)
  for (i in 1:nvoxels) {
    Fp.samples[i,] <- fit[[k]]$Fp.samples
  }
  if (samples) {
    Fp <- array(NA, c(I,J,K,samplesize))
    for (j in 1:samplesize) {
      Fp[,,,j][img.mask] <- Fp.samples[,j]
    }
  }

  if (nlr) {
    ktrans <- ve <- rep(NA, nvoxels) # NULL
    if (model=="weinmann") {
      kep <- NULL
      if (samples) {
	ktrans.samples <- matrix(NA, nvoxels, samplesize)
      }
      for (k in 1:nvoxels) {
	ktrans[k] <- fit[[k]]$par$ktrans
      	if (samples) {
	  ktrans.samples[k,] <- fit[[k]]$par$ktrans.samples
        }
      }
      if (samples) {
	kep.samples <- array(NA, c(nvoxels,samplesize))
      }
      for (k in 1:nvoxels) {
	kep[k] <- fit[[k]]$par$kep
	if (samples) {
	  kep.samples[k,] <- fit[[k]]$par$kep.samples
        }
      }
      ve <- ktrans / kep
      if (samples) {
	ve.samples <- ktrans.samples / kep.samples
      }
    }
    if (model == "AATH") {
      E <- F <- TC <- rep(NA, nvoxels)
      if (samples) {
	E.samples <- matrix(NA, nvoxels, samplesize)
      }
      for (k in 1:nvoxels) {
	E[k] <- fit[[k]]$par$E
	if (samples) {
	  E.samples[k,] <- fit[[k]]$par$E.samples
        }
      }
      if (samples) {
	F.sample <- matrix(NA, nvoxels, samplesize)
      }
      for (k in 1:nvoxels) {
	F[k] <- fit[[k]]$par$F
	if (samples) {
	  F.samples[k,] <- fit[[k]]$par$F.samples
        }
      }
      if (samples) {
	TC.samples <- matrix(NA, nvoxels, samplesize)
      }
      for (k in 1:nvoxels) {
	TC[k] <- fit[[k]]$par$TC
	if (samples) {
	  TC.samples[k,] <- fit[[k]]$par$TC.samples
        }
      }
      if (samples) {
	ve.samples <- matrix(NA, nvoxels, samplesize)
      }
      for (k in 1:nvoxels) {
	ve[k] <- fit[[k]]$par$ve
	if (samples) {
	  ve[k,] <- fit[[k]]$par$ve.samples
        }
      }
      ktrans <- E*F
      if (samples) {
	ktrans.samples <- E.samples / F.samples
      }
    }
  }

  beta.sample <- array(NA, c(nvoxels,p,samplesize))
  response.sample <- array(NA, c(nvoxels,Ti,samplesize))
  fitted.sample <- array(NA, c(nvoxels,T,samplesize))
  for (k in 1:nvoxels) {
    beta.sample[k,,] <- fit[[k]]$beta
    for (j in 1:samplesize) {
      response.sample[k,,j] <- fit[[k]]$fitted[[j]]
      fitted.sample[k,,j] <- D %*% beta.sample[k,,j]
    }
  }
  
  if ((response || fitted) && samples) {
    beta <- array(NA, c(I,J,K,p,samplesize))     # a 5-dimensional array!!!
  }
  if (response || fitted) {
    beta.med <- array(NA, c(I,J,K,p))
  }
  if (fitted && samples) {
    fitted <- array(NA, c(I,J,K,T,samplesize))   # a 5-dimensional array!!!
  }
  if (fitted) {
    fitted.med <- array(NA, c(I,J,K,T))
  }
  if (response && samples) {
    response <- array(NA, c(I,J,K,T,samplesize)) # a 5-dimensional array!!!
  }
  if (response) {
    response.med <- array(NA, c(I,J,K,T))
  }
  ktrans.med <- ve.med <- array(NA, c(I,J,K))
  ktrans.med[img.mask] <- ktrans
  ve.med[img.mask] <- ve
  if (model=="weinmann") {
    kep.med <- array(NA, c(I,J,K))
    kep.med[img.mask] <- kep
  }
  if (model=="AATH") {
    E.med <- F.med <- TC.med <- array(NA, c(I,J,K))
    E.med[img.mask] <- E
    F.med[img.mask] <- F
    TC.med[img.mask] <- TC
  } 
  if (samples) {
    ktrans <- ve <- array(NA, c(I,J,K,samplesize))
    for (i in 1:samplesize) {
      ktrans[,,,i][img.mask] <- ktrans.samples[,i]
      ve[,,,i][img.mask] <- ve.samples[,i]
    }
    if (model=="weinmann") {
      kep <- array(NA,c(I,J,K,samplesize))
      for (i in 1:samplesize) {
        kep[,,,i][img.mask] <- kep.samples[,i]
      }
    }
    if (model=="AATH") {
      E <- F <- TC <- array(NA, c(I,J,K,samplesize))
      for (i in 1:samplesize) {
        E[,,,i][img.mask] <- E.samples[,i]
        F[,,,i][img.mask] <- F.samples[,i]
        TC[,,,i][img.mask] <- TC.samples[,i]
      }
    }
  }
# }
  if (fitted || response) {
    if (I > 1) {
      for (j in 1:p) {
        beta.med[,,,j][img.mask] <- apply(beta.sample[,j,], 1, median)
        if (samples)
          for (i in 1:samplesize) {
            beta[,,,j,i][img.mask] <- beta.sample[,j,i]
          }
      }
    } else {
      for (j in 1:p) {
        beta.med[,,,j][img.mask] <- median(beta.sample[,j,])
        for (i in 1:samplesize) {
          beta[,,,j,i][img.mask] <- beta.sample[,j,i]
        }
      }
    }
  }

if (fitted || response) {
    if (I == 1)
    for (j in 1:T) {
      if (response)
        response.med[,,,j][img.mask] <- median(response.sample[,j,])
      if (fitted)
        fitted.med[,,,j][img.mask] <- median(fitted.sample[,j,])
      }
    if (I > 1) 
    for (j in 1:T) {
      if (response)
        response.med[,,,j][img.mask] <- apply(response.sample[,j,], 1, median)
      if (fitted)
        fitted.med[,,,j][img.mask] <- apply(fitted.sample[,j,], 1, median)
      }
    if (samples)
    for (i in 1:samplesize) 
    for (j in 1:T) {
      if (response)
        response[,,,j,i][img.mask] <- response.sample[,j,i]
      if (fitted)
        fitted[,,,j,i][img.mask] <- fitted.sample[,j,i]
      }
  }

  return.list <- vector("list")
  if (response || fitted) {
    return.list[["beta"]] <- beta.med
  }
  if ((response || fitted) && samples) {
    return.list[["beta.sample"]] <- beta
  }
  if (fitted) {
    return.list[["fit"]] <- fitted.med
  }
  if (fitted && samples) {
    return.list[["fit.sample"]] <- fitted
  }
  if (response) {
    return.list[["response"]] <- response.med
  }
  if (response && samples) {
    return.list[["response.sample"]] <- response
  }
  return.list[["Fp"]] <- Fp.img
  return.list[["A"]] <- A
  return.list[["B"]] <- B
  return.list[["D"]] <- D
  if (t0.compute) {
    return.list[["t0"]] <- t0
  }
  if (nlr) {
    if (model=="weinmann") {
      return.list[["kep"]] <- kep.med
    }
    if (model=="AATH") {
      return.list[["E"]] <- E.med
      return.list[["F"]] <- F.med
      return.list[["TC"]] <- TC.med
    }
    return.list[["ktrans"]] <- ktrans.med
    return.list[["ve"]] <- ve.med
    if (samples) {
      return.list[["ktrans.sample"]] <- ktrans
      return.list[["ve.sample"]] <- ve
      if (model=="weinmann") {
        return.list[["kep.samples"]] <- kep
      }
      if (model=="AATH") {
        return.list[["E.samples"]] <- E
        return.list[["F.samples"]] <- F
        return.list[["TC.samples"]] <- TC
      }
    }
  }
  if (samples) {
    return.list[["Fp.samples"]] <- Fp
  }

  return(return.list)
}