Nothing
#' Create Isosurface from 3D-array
#'
#' Create Isosurface from 3D-array using Marching Cubes algorithm
#'
#' @param vol an integer valued 3D-array
#' @param threshold threshold for creating the surface
#' @param from numeric: the lower threshold of a range (overrides \code{threshold})
#' @param to numeric: the upper threshold of a range (overrides \code{threshold})
#' @param spacing numeric 3D-vector: specifies the voxel dimensons in x,y,z direction.
#' @param origin numeric 3D-vector: origin of the original data set, will transpose the mesh onto that origin.
#' @param direction a 3x3 direction matrix
#' @param IJK2RAS 4x4 IJK2RAS transformation matrix
#' @param as.int logical: if TRUE, the array will be stored as integer (might decrease RAM usage)
#'
#' @return returns a triangular mesh of class "mesh3d"
#' @examples
#' #this is the example from the package "misc3d"
#' x <- seq(-2,2,len=50)
#' g <- expand.grid(x = x, y = x, z = x)
#' v <- array(g$x^4 + g$y^4 + g$z^4, rep(length(x),3))
#' storage.mode(v) <- "integer"
#' \dontrun{
#' mesh <- vcgIsosurface(v,threshold=10)
#' require(rgl)
#' wire3d(mesh)
#' ##now smooth it a little bit
#' wire3d(vcgSmooth(mesh,"HC",iteration=3),col=3)
#' }
#' @export
vcgIsosurface <- function(vol,threshold,from=NULL,to=NULL,spacing=NULL, origin=NULL,direction=NULL,IJK2RAS=diag(c(-1,-1,1,1)),as.int=FALSE) {
if (length(dim(vol)) != 3)
stop("3D array needed")
if (!is.null(from) && !is.null(to)) {
if (from > to)
stop("from must be > to")
tmpvol <- vol*0L+255L
storage.mode(tmpvol) <- "integer"
tmpvol[which(vol < from)] <- 0L
tmpvol[which(vol > to)] <- 0L
vol <- tmpvol
as.int <- TRUE
threshold <- 255
}
mirr <- FALSE
mvol <- max(vol)
minvol <- min(vol)
if (threshold == mvol)
threshold <- threshold-1e-5
else if (threshold > mvol || threshold < minvol)
stop("threshold is outside volume values")
if (as.int)
storage.mode(vol) <- "integer"
volmesh <- .Call("RMarchC",vol,threshold)
rm(vol)
gc()
volmesh$vb <- rbind(volmesh$vb,1)
volmesh$it <- volmesh$it
if (!is.null(origin))
origin <- as.vector(applyTransform(t(origin),IJK2RAS))
class(volmesh) <- "mesh3d"
if (!is.null(spacing))
volmesh$vb[1:3,] <- volmesh$vb[1:3,]*spacing
if (!is.null(direction)) {
IJK2RAS <- cbind(rbind(direction,0),c(0,0,0,1))%*%IJK2RAS
if (det(direction) < 0)
mirr <- TRUE
}
volmesh <- applyTransform(volmesh,IJK2RAS)
if (!is.null(origin))
volmesh$vb[1:3,] <- volmesh$vb[1:3,]+origin
if (mirr)
volmesh <- invertFaces(volmesh)
if (!checkFaceOrientation(volmesh))
volmesh <- invertFaces(volmesh)
return(volmesh)
}
###helpers imported from Morpho
invertFaces <- function (mesh) {
mesh$it <- mesh$it[c(3, 2, 1), ]
mesh <- vcgUpdateNormals(mesh)
return(mesh)
}
applyTransform <- function(x,trafo,inverse)UseMethod("applyTransform")
applyTransform.matrix <- function(x,trafo,inverse=FALSE) {
if (is.matrix(trafo)) {
if (ncol(trafo) == 3 && ncol(x) ==3)
trafo <- mat3x3tomat4x4(trafo)
if (inverse)
trafo <- solve(trafo)
out <-homg2mat(trafo%*%mat2homg(x))
} else {
stop("trafo must be a matrix")
}
return(out)
}
applyTransform.mesh3d <- function(x,trafo,inverse=FALSE) {
x$vb[1:3,] <- t(applyTransform(t(x$vb[1:3,]),trafo,inverse = inverse))
## case affine transformation
reflect <- FALSE
if (is.matrix(trafo)) {
if (det(trafo) < 0)
reflect <- TRUE
} else {
stop("trafo must be a matrix")
}
if (reflect) {
x <- invertFaces(x)
}
if (!is.null(x$normals))
x <- vcgUpdateNormals(x)
return(x)
}
mat3x3tomat4x4 <- function(x) {
n <- ncol(x)
x <- rbind(cbind(x,0),0);x[n+1,n+1] <-1
return(x)
}
mat2homg <- function(x) {
x <- rbind(t(x),1)
return(x)
}
homg2mat <- function(x) {
m <- nrow(x)
x <- t(x[1:(m-1),])
return(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.