inst/rcode/gen_mixtens.r In dti: Analysis of Diffusion Weighted Imaging (DWI) Data

```tdatamixtens <- function(w,th,alpha,beta,grad,bvalue,sigma,ns0){
#
#  w  - matrix of wghts (m+1,n1,n2,n3)
#  th - ev-parameters (2,n1,n2,n3) lambda1 = th[1]+th[2]
#                                  lambda2 = th[2]
#  alpha, beta (m,n1,n2,n3)  spherical coordinates of directions
#
args <- list(sys.call())
ddim <- dim(w)[2:4]
m <- dim(w)[1]-1
#  check array dimensions
if(any(dim(th)!=c(2,ddim))) return("wrong dimension of th")
if(any(dim(alpha)!=c(m,ddim))) return("wrong dimension of alpha")
if(any(dim(beta)!=c(m,ddim))) return("wrong dimension of beta")
d <- array(0,c(3,m,ddim))
sal <- sin(alpha)
cal <- cos(alpha)
sbe <- sin(beta)
cbe <- cos(beta)
d[1,,,,] <- sal*cbe
d[2,,,,] <- sal*sbe
d[3,,,,] <- cal
s0 <- array(10000,ddim)
for(i in 1:ns0) sb[,,,i] <- s0
sb[,,,i] <- w[1,,,]*exp(-bvalue[i]*(th[1,,,]+th[2,,,]))
for(j in 1:m){
sb[,,,i] <- sb[,,,i]+w[j+1,,,]*exp(-bvalue[i]*(th[2,,,]+th[1,,,]*dg^2))
}
sb[,,,i] <- sb[,,,i]*s0
}
sb <- sqrt(rnorm(sb,sb,sigma*s0)^2+rnorm(sb,0,sigma*s0)^2)
ddim0 <- ddim
invisible(new("dtiData",
call = args,
si     = sb,
btb    = sweep( create.designmatrix.dti(grad), 2, bvalue, "*"),
bvalue = bvalue,
s0ind  = as.integer(1:ns0), # indices of S_0 images
ddim   = ddim,
ddim0  = ddim0,
xind   = 1:ddim[1],
yind   = 1:ddim[2],
zind   = 1:ddim[3],
level  = 1,
sdcoef = c(5.e-01, 0, 1, 2.e+04),
voxelext = c(1,1,1),
orientation = as.integer(c(0,2,5)),
rotation = diag(3),
source = "artificial")
)
}
#
#  w  - matrix of wghts (m+1,n1,n2,n3)
#  th - ev-parameters (2,n1,n2,n3) lambda1 = th[1]+th[2]
#                                  lambda2 = th[2]
#  alpha, beta (m,n1,n2,n3)  spherical coordinates of directions
#
args <- list(sys.call())
ddim <- dim(w)[2:4]
m <- dim(w)[1]-1
#  check array dimensions
if(any(dim(th)!=c(2,ddim))) return("wrong dimension of th")
if(any(dim(alpha)!=c(m,ddim))) return("wrong dimension of alpha")
if(any(dim(beta)!=c(m,ddim))) return("wrong dimension of beta")
ddim0 <- ddim
orient <- array(0,c(2,m,ddim))
orient[1,,,,] <- alpha
orient[2,,,,] <- beta

invisible(new("dwiMixtensor",
model  = "homogeneous_prolate",
call   = args,
ev     = th,
mix    = w[-1,,,,drop=FALSE],
orient = orient,
order  = array(m,ddim),
p      = 1,
th0    = array(10000,ddim),
sigma  = array(10000*sigma,ddim),
scorr  = array(0,c(1,1,1)),
bw     = c(0,0,0),
hmax   = 1,
btb    = sweep( create.designmatrix.dti(grad), 2, bvalue, "*"),
bvalue = bvalue,
s0ind  = as.integer(1:ns0),
ddim   = ddim,
ddim0  = ddim0,
xind   = 1:ddim[1],
yind   = 1:ddim[2],
zind   = 1:ddim[3],
voxelext = c(1,1,1),
level = 1,
orientation = as.integer(c(0,2,5)),
rotation = diag(3),
source = "artificial",
outlier = numeric(0),
scale = 1,
method = "mixtensor")
)
}

btb
}
if(any(obj1@ddim!=obj2@ddim)) return(warning("incompatible dimensions"))
if(!class(obj1)%in%c("dwiMixtensor","dtiTensor","dwiQball")) return(warning("incompatible class of obj1"))
if(!class(obj2)%in%c("dwiMixtensor","dtiTensor","dwiQball")) return(warning("incompatible class of obj2"))
if(!exists("icosa0")) data("polyeders")
n <- prod(obj1@ddim)
if(n>419103) poly <- min(poly,3)
if(n>1672495) poly <- min(poly,2)
if(n>6628036) poly <- min(poly,1)
polyeder <- switch(poly+1,icosa0,icosa1,icosa2,icosa3,icosa4)
if(class(obj1)=="dwiQball"){
sphdesign <- design.spheven(obj1@order,polyeder\$vertices,obj1@lambda)\$design
sphcoef <- obj1@sphcoef
dim(sphcoef) <- c(dim(sphcoef)[1],prod(dim(sphcoef)[-1]))
}
as.double(polyeder\$vertices),
as.integer(polyeder\$nv),
as.double(obj1@ev),
as.double(obj1@orient),
as.double(obj1@mix),
as.integer(obj1@order),
as.integer(dim(obj1@mix)[1]),
as.integer(n),
DUP=FALSE,
} else if(class(obj1)=="dtiTensor"){
as.double(polyeder\$vertices),
as.integer(polyeder\$nv),
as.double(obj1@D),
as.integer(n),
DUP=FALSE,
} else t(sphdesign)%*%sphcoef
if(class(obj2)=="dwiQball") {
sphdesign <- design.spheven(obj2@order,polyeder\$vertices,obj2@lambda)\$design
sphcoef <- obj2@sphcoef
dim(sphcoef) <- c(dim(sphcoef)[1],prod(dim(sphcoef)[-1]))
}
as.double(polyeder\$vertices),
as.integer(polyeder\$nv),
as.double(obj2@ev),
as.double(obj2@orient),
as.double(obj2@mix),
as.integer(obj2@order),
as.integer(dim(obj2@mix)[1]),
as.integer(n),
DUP=FALSE,
} else if(class(obj2)=="dtiTensor"){
as.double(polyeder\$vertices),
as.integer(polyeder\$nv),
as.double(obj2@D),
as.integer(n),