R/varest.r
In WIAS-BERLIN/aws: Adaptive Weights Smoothing
Defines functions
awsgsigmasd
awslinsd
estGlobalSigma
AFLocalSigma
iniSigmaNCchi
IQQdiff
awsLocalSigma
estimateSigmaCompl
Documented in
AFLocalSigma
awslinsd
awsLocalSigma
estGlobalSigma
estimateSigmaCompl
estimateSigmaCompl <- function(magnitude,phase,mask,kstar=20,kmin=8,hsig=5,lambda=12,verbose=TRUE){
## kmin = 10 corresponds to an initial bandwidth of 1.47 giving positive weight to direct neighbors and
## 2D diagonal neigbors
args <- sys.call(-1)
sdim <- dim(mask)
# if(!is.numeric(magnitude)){
# if (verbose) cat("reading Magnitude file ... ")
# R <- readNIfTI(magnitude, reorient = FALSE)
# } else {
R <- magnitude
# }
# if(!is.numeric(phase)){
# if (verbose) cat("reading Phase file ... ")
# Ph <- readNIfTI(phase, reorient = FALSE)
# } else {
Ph <- phase
# }
ComplImg <- array(0,c(2,sdim))
ComplImg[1,,,] <- R*cos(Ph)
ComplImg[2,,,] <- R*sin(Ph)
## find the number of usable cores
mc.cores <- setCores(, reprt = FALSE)
##
## start smoothing and variance estimation
##
n <- prod(sdim)
lambda0 <- 1e40
sigma2 <- array(1e10,sdim)
# just inilitialize with something large, first step is nonadaptive due to lambda0
k <- kmin
hmax <- 1.25^(kstar/3)
## preparations for median smoothing
parammd <- getparam3d(hsig,c(1,1))
nwmd <- length(parammd$w)
if (verbose) pb <- txtProgressBar(min = 0, max = kstar-kmin+1, style = 3)
bi <- array(1,sdim)
zobj <- list(theta=ComplImg, bi=bi)
if (verbose) {
mae <- NULL
protocol <- matrix("", kstar-kmin+1, 1, dimnames = list(paste("step", kmin:kstar), "protocol"))
}
while (k <= kstar) {
## determine the actual bandwidth for this step
hakt <- gethani(1, 1.25*hmax, 2, 1.25^k, c(1,1), 1e-4)
## we need the (approx.) size of the weigthing scheme array
dlw <- (2*trunc(hakt/c(1, 1, 1))+1)[1:3]
## perform the actual adaptive smoothing
zobj <- .Fortran(C_cplxawss,
as.double(ComplImg),
as.logical(mask),
as.integer(2),
as.integer(sdim[1]),
as.integer(sdim[2]),
as.integer(sdim[3]),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(sigma2),
bi = as.double(zobj$bi),
theta = double(2*n),
sigma2 = double(n),
as.integer(mc.cores),
double(prod(dlw)),
as.double(c(1,1)),
double(2 * mc.cores))[c("bi", "theta", "hakt","sigma2")]
##
## now get local median variance estimates
##
dim(zobj$sigma2) <- sdim
sigma2 <- .Fortran(C_mediansm,
as.double(zobj$sigma2),
as.logical(mask),
as.integer(sdim[1]),
as.integer(sdim[2]),
as.integer(sdim[3]),
as.integer(parammd$ind),
as.integer(nwmd),
double(nwmd*mc.cores), # work(nw,nthreds)
as.integer(mc.cores),
sigma2n = double(n))$sigma2n/0.6931
# sigma2n containes sum of 2 independent squared residuals
# 0.6931 approximates median \chi_2 /2
# needed to get correct results
## use maximum ni
bi <- zobj$bi <- pmax(bi, zobj$bi)
## some verbose stuff
if (verbose) {
protocol[k-kmin+1,1] <- paste("bandwidth: ", signif(hakt, 3),
"sigma: mean: ", signif(sqrt(mean(sigma2[mask])),3),
"median: ", signif(sqrt(median(sigma2[mask])),3),
"sd: ", signif(sd(sqrt(sigma2[mask])),3),
"median(bi):", signif(median(zobj$bi[mask]),3),
"max(bi):", signif(max(zobj$bi[mask]),3))
setTxtProgressBar(pb, k-kmin+1)
}
## go for next iteration
k <- k+1
lambda0 <- lambda
gc()
}
dim(zobj$theta) <- c(2,sdim)
# return estimated parameters of rician distribution
list(sigma=array(sqrt(sigma2),sdim),
theta=array(sqrt(zobj$theta[1,,,]^2+zobj$theta[2,,,]^2),sdim),
sigmal=array(sqrt(zobj$sigma2),sdim),mask=mask,
protocol=protocol,args=args)
}
awsLocalSigma <- function(y, steps, mask, ncoils, vext=c(1,1),
lambda = 5, minni=2, hsig=5, sigma=NULL, family=c("NCchi","Gauss"),
verbose = FALSE, trace = FALSE, u=NULL){
if(trace) tergs <- array(0,c(steps,4,sum(mask))) else tergs <- NULL
family <- match.arg(family[1],c("NCchi","Gauss","Gaussian"))
if(family == "NCchi"){
varstats <- sofmchi(ncoils)
th <- seq(0,30,.01)
z <- fncchiv(th,varstats)
minz <- min(z)
th <- th[z>min(z)]
minth <- min(th)
z <- z[z>min(z)]
nz <- nls(z~(a*th+b*th*th+c*th*th*th+d)/(a*th+b*th*th+c*th*th*th+1+d),data=list(th=th,z=z),start=list(a=1,b=1,c=1,d=1))
vpar <- c(minth,minz,coefficients(nz))
## this provides an excellent approximation for the variance reduction in case of low ncp
}
if(length(vext)==3) vext <- vext[2:3]/vext[1]
## dimension and size of cubus
ddim <- dim(y)
n <- prod(ddim)
## test dimension
if (length(ddim) != 3) stop("first argument should be a 3-dimentional array")
## check mask
if (is.null(mask)) mask <- array(TRUE, ddim)
if(length(mask) != n) stop("dimensions of data array and mask should coincide")
## initial value for sigma_0
if(is.null(sigma)){
# sigma <- sqrt( mean( y[mask]^2) / 2 / ncoils)
sigma <- IQQdiff( y, mask, .25, verbose=verbose)
if("NCchi" == family){
## iterative improvement for NC chi-distribution
sigma <- iniSigmaNCchi( y, mask, .25, ncoils, sigma)
sigma <- iniSigmaNCchi( y, mask, .25, ncoils, sigma)
sigma <- iniSigmaNCchi( y, mask, .25, ncoils, sigma)
}
}
##
## Prepare for diagnostics plots
##
if(verbose){
mslice <- (ddim[3]+1)/2
ymslice <- y[,,mslice]
ymslice[!mask[,,mslice]] <- 0
if(!is.null(u)&&"NCchi" == family){
par(mfrow=c(2,4),mar=c(3,3,3,1),mgp=c(2,1,0))
} else {
par(mfrow=c(2,3),mar=c(3,3,3,1),mgp=c(2,1,0))
}
} else {
cat("step")
}
## define initial arrays for parameter estimates and sum of weights (see PS)
th <- y
ksi <- array( y^2, ddim)
ni <- array( 1, ddim)
sigma <- array(sigma, ddim)
# initialize array for local sigma by global estimate
mc.cores <- setCores(,reprt=FALSE)
## preparations for median smoothing
parammd <- getparam3d(hsig, vext)
## iterate PS starting with bandwidth h0
if(family=="NCchi"){
# precompute values of log(besselI) for interpolation
nfb <- 200 * ncoils # for arguments > nfb asymp. approx can be used
x <- 1:nfb
flb <- x + log(besselI(x, ncoils-1, TRUE))
}
for (i in 1:steps) {
h <- 1.25^((i-1)/3)
param <- getparam3d(h, vext)
nw <- length(param$w)
## perform one step PS with bandwidth h
if(family=="NCchi"){
z <- .Fortran(C_awslchi2,
as.double(y),# data
as.double(ksi),# \sum w_ij S_j^2
ni = as.double(ni),
as.double(sigma),
as.double(vpar),# parameters for var. reduction
as.double(ncoils),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(param$ind),
as.double(param$w),
as.integer(nw),
as.double(minni),
double(nw*mc.cores), # wad(nw,nthreds)
double(nw*mc.cores), # sad(nw,nthreds)
as.double(lambda),
as.integer(mc.cores),
as.integer(floor(ncoils)),
double(floor(ncoils)*mc.cores), # work(L,nthreds)
th = double(n),
sigman = double(n),
ksi = double(n),
as.double(flb),
as.integer(nfb))[c("ni","ksi","th","sigman")]
thchi <- z$th
ksi <- z$ksi
thchi[!mask] <- 0
} else{
## assume Gaussian data
z <- .Fortran(C_awslgaus,
as.double(y),# data
as.double(th),# \sum w_ij S_j
ni = as.double(ni),
as.double(sigma),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(param$ind),
as.double(param$w),
as.integer(nw),
as.double(minni),
as.double(lambda),
th = double(n),
sigman = double(n))[c("ni","th","sigman")]
}
## extract sum of weigths (see PS) and consider only voxels with ni larger then mean
th <- array(z$th,ddim)
ni <- array(z$ni,ddim)
z$sigman[z$sigman==0] <- median(z$sigman[z$sigman>0])
nmask <- sum(mask)
if(verbose) cat("local estimation in step ",i," completed",format(Sys.time()),"\n")
##
## nonadaptive smoothing of estimated standard deviations
##
if(any(ni[mask]>minni)){
z$sigman[mask] <- z$sigman[mask]*(1+runif(nmask,-.0001,.0001))
##
## avoid ties in local neighborhood cubes of z$sigman
##
nwmd <- length(parammd$w)
sigma <- .Fortran(C_mediansm,
as.double(z$sigman),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(parammd$ind),
as.integer(nwmd),
double(nwmd*mc.cores), # work(nw,nthreds)
as.integer(mc.cores),
sigman = double(n))$sigman
}
dim(sigma) <- ddim
mask[sigma==0] <- FALSE
if(verbose) cat("local median smoother in step ",i," completed",format(Sys.time()),"\n")
##
## diagnostics
##
if(verbose){
meds <- median(sigma[mask])
means <- mean(sigma[mask])
image(ymslice,col=grey(0:255/255))
title(paste("S max=",signif(max(y[mask]),3)," median=",signif(median(y[mask]),3)))
image(th[,,mslice],col=grey(0:255/255))
title(paste("E(S) max=",signif(max(th[mask]),3)," median=",signif(median(th[mask]),3)))
image(sigma[,,mslice],col=grey(0:255/255),zlim=c(0,max(sigma[mask])))
title(paste("sigma max=",signif(max(sigma[mask]),3)," median=",signif(meds,3)))
image(ni[,,mslice],col=grey(0:255/255))
title(paste("Ni max=",signif(max(ni[mask]),3)," median=",signif(median(ni[mask]),3)))
plot(density(sigma[mask]),main="density of sigma")
plot(density(ni[mask]),main="density of Ni")
cat("mean sigma",means,"median sigma",meds,"sd sigma",sd(sigma[mask]),"\n")
if(!is.null(u)&&"NCchi" == family){
thchims <- fncchir(th/sigma,varstats)*sigma
thchims[!mask] <- u[!mask]
image(abs(thchims-u)[,,mslice],col=grey(0:255/255))
title("abs Error in thchims")
plot(density((thchims-u)[mask]),main="density of thchims-u")
cat("MAE(th)",mean(abs(thchi-u)[mask]),"RMSE(th)",sqrt(mean((thchi-u)[mask]^2)),"MAE(thms)",mean(abs(thchims-u)[mask]),"RMSE(thms)",sqrt(mean((thchims-u)[mask]^2)),"\n")
}
} else {
cat(" ",i)
}
if(trace) {
tergs[i,1,] <- ni[mask]
tergs[i,2,] <- th[mask]
tergs[i,3,] <- z$sigma[mask]
tergs[i,4,] <- sigma[mask]
}
}
## END PS iteration
if(!verbose) cat("\n")
if(family=="NCchi"){
thchi <- fncchir(th/sigma,varstats)*sigma
thchi[!mask] <- 0
} else {
thchi <- NULL
}
## this is the result (th is expectation, not the non-centrality parameter !!!)
invisible(list(sigma = sigma,
sigmal = array(z$sigman,ddim),
theta = th,
thchi = thchi,
ni = ni,
tergs = tergs,
mask = mask))
}
IQQ <- function (x, q = .25, na.rm = FALSE, type = 7){
cqz <- qnorm(.05)/qnorm(q)
x <- as.numeric(x)
z <- diff(quantile(x, c(q, 1-q), na.rm = na.rm, names = FALSE, type = type))
z0 <- 0
while(abs(z-z0)>1e-5*z0){
# outlier removal
z0 <- z
# cat(sum(x>(z0*cqz))," ")
x <- x[x<(z0*cqz)]
z <- diff(quantile(x, c(q, 1-q), na.rm = na.rm, names = FALSE, type = type))
}
z
}
IQQdiff <- function(y, mask, q = .25, verbose = FALSE) {
cq <- qnorm(1-q)*sqrt(2)*2
sx <- IQQ( diff(y[mask]), q)/cq
sy <- IQQ( diff(aperm(y,c(2,1,3))[aperm(mask,c(2,1,3))]), q)/cq
sz <- IQQ( diff(aperm(y,c(3,1,2))[aperm(mask,c(3,1,2))]), q)/cq
if(verbose) cat( "Pilot estimates of sigma", sx, sy, sz, "\n")
min( sx, sy, sz)
}
iniSigmaNCchi <- function(y, mask, q, L, sigma, verbose=FALSE){
#
# robust improvement for sigma using moment estimates and IQQdiff
#
meany <- mean(y[mask]^2)
eta <- sqrt(max( 0, meany/sigma^2-2*L))
m <- sqrt(pi/2)*gamma(L+1/2)/gamma(L)/gamma(3/2)*hg1f1(-1/2,L,-eta^2/2)
v <- max( .01, 2*L+eta^2-m^2)
if(verbose) cat( eta, m, v, "\n")
IQQdiff( y, mask, q)/sqrt(v)
}
AFLocalSigma <- function(y,ncoils,level=NULL,mask=NULL,h=2,hadj=1,vext = c( 1, 1)){
##
## estimate effective sigma and effective ncoils (L) according to Aja-Fernandez MRI 2013
##
# sh2B stands for \hat{\sigma}^2_{n,B}
# sh2L stands for \hat{\sigma}^2_{nL}
# sh2eB stands for \hat{\sigma}^2_{eff,B}
# sh2eS stands for \hat{\sigma}^2_{eff,S}
# LeB stands for \hat{L}_{eff,B}
# vmlb - Var(M_L(x)|x \in B)
# vmlbx - Var(M_L(x)|x \in B)_x
# phix - \hat{\Phi}_n(x)
# seff - 3D array of local variance estimates
# Leff - 3D array of effective L
ddim <- dim(y)
n <- prod(ddim)
if(is.null(level)&is.null(mask)){
warning("no background definition, need either level or mask")
return(invisible(NULL))
}
if(is.null(level)) indB <- !mask else indB <- y<level
sh2Lsimple <- mean(y[indB]^2)/2
mask1 <- array(TRUE,ddim)
##
## local variance estimates in vx
##
vx <- .Fortran(C_afmodevn,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask1),
as.double(h),
as.double(vext),
sigma = double(n))$sigma
dim(vx) <- ddim
vxb <- vx[indB]
vxs <- vx[!indB]
dv2b <- density( vxb[vxb>0], n = 4092, adjust = hadj, to = min( max(vxb[vxb>0]), median(vxb[vxb>0])*5) )
sh2B <- dv2b$x[dv2b$y == max(dv2b$y)][1]
dv2s <- density( vxs[vxs>0], n = 4092, adjust = hadj, to = min( max(vxs[vxs>0]), median(vxs[vxs>0])*5) )
sh2eS <- dv2s$x[dv2s$y == max(dv2s$y)][1]
m2 <- .Fortran(C_afmodem2,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask1),
as.double(h),
as.double(vext),
sm = double(n))$sm
dim(m2) <- ddim
m2b <- m2[indB]
dm2b <- density( m2b[m2b>0], n = 4092, adjust = hadj, to = min( max(m2b[m2b>0]), median(m2b[m2b>0])*5) )
sh2L <- dm2b$x[dm2b$y == max(dm2b$y)][1]/2
cat("sh2B",sh2B,"sh2Lsimple",sh2Lsimple,"sh2L",sh2L,"\n")
##
## now get sh2eB and LeB
##
LeBn <- ncoils
sh2eBn <- sh2B/2/(LeBn-gamma(LeBn+.5)^2/gamma(LeBn)^2)
LeB <- 1
sh2eB <- 0
while(abs(sh2eB-sh2eBn)>1e-5||abs(LeB-LeBn)>1e-5){
LeB <- LeBn
sh2eB <- sh2eBn
LeBn <- sh2L/sh2eB
sh2eBn <- sh2B/2/(LeBn-gamma(LeBn+.5)^2/gamma(LeBn)^2)
cat("LeB",LeB,"LeBn",LeBn,"sh2eB",sh2eB,"sh2eBn",sh2eBn,"\n")
}
sh2eB <- sh2eBn
LeB <- sh2L/sh2eB
phix <- pmax(0,pmin(1,sh2L/(m2-sh2L)))
dim(phix) <- dim(m2)
seff <- sqrt((1-phix)*sh2eS+phix*sh2eB)
Leff <- sh2L/seff^2
list(sigma=seff,Leff=Leff,sh2L=sh2L,sh2B=sh2B,sh2eS=sh2eS,sh2eB=sh2eB,LeB=LeB,m2=m2,vx=vx,indB=indB)
}
estGlobalSigma <- function(y, mask=NULL, ncoils=1, steps=16, vext=c(1,1),
lambda=20, hinit=2, hadj=1, q=.25, level=NULL,
sequence=FALSE,method=c("awsVar","awsMAD",
"AFmodevn","AFmodem1chi","AFbkm2chi","AFbkm1chi")){
##
## estimate global scale parameter sigma for NCchi distributed data
##
method <- match.arg(method)
qni <- .8
## methods using the propagation-separation (PS) approach
if(method %in% c("awsVar","awsMAD")){
varstats <- sofmchi(ncoils)
if(length(vext)==3) vext <- vext[2:3]/vext[1]
## dimension and size of cubus
ddim <- dim(y)
n <- prod(ddim)
## test dimension
if (length(ddim) != 3) stop("first argument should be a 3-dimentional array")
## check mask
if (is.null(mask)) mask <- array(TRUE, ddim)
if(length(mask) != n) stop("dimensions of data array and mask should coincide")
## initial value for sigma_0
# sigma <- sqrt( mean( y[mask]^2) / 2 / ncoils)
sigma <- IQQdiff( y, mask, q)
# cat( "sigmahat1", sigma, "\n")
sigma <- iniSigmaNCchi( y, mask, q, ncoils, sigma)
# cat( "sigmahat2", sigma, "\n")
sigma <- iniSigmaNCchi( y, mask, q, ncoils, sigma)
# cat( "sigmahat3", sigma, "\n")
sigma <- iniSigmaNCchi( y, mask, q, ncoils, sigma)
# cat( "sigmahat4", sigma,"\n")
## define initial arrays for parameter estimates and sum of weights (see PS)
th <- array( 1, ddim)
ni <- array( 1, ddim)
# y <- y/sigma # rescale to avoid passing sigma to awsvchi
minlev <- sqrt(2)*gamma(ncoils+.5)/gamma(ncoils)
if (sequence) sigmas <- lind <- minni <- numeric(steps)
mc.cores <- setCores(,reprt=FALSE)
## iterate PS starting with bandwidth hinit
for (i in 1:steps) {
h <- hinit * 1.25^((i-1)/3)
param <- getparam3d(h,vext)
nw <- length(param$w)
fncchi <- fncchiv(th/sigma,varstats)
## correction factor for variance of NC Chi distribution
## perform one step PS with bandwidth h
if(method=="awsVAR"){
z <- .Fortran(C_awsvchi,
as.double(y), # data
as.double(th), # previous estimates
ni = as.double(ni),
as.double(fncchi/2),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(param$ind),
as.double(param$w),
as.integer(nw),
as.double(lambda),
as.double(sigma),
th = double(n),
sy = double(n))[c("ni","th","sy")]
} else { # method=="awsMAD"
z <- .Fortran(C_awsadchi,
as.double(y), # y(n1,n2,n3)
as.double(th), # th(n1,n2,n3)
ni = as.double(ni), # ni(n1,n2,n3)
as.double(fncchi/2), # fns(n1,n2,n3)
as.integer(mask), # mask(n1,n2,n3)
as.integer(ddim[1]), # n1
as.integer(ddim[2]), # n2
as.integer(ddim[3]), # n3
as.integer(param$ind), # ind(3,nw)
as.double(param$w), # w(nw)
as.integer(nw), # nw
as.double(lambda), # lambda
as.double(sigma), # sigma
double(nw*mc.cores), # wad(nw,nthreds)
as.integer(mc.cores), # nthreds
th = double(n), # thn(n1*n2*n3)
sy = double(n))[c("ni","th","sy")]
}
## extract sum of weigths (see PS) and consider only voxels with ni larger then mean
th <- z$th
ni <- z$ni
ni[!mask]<-1
ind <- (ni > .9999*quantile(ni[ni>1],qni))#&(z$th>sigma*minlev)
## use correction factor for sd of NC Chi distribution
sy1 <- z$sy[ind]
th1 <- th[ind]
sy1 <- sy1/fncchis(th1/sigma,varstats)
## use the maximal mode of estimated local sd parameters, exclude largest values for better precision
dsigma <- density( sy1[sy1>0], n = 4092, adjust = hadj, to = min( max(sy1[sy1>0]), median(sy1[sy1>0])*5) )
sigma <- dsigma$x[dsigma$y == max(dsigma$y)][1]
if (sequence) {
sigmas[i] <- sigma
lind[i] <- sum(ind)
minni[i] <- min(ni[ind])
}
}
## END PS iteration
## this is the result (th is expectation, not the non-centrality parameter !!!)
result <- list(sigma = if(sequence) sigmas else sigma,
theta = th,
lind = if(sequence) lind else sum(ind),
minni = if(sequence) minni else min(ni[ind]))
}
##
## estimation of sigma from background or brain area
## various methods from table 2 in Aja-Fernandez 2009 following Aja-Fernandez
##
if(method%in%c("AFmodevn","AFmodem1chi","AFbkm2chi","AFbkm1chi")){
if(method!="AFmodevn"&is.null(level)&is.null(mask)) {
stop("need information on background using either level or mask")
}
## dimension and size of cubus
ddim <- dim(y)
n <- prod(ddim)
## test dimension
if (length(ddim) != 3) stop("first argument should be a 3-dimentional array")
## check mask
if (is.null(mask)) mask <- array(TRUE, ddim)
if(!is.null(level)){
if (method=="AFmodevn"){
mask[y<level] <- FALSE
} else {
mask[y>level] <- FALSE
}
}
if(length(mask) != n) stop("dimensions of data array and mask should coincide")
## let FORTRAN do the calculation
if(method%in%c("AFmodevn","AFmodem1chi")){
if(method=="AFmodevn"){
sigma <- .Fortran(C_afmodevn,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask),
as.double(hinit),
as.double(vext),
sigma = double(n))$sigma
sigma <- sigma/2/(ncoils-gamma(ncoils+.5)^2/gamma(ncoils)^2)
sigma <- array( sqrt(sigma), ddim)
} else {
afactor <- sqrt(1/2)*gamma(ncoils)/gamma(ncoils+.5)
sigma <- .Fortran(C_afmodem1,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask),
as.double(hinit),
as.double(vext),
sigma = double(n))$sigma
sigma <- array( afactor*sigma, ddim)
}
## use the maximal mode of estimated local variance parameters, exclude largest values for better precision
dsigma <- density( sigma[sigma>0], n = 4092, adjust = hadj, to = min( max(sigma[sigma>0]), median(sigma[sigma>0])*5) )
sigmag <- dsigma$x[dsigma$y == max(dsigma$y)][1]
} else {
if(method=="AFbkm2chi"){
sigmag <- sqrt(mean(y[mask]^2)/2/ncoils)
} else {
sigmag <- mean(y[mask])*sqrt(ncoils/2)*gamma(ncoils)/gamma(ncoils+.5)/sqrt(ncoils)
}
}
result <- list(sigma = sigmag)
}
class(result) <- "NCchiGlobalSigma"
result
}
awslinsd <- function(y,hmax=NULL,hpre=NULL,h0=NULL,mask=NULL,
ladjust=1,wghts=NULL,varprop=.1,A0,A1)
{
#
# first check arguments and initialize
#
homogen <- TRUE
wghts <- NULL
args <- match.call()
dy<-dim(y)
if(length(dy)!=3) stop("Image should be 3D")
n1 <- dy[1]
n2 <- dy[2]
n3 <- dy[3]
#
# set appropriate defaults
#
if(is.null(wghts)) wghts <- c(1,1,1)
wghts <- wghts[1]/wghts[2:3]
cpar<-setawsdefaults(mean(y),ladjust,hmax,wghts)
qlambda <- .98
if(is.null(hmax)) hmax <- 5
lambda <- ladjust*qchisq(qlambda,1)*2
maxvol <- getvofh(hmax,c(1,0,0,1,0,1),c(1,wghts))
kstar <- as.integer(log(maxvol)/log(1.25))
k <- 6
cat("Estimating variance model using PS with with lambda=",signif(lambda,3)," hmax=",hmax,"number of iterations:",kstar-k+1,"\n")
if(is.null(mask)) {
if(length(dy)==0) mask <- rep(TRUE,length(y)) else mask <- array(TRUE,dy)
}
n<-length(y)
dmask <- dim(mask)
nvoxel <- sum(mask)
position <- array(0,dmask)
position[mask] <- 1:nvoxel
y <- y[mask]
#
# family dependent transformations
#
sigma2 <- max(1,IQRdiff(as.vector(y))^2)
if(any(h0)>0) sigma2<-sigma2*Varcor.gauss(h0)
# cat("Estimated variance: ", signif(sigma2,4),"\n")
sigma2 <- rep(sigma2, nvoxel)
sigma2 <- 1/sigma2
# taking the invers yields simpler formulaes
# now check which procedure is appropriate
## this is the version on a grid
#
# Initialize for the iteration
#
zobj<-list(ai=y, bi= rep(1,nvoxel), theta= y)
mae<-NULL
lambda0<-1e50 # that removes the stochstic term for the first step, initialization by kernel estimates
#
# produce a presmoothed estimate to stabilze variance estimates
#
if(is.null(hpre)) hpre<-20^(1/3)
dlw<-(2*trunc(hpre/c(1,wghts))+1)[1:3]
hobj <- .Fortran(C_smooth3d,
as.double(y[mask]),
as.double(rep(1,nvoxel)),
as.integer(position),
as.integer(0L),
as.integer(nvoxel),
as.integer(n1),
as.integer(n2),
as.integer(n3),
as.integer(1L),
as.double(hpre),
theta=as.double(zobj$theta),
bi=as.double(zobj$bi),
as.integer(2L), # lkern
double(prod(dlw)),
as.double(wghts),
double(1))[c("bi","theta")]
theta <- bi <- array(0,dy)
theta[mask] <- hobj$theta
bi[mask] <- hobj$bi
hobj <- list(bi=bi, theta=theta)
#
# iteratate until maximal bandwidth is reached
#
# cat("Progress:")
# total <- cumsum(1.25^(1:kstar))/sum(1.25^(1:kstar))
pb <- txtProgressBar(0, kstar, style = 3)
while (k<=kstar) {
hakt0 <- gethani(1,10,1.25^(k-1),c(1,0,0,1,0,1),c(1,wghts),1e-4)
hakt <- gethani(1,10,1.25^k,c(1,0,0,1,0,1),c(1,wghts),1e-4)
dlw<-(2*trunc(hakt/c(1,wghts))+1)[1:3]
if(any(h0>0)) lambda0<-lambda0*Spatialvar.gauss(hakt0/0.42445/4,h0,3)/Spatialvar.gauss(hakt0/0.42445/4,1e-5,3)
# Correction for spatial correlation depends on h^{(k)}
hakt0<-hakt
# heteroskedastic Gaussian case
zobj <- .Fortran(C_cgaws,
as.double(y),
as.integer(position),
as.double(sigma2),
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
bi = as.double(zobj$bi),
bi2 = double(n),
bi0 = as.double(zobj$bi0),
gi = double(n),
gi2 = double(n),
ai = as.double(zobj$ai),
as.integer(2L),
as.double(0.25),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi0", "bi2", "gi2", "ai", "gi", "hakt")]
zobj$theta <- (zobj$ai/zobj$bi)
#
# Calculate MAE and MSE if true parameters are given in u
# this is for demonstration and testing for propagation (parameter adjustments)
# only.
#
# Prepare for next iteration
#
#
# Create new variance estimate
#
vobj <- awsgsigmasd(y,hobj,zobj,varprop,h0,A0,A1)
sigma2 <- vobj$sigma2inv
coef <- vobj$coef
rm(vobj)
lambda0<-lambda
setTxtProgressBar(pb, k)
# if (max(total) >0) {
# cat(signif(total[k],2)*100,"% . ",sep="")
# }
k <- k+1
gc()
}
close(pb)
# cat("\n")
theta <-array(0, dmask)
theta[mask] <- zobj$theta
###
### end iterations now prepare results
###
list(theta = theta, vcoef=coef, mask=mask)
}
############################################################################
#
# estimate inverse of variances, uses nonadaptive hobj to stabilize,
# based on absolute residuals, linear
# expects only voxel within mask
#
############################################################################
awsgsigmasd <- function(y,hobj,tobj,varprop,h0,thmin,thmax){
## specify sd to be linear to mean
thrange <- range(y)
# thmax <- thrange[2]-.2*diff(thrange)
# thmin <- thrange[1]+.1*diff(thrange)
ind <- tobj$gi>1.5&tobj$theta>thmin&tobj$theta<thmax
absresid <- abs((y-tobj$theta)[ind]*tobj$gi[ind]/sqrt(tobj$gi[ind]^2-tobj$gi2[ind]))/.8
# absresid <- abs((y-tobj$theta)[ind])/.8
theta <- tobj$theta[ind]
wght <- (tobj$gi[ind]^2-tobj$gi2[ind])/tobj$gi[ind]^2
coef <- coefficients(lm(absresid~theta,weights=wght^2))
# force positive variance for positive mean by increasing variance estimate
if(coef[2] < 0){
coef[1] <- coefficients(lm(absresid~1,weights=wght^2))
coef[2] <- 0
}
if(coef[1] <0.5){
coef[2] <- coefficients(lm(absresid~theta-1,weights=wght^2))
coef[1] <- 0
}
gamma <- pmin(tobj$gi/hobj$bi,1)
theta <- gamma*tobj$theta+(1-gamma)*hobj$theta
sigma2 <- (coef[1]+coef[2]*theta)^2
varquantile <- varprop*mean(sigma2)
sigma2 <- pmax(sigma2,varquantile)
# cat("Estimated mean variance",signif(mean(sigma2),3)," Variance parameters:",signif(coef,3),"\n")
list(sigma2inv=1/sigma2,coef=coef)
}
WIAS-BERLIN/aws documentation built on Sept. 10, 2023, 6:20 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions awsgsigmasd awslinsd estGlobalSigma AFLocalSigma iniSigmaNCchi IQQdiff awsLocalSigma estimateSigmaCompl
Documented in AFLocalSigma awslinsd awsLocalSigma estGlobalSigma estimateSigmaCompl
estimateSigmaCompl <- function(magnitude,phase,mask,kstar=20,kmin=8,hsig=5,lambda=12,verbose=TRUE){
## kmin = 10 corresponds to an initial bandwidth of 1.47 giving positive weight to direct neighbors and
## 2D diagonal neigbors
args <- sys.call(-1)
sdim <- dim(mask)
# if(!is.numeric(magnitude)){
# if (verbose) cat("reading Magnitude file ... ")
# R <- readNIfTI(magnitude, reorient = FALSE)
# } else {
R <- magnitude
# }
# if(!is.numeric(phase)){
# if (verbose) cat("reading Phase file ... ")
# Ph <- readNIfTI(phase, reorient = FALSE)
# } else {
Ph <- phase
# }
ComplImg <- array(0,c(2,sdim))
ComplImg[1,,,] <- R*cos(Ph)
ComplImg[2,,,] <- R*sin(Ph)
## find the number of usable cores
mc.cores <- setCores(, reprt = FALSE)
##
## start smoothing and variance estimation
##
n <- prod(sdim)
lambda0 <- 1e40
sigma2 <- array(1e10,sdim)
# just inilitialize with something large, first step is nonadaptive due to lambda0
k <- kmin
hmax <- 1.25^(kstar/3)
## preparations for median smoothing
parammd <- getparam3d(hsig,c(1,1))
nwmd <- length(parammd$w)
if (verbose) pb <- txtProgressBar(min = 0, max = kstar-kmin+1, style = 3)
bi <- array(1,sdim)
zobj <- list(theta=ComplImg, bi=bi)
if (verbose) {
mae <- NULL
protocol <- matrix("", kstar-kmin+1, 1, dimnames = list(paste("step", kmin:kstar), "protocol"))
}
while (k <= kstar) {
## determine the actual bandwidth for this step
hakt <- gethani(1, 1.25*hmax, 2, 1.25^k, c(1,1), 1e-4)
## we need the (approx.) size of the weigthing scheme array
dlw <- (2*trunc(hakt/c(1, 1, 1))+1)[1:3]
## perform the actual adaptive smoothing
zobj <- .Fortran(C_cplxawss,
as.double(ComplImg),
as.logical(mask),
as.integer(2),
as.integer(sdim[1]),
as.integer(sdim[2]),
as.integer(sdim[3]),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(sigma2),
bi = as.double(zobj$bi),
theta = double(2*n),
sigma2 = double(n),
as.integer(mc.cores),
double(prod(dlw)),
as.double(c(1,1)),
double(2 * mc.cores))[c("bi", "theta", "hakt","sigma2")]
##
## now get local median variance estimates
##
dim(zobj$sigma2) <- sdim
sigma2 <- .Fortran(C_mediansm,
as.double(zobj$sigma2),
as.logical(mask),
as.integer(sdim[1]),
as.integer(sdim[2]),
as.integer(sdim[3]),
as.integer(parammd$ind),
as.integer(nwmd),
double(nwmd*mc.cores), # work(nw,nthreds)
as.integer(mc.cores),
sigma2n = double(n))$sigma2n/0.6931
# sigma2n containes sum of 2 independent squared residuals
# 0.6931 approximates median \chi_2 /2
# needed to get correct results
## use maximum ni
bi <- zobj$bi <- pmax(bi, zobj$bi)
## some verbose stuff
if (verbose) {
protocol[k-kmin+1,1] <- paste("bandwidth: ", signif(hakt, 3),
"sigma: mean: ", signif(sqrt(mean(sigma2[mask])),3),
"median: ", signif(sqrt(median(sigma2[mask])),3),
"sd: ", signif(sd(sqrt(sigma2[mask])),3),
"median(bi):", signif(median(zobj$bi[mask]),3),
"max(bi):", signif(max(zobj$bi[mask]),3))
setTxtProgressBar(pb, k-kmin+1)
}
## go for next iteration
k <- k+1
lambda0 <- lambda
gc()
}
dim(zobj$theta) <- c(2,sdim)
# return estimated parameters of rician distribution
list(sigma=array(sqrt(sigma2),sdim),
theta=array(sqrt(zobj$theta[1,,,]^2+zobj$theta[2,,,]^2),sdim),
sigmal=array(sqrt(zobj$sigma2),sdim),mask=mask,
protocol=protocol,args=args)
}
awsLocalSigma <- function(y, steps, mask, ncoils, vext=c(1,1),
lambda = 5, minni=2, hsig=5, sigma=NULL, family=c("NCchi","Gauss"),
verbose = FALSE, trace = FALSE, u=NULL){
if(trace) tergs <- array(0,c(steps,4,sum(mask))) else tergs <- NULL
family <- match.arg(family[1],c("NCchi","Gauss","Gaussian"))
if(family == "NCchi"){
varstats <- sofmchi(ncoils)
th <- seq(0,30,.01)
z <- fncchiv(th,varstats)
minz <- min(z)
th <- th[z>min(z)]
minth <- min(th)
z <- z[z>min(z)]
nz <- nls(z~(a*th+b*th*th+c*th*th*th+d)/(a*th+b*th*th+c*th*th*th+1+d),data=list(th=th,z=z),start=list(a=1,b=1,c=1,d=1))
vpar <- c(minth,minz,coefficients(nz))
## this provides an excellent approximation for the variance reduction in case of low ncp
}
if(length(vext)==3) vext <- vext[2:3]/vext[1]
## dimension and size of cubus
ddim <- dim(y)
n <- prod(ddim)
## test dimension
if (length(ddim) != 3) stop("first argument should be a 3-dimentional array")
## check mask
if (is.null(mask)) mask <- array(TRUE, ddim)
if(length(mask) != n) stop("dimensions of data array and mask should coincide")
## initial value for sigma_0
if(is.null(sigma)){
# sigma <- sqrt( mean( y[mask]^2) / 2 / ncoils)
sigma <- IQQdiff( y, mask, .25, verbose=verbose)
if("NCchi" == family){
## iterative improvement for NC chi-distribution
sigma <- iniSigmaNCchi( y, mask, .25, ncoils, sigma)
sigma <- iniSigmaNCchi( y, mask, .25, ncoils, sigma)
sigma <- iniSigmaNCchi( y, mask, .25, ncoils, sigma)
}
}
##
## Prepare for diagnostics plots
##
if(verbose){
mslice <- (ddim[3]+1)/2
ymslice <- y[,,mslice]
ymslice[!mask[,,mslice]] <- 0
if(!is.null(u)&&"NCchi" == family){
par(mfrow=c(2,4),mar=c(3,3,3,1),mgp=c(2,1,0))
} else {
par(mfrow=c(2,3),mar=c(3,3,3,1),mgp=c(2,1,0))
}
} else {
cat("step")
}
## define initial arrays for parameter estimates and sum of weights (see PS)
th <- y
ksi <- array( y^2, ddim)
ni <- array( 1, ddim)
sigma <- array(sigma, ddim)
# initialize array for local sigma by global estimate
mc.cores <- setCores(,reprt=FALSE)
## preparations for median smoothing
parammd <- getparam3d(hsig, vext)
## iterate PS starting with bandwidth h0
if(family=="NCchi"){
# precompute values of log(besselI) for interpolation
nfb <- 200 * ncoils # for arguments > nfb asymp. approx can be used
x <- 1:nfb
flb <- x + log(besselI(x, ncoils-1, TRUE))
}
for (i in 1:steps) {
h <- 1.25^((i-1)/3)
param <- getparam3d(h, vext)
nw <- length(param$w)
## perform one step PS with bandwidth h
if(family=="NCchi"){
z <- .Fortran(C_awslchi2,
as.double(y),# data
as.double(ksi),# \sum w_ij S_j^2
ni = as.double(ni),
as.double(sigma),
as.double(vpar),# parameters for var. reduction
as.double(ncoils),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(param$ind),
as.double(param$w),
as.integer(nw),
as.double(minni),
double(nw*mc.cores), # wad(nw,nthreds)
double(nw*mc.cores), # sad(nw,nthreds)
as.double(lambda),
as.integer(mc.cores),
as.integer(floor(ncoils)),
double(floor(ncoils)*mc.cores), # work(L,nthreds)
th = double(n),
sigman = double(n),
ksi = double(n),
as.double(flb),
as.integer(nfb))[c("ni","ksi","th","sigman")]
thchi <- z$th
ksi <- z$ksi
thchi[!mask] <- 0
} else{
## assume Gaussian data
z <- .Fortran(C_awslgaus,
as.double(y),# data
as.double(th),# \sum w_ij S_j
ni = as.double(ni),
as.double(sigma),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(param$ind),
as.double(param$w),
as.integer(nw),
as.double(minni),
as.double(lambda),
th = double(n),
sigman = double(n))[c("ni","th","sigman")]
}
## extract sum of weigths (see PS) and consider only voxels with ni larger then mean
th <- array(z$th,ddim)
ni <- array(z$ni,ddim)
z$sigman[z$sigman==0] <- median(z$sigman[z$sigman>0])
nmask <- sum(mask)
if(verbose) cat("local estimation in step ",i," completed",format(Sys.time()),"\n")
##
## nonadaptive smoothing of estimated standard deviations
##
if(any(ni[mask]>minni)){
z$sigman[mask] <- z$sigman[mask]*(1+runif(nmask,-.0001,.0001))
##
## avoid ties in local neighborhood cubes of z$sigman
##
nwmd <- length(parammd$w)
sigma <- .Fortran(C_mediansm,
as.double(z$sigman),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(parammd$ind),
as.integer(nwmd),
double(nwmd*mc.cores), # work(nw,nthreds)
as.integer(mc.cores),
sigman = double(n))$sigman
}
dim(sigma) <- ddim
mask[sigma==0] <- FALSE
if(verbose) cat("local median smoother in step ",i," completed",format(Sys.time()),"\n")
##
## diagnostics
##
if(verbose){
meds <- median(sigma[mask])
means <- mean(sigma[mask])
image(ymslice,col=grey(0:255/255))
title(paste("S max=",signif(max(y[mask]),3)," median=",signif(median(y[mask]),3)))
image(th[,,mslice],col=grey(0:255/255))
title(paste("E(S) max=",signif(max(th[mask]),3)," median=",signif(median(th[mask]),3)))
image(sigma[,,mslice],col=grey(0:255/255),zlim=c(0,max(sigma[mask])))
title(paste("sigma max=",signif(max(sigma[mask]),3)," median=",signif(meds,3)))
image(ni[,,mslice],col=grey(0:255/255))
title(paste("Ni max=",signif(max(ni[mask]),3)," median=",signif(median(ni[mask]),3)))
plot(density(sigma[mask]),main="density of sigma")
plot(density(ni[mask]),main="density of Ni")
cat("mean sigma",means,"median sigma",meds,"sd sigma",sd(sigma[mask]),"\n")
if(!is.null(u)&&"NCchi" == family){
thchims <- fncchir(th/sigma,varstats)*sigma
thchims[!mask] <- u[!mask]
image(abs(thchims-u)[,,mslice],col=grey(0:255/255))
title("abs Error in thchims")
plot(density((thchims-u)[mask]),main="density of thchims-u")
cat("MAE(th)",mean(abs(thchi-u)[mask]),"RMSE(th)",sqrt(mean((thchi-u)[mask]^2)),"MAE(thms)",mean(abs(thchims-u)[mask]),"RMSE(thms)",sqrt(mean((thchims-u)[mask]^2)),"\n")
}
} else {
cat(" ",i)
}
if(trace) {
tergs[i,1,] <- ni[mask]
tergs[i,2,] <- th[mask]
tergs[i,3,] <- z$sigma[mask]
tergs[i,4,] <- sigma[mask]
}
}
## END PS iteration
if(!verbose) cat("\n")
if(family=="NCchi"){
thchi <- fncchir(th/sigma,varstats)*sigma
thchi[!mask] <- 0
} else {
thchi <- NULL
}
## this is the result (th is expectation, not the non-centrality parameter !!!)
invisible(list(sigma = sigma,
sigmal = array(z$sigman,ddim),
theta = th,
thchi = thchi,
ni = ni,
tergs = tergs,
mask = mask))
}
IQQ <- function (x, q = .25, na.rm = FALSE, type = 7){
cqz <- qnorm(.05)/qnorm(q)
x <- as.numeric(x)
z <- diff(quantile(x, c(q, 1-q), na.rm = na.rm, names = FALSE, type = type))
z0 <- 0
while(abs(z-z0)>1e-5*z0){
# outlier removal
z0 <- z
# cat(sum(x>(z0*cqz))," ")
x <- x[x<(z0*cqz)]
z <- diff(quantile(x, c(q, 1-q), na.rm = na.rm, names = FALSE, type = type))
}
z
}
IQQdiff <- function(y, mask, q = .25, verbose = FALSE) {
cq <- qnorm(1-q)*sqrt(2)*2
sx <- IQQ( diff(y[mask]), q)/cq
sy <- IQQ( diff(aperm(y,c(2,1,3))[aperm(mask,c(2,1,3))]), q)/cq
sz <- IQQ( diff(aperm(y,c(3,1,2))[aperm(mask,c(3,1,2))]), q)/cq
if(verbose) cat( "Pilot estimates of sigma", sx, sy, sz, "\n")
min( sx, sy, sz)
}
iniSigmaNCchi <- function(y, mask, q, L, sigma, verbose=FALSE){
#
# robust improvement for sigma using moment estimates and IQQdiff
#
meany <- mean(y[mask]^2)
eta <- sqrt(max( 0, meany/sigma^2-2*L))
m <- sqrt(pi/2)*gamma(L+1/2)/gamma(L)/gamma(3/2)*hg1f1(-1/2,L,-eta^2/2)
v <- max( .01, 2*L+eta^2-m^2)
if(verbose) cat( eta, m, v, "\n")
IQQdiff( y, mask, q)/sqrt(v)
}
AFLocalSigma <- function(y,ncoils,level=NULL,mask=NULL,h=2,hadj=1,vext = c( 1, 1)){
##
## estimate effective sigma and effective ncoils (L) according to Aja-Fernandez MRI 2013
##
# sh2B stands for \hat{\sigma}^2_{n,B}
# sh2L stands for \hat{\sigma}^2_{nL}
# sh2eB stands for \hat{\sigma}^2_{eff,B}
# sh2eS stands for \hat{\sigma}^2_{eff,S}
# LeB stands for \hat{L}_{eff,B}
# vmlb - Var(M_L(x)|x \in B)
# vmlbx - Var(M_L(x)|x \in B)_x
# phix - \hat{\Phi}_n(x)
# seff - 3D array of local variance estimates
# Leff - 3D array of effective L
ddim <- dim(y)
n <- prod(ddim)
if(is.null(level)&is.null(mask)){
warning("no background definition, need either level or mask")
return(invisible(NULL))
}
if(is.null(level)) indB <- !mask else indB <- y<level
sh2Lsimple <- mean(y[indB]^2)/2
mask1 <- array(TRUE,ddim)
##
## local variance estimates in vx
##
vx <- .Fortran(C_afmodevn,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask1),
as.double(h),
as.double(vext),
sigma = double(n))$sigma
dim(vx) <- ddim
vxb <- vx[indB]
vxs <- vx[!indB]
dv2b <- density( vxb[vxb>0], n = 4092, adjust = hadj, to = min( max(vxb[vxb>0]), median(vxb[vxb>0])*5) )
sh2B <- dv2b$x[dv2b$y == max(dv2b$y)][1]
dv2s <- density( vxs[vxs>0], n = 4092, adjust = hadj, to = min( max(vxs[vxs>0]), median(vxs[vxs>0])*5) )
sh2eS <- dv2s$x[dv2s$y == max(dv2s$y)][1]
m2 <- .Fortran(C_afmodem2,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask1),
as.double(h),
as.double(vext),
sm = double(n))$sm
dim(m2) <- ddim
m2b <- m2[indB]
dm2b <- density( m2b[m2b>0], n = 4092, adjust = hadj, to = min( max(m2b[m2b>0]), median(m2b[m2b>0])*5) )
sh2L <- dm2b$x[dm2b$y == max(dm2b$y)][1]/2
cat("sh2B",sh2B,"sh2Lsimple",sh2Lsimple,"sh2L",sh2L,"\n")
##
## now get sh2eB and LeB
##
LeBn <- ncoils
sh2eBn <- sh2B/2/(LeBn-gamma(LeBn+.5)^2/gamma(LeBn)^2)
LeB <- 1
sh2eB <- 0
while(abs(sh2eB-sh2eBn)>1e-5||abs(LeB-LeBn)>1e-5){
LeB <- LeBn
sh2eB <- sh2eBn
LeBn <- sh2L/sh2eB
sh2eBn <- sh2B/2/(LeBn-gamma(LeBn+.5)^2/gamma(LeBn)^2)
cat("LeB",LeB,"LeBn",LeBn,"sh2eB",sh2eB,"sh2eBn",sh2eBn,"\n")
}
sh2eB <- sh2eBn
LeB <- sh2L/sh2eB
phix <- pmax(0,pmin(1,sh2L/(m2-sh2L)))
dim(phix) <- dim(m2)
seff <- sqrt((1-phix)*sh2eS+phix*sh2eB)
Leff <- sh2L/seff^2
list(sigma=seff,Leff=Leff,sh2L=sh2L,sh2B=sh2B,sh2eS=sh2eS,sh2eB=sh2eB,LeB=LeB,m2=m2,vx=vx,indB=indB)
}
estGlobalSigma <- function(y, mask=NULL, ncoils=1, steps=16, vext=c(1,1),
lambda=20, hinit=2, hadj=1, q=.25, level=NULL,
sequence=FALSE,method=c("awsVar","awsMAD",
"AFmodevn","AFmodem1chi","AFbkm2chi","AFbkm1chi")){
##
## estimate global scale parameter sigma for NCchi distributed data
##
method <- match.arg(method)
qni <- .8
## methods using the propagation-separation (PS) approach
if(method %in% c("awsVar","awsMAD")){
varstats <- sofmchi(ncoils)
if(length(vext)==3) vext <- vext[2:3]/vext[1]
## dimension and size of cubus
ddim <- dim(y)
n <- prod(ddim)
## test dimension
if (length(ddim) != 3) stop("first argument should be a 3-dimentional array")
## check mask
if (is.null(mask)) mask <- array(TRUE, ddim)
if(length(mask) != n) stop("dimensions of data array and mask should coincide")
## initial value for sigma_0
# sigma <- sqrt( mean( y[mask]^2) / 2 / ncoils)
sigma <- IQQdiff( y, mask, q)
# cat( "sigmahat1", sigma, "\n")
sigma <- iniSigmaNCchi( y, mask, q, ncoils, sigma)
# cat( "sigmahat2", sigma, "\n")
sigma <- iniSigmaNCchi( y, mask, q, ncoils, sigma)
# cat( "sigmahat3", sigma, "\n")
sigma <- iniSigmaNCchi( y, mask, q, ncoils, sigma)
# cat( "sigmahat4", sigma,"\n")
## define initial arrays for parameter estimates and sum of weights (see PS)
th <- array( 1, ddim)
ni <- array( 1, ddim)
# y <- y/sigma # rescale to avoid passing sigma to awsvchi
minlev <- sqrt(2)*gamma(ncoils+.5)/gamma(ncoils)
if (sequence) sigmas <- lind <- minni <- numeric(steps)
mc.cores <- setCores(,reprt=FALSE)
## iterate PS starting with bandwidth hinit
for (i in 1:steps) {
h <- hinit * 1.25^((i-1)/3)
param <- getparam3d(h,vext)
nw <- length(param$w)
fncchi <- fncchiv(th/sigma,varstats)
## correction factor for variance of NC Chi distribution
## perform one step PS with bandwidth h
if(method=="awsVAR"){
z <- .Fortran(C_awsvchi,
as.double(y), # data
as.double(th), # previous estimates
ni = as.double(ni),
as.double(fncchi/2),
as.integer(mask),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(param$ind),
as.double(param$w),
as.integer(nw),
as.double(lambda),
as.double(sigma),
th = double(n),
sy = double(n))[c("ni","th","sy")]
} else { # method=="awsMAD"
z <- .Fortran(C_awsadchi,
as.double(y), # y(n1,n2,n3)
as.double(th), # th(n1,n2,n3)
ni = as.double(ni), # ni(n1,n2,n3)
as.double(fncchi/2), # fns(n1,n2,n3)
as.integer(mask), # mask(n1,n2,n3)
as.integer(ddim[1]), # n1
as.integer(ddim[2]), # n2
as.integer(ddim[3]), # n3
as.integer(param$ind), # ind(3,nw)
as.double(param$w), # w(nw)
as.integer(nw), # nw
as.double(lambda), # lambda
as.double(sigma), # sigma
double(nw*mc.cores), # wad(nw,nthreds)
as.integer(mc.cores), # nthreds
th = double(n), # thn(n1*n2*n3)
sy = double(n))[c("ni","th","sy")]
}
## extract sum of weigths (see PS) and consider only voxels with ni larger then mean
th <- z$th
ni <- z$ni
ni[!mask]<-1
ind <- (ni > .9999*quantile(ni[ni>1],qni))#&(z$th>sigma*minlev)
## use correction factor for sd of NC Chi distribution
sy1 <- z$sy[ind]
th1 <- th[ind]
sy1 <- sy1/fncchis(th1/sigma,varstats)
## use the maximal mode of estimated local sd parameters, exclude largest values for better precision
dsigma <- density( sy1[sy1>0], n = 4092, adjust = hadj, to = min( max(sy1[sy1>0]), median(sy1[sy1>0])*5) )
sigma <- dsigma$x[dsigma$y == max(dsigma$y)][1]
if (sequence) {
sigmas[i] <- sigma
lind[i] <- sum(ind)
minni[i] <- min(ni[ind])
}
}
## END PS iteration
## this is the result (th is expectation, not the non-centrality parameter !!!)
result <- list(sigma = if(sequence) sigmas else sigma,
theta = th,
lind = if(sequence) lind else sum(ind),
minni = if(sequence) minni else min(ni[ind]))
}
##
## estimation of sigma from background or brain area
## various methods from table 2 in Aja-Fernandez 2009 following Aja-Fernandez
##
if(method%in%c("AFmodevn","AFmodem1chi","AFbkm2chi","AFbkm1chi")){
if(method!="AFmodevn"&is.null(level)&is.null(mask)) {
stop("need information on background using either level or mask")
}
## dimension and size of cubus
ddim <- dim(y)
n <- prod(ddim)
## test dimension
if (length(ddim) != 3) stop("first argument should be a 3-dimentional array")
## check mask
if (is.null(mask)) mask <- array(TRUE, ddim)
if(!is.null(level)){
if (method=="AFmodevn"){
mask[y<level] <- FALSE
} else {
mask[y>level] <- FALSE
}
}
if(length(mask) != n) stop("dimensions of data array and mask should coincide")
## let FORTRAN do the calculation
if(method%in%c("AFmodevn","AFmodem1chi")){
if(method=="AFmodevn"){
sigma <- .Fortran(C_afmodevn,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask),
as.double(hinit),
as.double(vext),
sigma = double(n))$sigma
sigma <- sigma/2/(ncoils-gamma(ncoils+.5)^2/gamma(ncoils)^2)
sigma <- array( sqrt(sigma), ddim)
} else {
afactor <- sqrt(1/2)*gamma(ncoils)/gamma(ncoils+.5)
sigma <- .Fortran(C_afmodem1,
as.double(y),
as.integer(ddim[1]),
as.integer(ddim[2]),
as.integer(ddim[3]),
as.integer(mask),
as.double(hinit),
as.double(vext),
sigma = double(n))$sigma
sigma <- array( afactor*sigma, ddim)
}
## use the maximal mode of estimated local variance parameters, exclude largest values for better precision
dsigma <- density( sigma[sigma>0], n = 4092, adjust = hadj, to = min( max(sigma[sigma>0]), median(sigma[sigma>0])*5) )
sigmag <- dsigma$x[dsigma$y == max(dsigma$y)][1]
} else {
if(method=="AFbkm2chi"){
sigmag <- sqrt(mean(y[mask]^2)/2/ncoils)
} else {
sigmag <- mean(y[mask])*sqrt(ncoils/2)*gamma(ncoils)/gamma(ncoils+.5)/sqrt(ncoils)
}
}
result <- list(sigma = sigmag)
}
class(result) <- "NCchiGlobalSigma"
result
}
awslinsd <- function(y,hmax=NULL,hpre=NULL,h0=NULL,mask=NULL,
ladjust=1,wghts=NULL,varprop=.1,A0,A1)
{
#
# first check arguments and initialize
#
homogen <- TRUE
wghts <- NULL
args <- match.call()
dy<-dim(y)
if(length(dy)!=3) stop("Image should be 3D")
n1 <- dy[1]
n2 <- dy[2]
n3 <- dy[3]
#
# set appropriate defaults
#
if(is.null(wghts)) wghts <- c(1,1,1)
wghts <- wghts[1]/wghts[2:3]
cpar<-setawsdefaults(mean(y),ladjust,hmax,wghts)
qlambda <- .98
if(is.null(hmax)) hmax <- 5
lambda <- ladjust*qchisq(qlambda,1)*2
maxvol <- getvofh(hmax,c(1,0,0,1,0,1),c(1,wghts))
kstar <- as.integer(log(maxvol)/log(1.25))
k <- 6
cat("Estimating variance model using PS with with lambda=",signif(lambda,3)," hmax=",hmax,"number of iterations:",kstar-k+1,"\n")
if(is.null(mask)) {
if(length(dy)==0) mask <- rep(TRUE,length(y)) else mask <- array(TRUE,dy)
}
n<-length(y)
dmask <- dim(mask)
nvoxel <- sum(mask)
position <- array(0,dmask)
position[mask] <- 1:nvoxel
y <- y[mask]
#
# family dependent transformations
#
sigma2 <- max(1,IQRdiff(as.vector(y))^2)
if(any(h0)>0) sigma2<-sigma2*Varcor.gauss(h0)
# cat("Estimated variance: ", signif(sigma2,4),"\n")
sigma2 <- rep(sigma2, nvoxel)
sigma2 <- 1/sigma2
# taking the invers yields simpler formulaes
# now check which procedure is appropriate
## this is the version on a grid
#
# Initialize for the iteration
#
zobj<-list(ai=y, bi= rep(1,nvoxel), theta= y)
mae<-NULL
lambda0<-1e50 # that removes the stochstic term for the first step, initialization by kernel estimates
#
# produce a presmoothed estimate to stabilze variance estimates
#
if(is.null(hpre)) hpre<-20^(1/3)
dlw<-(2*trunc(hpre/c(1,wghts))+1)[1:3]
hobj <- .Fortran(C_smooth3d,
as.double(y[mask]),
as.double(rep(1,nvoxel)),
as.integer(position),
as.integer(0L),
as.integer(nvoxel),
as.integer(n1),
as.integer(n2),
as.integer(n3),
as.integer(1L),
as.double(hpre),
theta=as.double(zobj$theta),
bi=as.double(zobj$bi),
as.integer(2L), # lkern
double(prod(dlw)),
as.double(wghts),
double(1))[c("bi","theta")]
theta <- bi <- array(0,dy)
theta[mask] <- hobj$theta
bi[mask] <- hobj$bi
hobj <- list(bi=bi, theta=theta)
#
# iteratate until maximal bandwidth is reached
#
# cat("Progress:")
# total <- cumsum(1.25^(1:kstar))/sum(1.25^(1:kstar))
pb <- txtProgressBar(0, kstar, style = 3)
while (k<=kstar) {
hakt0 <- gethani(1,10,1.25^(k-1),c(1,0,0,1,0,1),c(1,wghts),1e-4)
hakt <- gethani(1,10,1.25^k,c(1,0,0,1,0,1),c(1,wghts),1e-4)
dlw<-(2*trunc(hakt/c(1,wghts))+1)[1:3]
if(any(h0>0)) lambda0<-lambda0*Spatialvar.gauss(hakt0/0.42445/4,h0,3)/Spatialvar.gauss(hakt0/0.42445/4,1e-5,3)
# Correction for spatial correlation depends on h^{(k)}
hakt0<-hakt
# heteroskedastic Gaussian case
zobj <- .Fortran(C_cgaws,
as.double(y),
as.integer(position),
as.double(sigma2),
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
bi = as.double(zobj$bi),
bi2 = double(n),
bi0 = as.double(zobj$bi0),
gi = double(n),
gi2 = double(n),
ai = as.double(zobj$ai),
as.integer(2L),
as.double(0.25),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi0", "bi2", "gi2", "ai", "gi", "hakt")]
zobj$theta <- (zobj$ai/zobj$bi)
#
# Calculate MAE and MSE if true parameters are given in u
# this is for demonstration and testing for propagation (parameter adjustments)
# only.
#
# Prepare for next iteration
#
#
# Create new variance estimate
#
vobj <- awsgsigmasd(y,hobj,zobj,varprop,h0,A0,A1)
sigma2 <- vobj$sigma2inv
coef <- vobj$coef
rm(vobj)
lambda0<-lambda
setTxtProgressBar(pb, k)
# if (max(total) >0) {
# cat(signif(total[k],2)*100,"% . ",sep="")
# }
k <- k+1
gc()
}
close(pb)
# cat("\n")
theta <-array(0, dmask)
theta[mask] <- zobj$theta
###
### end iterations now prepare results
###
list(theta = theta, vcoef=coef, mask=mask)
}
############################################################################
#
# estimate inverse of variances, uses nonadaptive hobj to stabilize,
# based on absolute residuals, linear
# expects only voxel within mask
#
############################################################################
awsgsigmasd <- function(y,hobj,tobj,varprop,h0,thmin,thmax){
## specify sd to be linear to mean
thrange <- range(y)
# thmax <- thrange[2]-.2*diff(thrange)
# thmin <- thrange[1]+.1*diff(thrange)
ind <- tobj$gi>1.5&tobj$theta>thmin&tobj$theta<thmax
absresid <- abs((y-tobj$theta)[ind]*tobj$gi[ind]/sqrt(tobj$gi[ind]^2-tobj$gi2[ind]))/.8
# absresid <- abs((y-tobj$theta)[ind])/.8
theta <- tobj$theta[ind]
wght <- (tobj$gi[ind]^2-tobj$gi2[ind])/tobj$gi[ind]^2
coef <- coefficients(lm(absresid~theta,weights=wght^2))
# force positive variance for positive mean by increasing variance estimate
if(coef[2] < 0){
coef[1] <- coefficients(lm(absresid~1,weights=wght^2))
coef[2] <- 0
}
if(coef[1] <0.5){
coef[2] <- coefficients(lm(absresid~theta-1,weights=wght^2))
coef[1] <- 0
}
gamma <- pmin(tobj$gi/hobj$bi,1)
theta <- gamma*tobj$theta+(1-gamma)*hobj$theta
sigma2 <- (coef[1]+coef[2]*theta)^2
varquantile <- varprop*mean(sigma2)
sigma2 <- pmax(sigma2,varquantile)
# cat("Estimated mean variance",signif(mean(sigma2),3)," Variance parameters:",signif(coef,3),"\n")
list(sigma2inv=1/sigma2,coef=coef)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.