R/clip.R

In MeshesTools: Some Tools for 3D Meshes

#' @title Clip a mesh
#' @description Clip a mesh to the volume bounded by another mesh.
#'
#' @param mesh a mesh given either as a list containing (at least) the fields
#'   \code{vertices} and \code{faces}, otherwise a \strong{rgl} mesh
#'   (i.e. a \code{mesh3d} object)
#' @param clipper a mesh given either as a list containing (at least) the fields
#'   \code{vertices} and \code{faces}, otherwise a \strong{rgl} mesh
#'   (i.e. a \code{mesh3d} object)
#' @param clipVolume Boolean, whether the clipping has to be done on the volume
#'   bounded by \code{mesh} rather than on its surface (i.e. \code{mesh} will be
#'   kept closed if it is closed)
#' @param normals Boolean, whether to compute the vertex normals of the
#'   output mesh
#'
#' @return A triangle mesh of class \code{cgalMesh} (see 
#'   \code{\link[PolygonSoup:Mesh]{Mesh}} for details).
#'
#' @export
#'
#' @note The clipping mesh (\code{clipper}) must be closed.
#'
#' @examples
#' # cube clipped to sphere
#' library(MeshesTools)
#' library(rgl)
#' mesh    <- cube3d()
#' clipper <- sphereMesh(r= sqrt(2))
#' clippedMesh <- clipMesh(mesh, clipper)
#' open3d(windowRect = c(50, 50, 562, 562))
#' view3d(zoom = 0.9)
#' shade3d(toRGL(clippedMesh), color = "purple")
#'
#' # Togliatti surface clipped to a ball ####
#' library(rmarchingcubes)
#' library(rgl)
#' library(MeshesTools)
#' 
#' # Togliatti surface equation: f(x,y,z) = 0
#' f <- function(x, y, z) {
#'   64*(x-1) *
#'     (x^4 - 4*x^3 - 10*x^2*y^2 - 4*x^2 + 16*x - 20*x*y^2 + 5*y^4 + 16 - 20*y^2) - 
#'     5*sqrt(5-sqrt(5))*(2*z - sqrt(5-sqrt(5))) * 
#'     (4*(x^2 + y^2 - z^2) + (1 + 3*sqrt(5)))^2
#' }
#' 
#' # grid
#' n <- 200L
#' x <- y <- seq(-5, 5, length.out = n)
#' z <- seq(-4, 4, length.out = n)
#' Grid <- expand.grid(X = x, Y = y, Z = z)
#' # calculate voxel
#' voxel <- array(with(Grid, f(X, Y, Z)), dim = c(n, n, n))
#' # calculate isosurface
#' contour_shape <- contour3d(
#'   griddata = voxel, level = 0, x = x, y = y, z = z
#' )
#' # make rgl mesh (plotted later)
#' mesh <- tmesh3d(
#'   vertices = t(contour_shape[["vertices"]]),
#'   indices  = t(contour_shape[["triangles"]]),
#'   normals  = contour_shape[["normals"]],
#'   homogeneous = FALSE
#' )
#' 
#' # clip to sphere of radius 4.8
#' clipper <- sphereMesh(r = 4.8)
#' \donttest{clippedMesh <- clipMesh(mesh, clipper, clipVolume = FALSE, normals = TRUE)
#' 
#' # plot
#' open3d(windowRect = c(50, 50, 950, 500))
#' mfrow3d(1L, 2L)
#' view3d(0, -70, zoom = 0.8)
#' shade3d(mesh, color = "firebrick")
#' next3d()
#' view3d(0, -70, zoom = 0.8)
#' shade3d(toRGL(clippedMesh), color = "firebrick")}
clipMesh <- function(mesh, clipper, clipVolume = TRUE, normals = FALSE){
	stopifnot(isBoolean(clipVolume))
	stopifnot(isBoolean(normals))
	if(inherits(mesh, "mesh3d")){
		vft  <- getVFT(mesh, beforeCheck = TRUE)
		mesh <- vft[["rmesh"]]
	}
	vertices <- mesh[["vertices"]]
	faces    <- mesh[["faces"]]
	checkedMesh <- checkMesh(vertices, faces, gmp = FALSE, aslist = TRUE)
	vertices         <- checkedMesh[["vertices"]]
	faces            <- checkedMesh[["faces"]]
	isTriangle       <- checkedMesh[["isTriangle"]]
	rmesh <- list("vertices" = vertices, "faces" = faces)
	triangulate1 <- !isTriangle
	if(inherits(clipper, "mesh3d")){
		vft  <- getVFT(clipper, beforeCheck = TRUE)
		clipper <- vft[["rmesh"]]
	}
	vertices <- clipper[["vertices"]]
	faces    <- clipper[["faces"]]
	checkedMesh <- checkMesh(vertices, faces, gmp = FALSE, aslist = TRUE)
	vertices         <- checkedMesh[["vertices"]]
	faces            <- checkedMesh[["faces"]]
	isTriangle       <- checkedMesh[["isTriangle"]]
	rclipper <- list("vertices" = vertices, "faces" = faces)
	triangulate2 <- !isTriangle
	mesh <- clipMeshEK(
		rmesh, rclipper, clipVolume, triangulate1, triangulate2, normals
	)
	mesh[["vertices"]] <- t(mesh[["vertices"]])
	mesh[["faces"]] <- t(mesh[["faces"]])
	edgesDF <- mesh[["edges"]]
	mesh[["edgesDF"]] <- edgesDF
	mesh[["edges"]] <- as.matrix(edgesDF[, c("i1", "i2")])
	exteriorEdges <- as.matrix(subset(edgesDF, exterior)[, c("i1", "i2")])
	mesh[["exteriorEdges"]] <- exteriorEdges
	mesh[["exteriorVertices"]] <- which(table(exteriorEdges) != 2L)
	if(normals){
		mesh[["normals"]] <- t(mesh[["normals"]])
	}
	attr(mesh, "toRGL") <- 3L
	class(mesh) <- "cgalMesh"
	mesh
}