#' Euler
#'
#' Calculates the euler characteristic at a threshold level
#'
#' @param u Statistical value (typically the maxima of a cluster or statistical field).
#' @param df Degrees of freedom expressed as c(degrees of interest, degrees of error).
#' @param fieldType:
#' \itemize{
#' \item{'T'}{T-field}
#' \item{'F'}{F-field}
#' \item{'X'}{Chi-square field'}
#' \item{'Z'}{Gaussian field}
#' }
#' @return A vector of estimated euler characteristics for dimensions 0:D.
#'
#' @references
#' Worlsey K.J., (1996) A Unified Statistical Approach for Determining Significant Signals in Images of Cerebral Activation.
#' @author Zachary P. Christensen
#'
#' @seealso \code{\link{rftPval}}, \code{\link{resels}}
#' @examples
#'
#' ## generate some data as if we just fitted a linear regression
#' outimg1 <- makeImage(c(10,10,10), rt(1000))
#' maskimg <- getMask(outimg1)
#' fwhm <- estSmooth(outimg1, maskimg)
#' resels <- resels(maskimg, fwhm$fwhm)
#' ec <- euler(max(outimg1), c(1,10), fieldType='T')
#' pvox <- sum(ec*resels)
#'
#' @export euler
euler <- function(u, df, fieldType) {
if (missing(fieldType))
stop("Must specify fieldType") else if (missing(df))
stop("Must specify df") else if (missing(u))
stop("Must specify u")
ec <- c(0, 0, 0, 0)
if (fieldType == "T") {
ec[1] <- 1 - pt(u, df[2])
ec[2] <- (((4 * log(2))^(1/2))/(2 * pi)) * ((1 + ((u^2)/df[2]))^(-1/2 * (df[2] - 1)))
ec[3] <- (4 * log(2))/((2 * pi)^(3/2)) * ((1 + u^2/df[2])^((1 - df[2])/2)) * u/((df[2]/2)^(1/2)) *
exp(lgamma((df[2] + 1)/2) - lgamma(df[2]/2))
ec[4] <- (((4 * log(2))^(3/2))/((2 * pi)^2)) * ((1 + ((u^2)/df[2]))^(-1/2 * (df[2] - 1))) *
((((df[2] - 1)/df[2]) * (u^2)) - 1)
} else if (fieldType == "F") {
ec[1] <- 1 - pf(u, df[1], df[2])
ec[2] <- ((4 * log(2))/(2 * pi))^(1/2) * exp(lgamma((df[2] + df[1] - 1)/2) - (lgamma(df[2]/2) +
lgamma(df[1]/2))) * 2^(1/2) * (df[1] * u/df[2])^(1/2 * (df[1] - 1)) * (1 + df[1] *
u/df[2])^(-1/2 * (df[2] + df[1] - 2))
ec[3] <- ((4 * log(2))/(2 * pi)) * exp(lgamma((df[2] + df[1] - 2)/2) - (lgamma(df[2]/2) +
lgamma(df[1]/2))) * (df[1] * u/df[2])^(1/2 * (df[1] - 2)) * (1 + df[1] * u/df[2])^(-1/2 *
(df[2] + df[1] - 2)) * ((df[2] - 1) * df[1] * u/df[2] - (df[1] - 1))
ec[4] <- ((4 * log(2))/(2 * pi))^(3/2) * exp(lgamma((df[2] + df[1] - 3)/2) - (lgamma(df[2]/2) +
lgamma(df[1]/2))) * 2^(-1/2) * (df[1] * u/df[2])^(1/2 * (df[1] - 3)) * (1 + df[1] *
u/df[2])^(-1/2 * (df[2] + df[1] - 2)) * ((df[2] - 1) * (df[2] - 2) * (df[1] * u/df[2])^2 -
(2 * df[2] * df[1] - df[2] - df[1] - 1) * (df[1] * u/df[2]) + (df[1] - 1) * (df[1] -
2))
} else if (fieldType == "X") {
ec[1] <- 1 - pchisq(u, df[2])
ec[2] <- ((4 * log(2))/(2 * pi))^(1/2) * (u^(1/2 * (df[2] - 1)) * exp(-u/2 - lgamma(df[2]/2))/2^((df[2] -
2)/2))
ec[3] <- ((4 * log(2))/(2 * pi)) * (u^(1/2 * (df[2] - 1)) * exp(-u/2 - lgamma(df[2]/2))/2^((df[2] -
2)/2)) * (u - (df[2] - 1))
ec[4] <- ((4 * log(2))/(2 * pi))^(3/2) * (u^(1/2 * (df[2] - 1)) * exp(-u/2 - lgamma(df[2]/2))/2^((df[2] -
2)/2)) * (u^2 - (2 * df[2] - 1) * u + (df[2] - 1) * (df[2] - 2))
} else if (fieldType == "Z") {
ec[1] <- 1 - pnorm(u, df[2])
ec[2] <- (4 * log(2))^(1/2)/(2 * pi) * exp(-u^2/2)
ec[3] <- (4 * log(2))/((2 * pi)^(3/2)) * exp(-u^2/2) * u
ec[4] <- (4 * log(2))^(3/2)/((2 * pi)^2) * exp(-u^2/2) * (u^2 - 1)
}
ec
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.