Title: | Mesh Generation and Surface Tessellation |
---|---|
Description: | Makes the 'Qhull' library <http://www.qhull.org> available in R, in a similar manner as in Octave and MATLAB. Qhull computes convex hulls, Delaunay triangulations, halfspace intersections about a point, Voronoi diagrams, furthest-site Delaunay triangulations, and furthest-site Voronoi diagrams. It runs in 2D, 3D, 4D, and higher dimensions. It implements the Quickhull algorithm for computing the convex hull. Qhull does not support constrained Delaunay triangulations, or mesh generation of non-convex objects, but the package does include some R functions that allow for this. |
Authors: | Jean-Romain Roussel [cph, ctb] (wrote tsearch function with
QuadTrees), C. B. Barber [cph], Kai Habel [cph, aut], Raoul
Grasman [cph, aut], Robert B. Gramacy [cph, aut], Pavlo
Mozharovskyi [cph, aut], David C. Sterratt [cph, aut, cre]
|
Maintainer: | David C. Sterratt <[email protected]> |
License: | GPL (>= 3) |
Version: | 0.5.0 |
Built: | 2024-09-01 23:20:21 UTC |
Source: | CRAN |
Given the barycentric coordinates of one or more points with respect to a simplex, compute the Cartesian coordinates of these points.
bary2cart(X, Beta)
bary2cart(X, Beta)
X |
Reference simplex in |
Beta |
|
-by-
matrix in which each row is the
Cartesian coordinates of corresponding row of
Beta
David Sterratt
## Define simplex in 2D (i.e. a triangle) X <- rbind(c(0, 0), c(0, 1), c(1, 0)) ## Cartesian cooridinates of points beta <- rbind(c(0, 0.5, 0.5), c(0.1, 0.8, 0.1)) ## Plot triangle and points trimesh(rbind(1:3), X) text(X[,1], X[,2], 1:3) # Label vertices P <- bary2cart(X, beta) points(P)
## Define simplex in 2D (i.e. a triangle) X <- rbind(c(0, 0), c(0, 1), c(1, 0)) ## Cartesian cooridinates of points beta <- rbind(c(0, 0.5, 0.5), c(0.1, 0.8, 0.1)) ## Plot triangle and points trimesh(rbind(1:3), X) text(X[,1], X[,2], 1:3) # Label vertices P <- bary2cart(X, beta) points(P)
Given the Cartesian coordinates of one or more points, compute the barycentric coordinates of these points with respect to a simplex.
cart2bary(X, P)
cart2bary(X, P)
X |
Reference simplex in |
P |
|
Given a reference simplex in dimensions represented by a
-by-
matrix an arbitrary point
in
Cartesian coordinates, represented by a 1-by-
row vector, can be
written as
where is an
vector of the barycentric coordinates.
A criterion on
is that
Now partition the simplex into its first rows
and
its
th row
. Partition the barycentric
coordinates into the first
columns
and the
th column
. This allows us to write
which can be written
where is an
-by-1 matrix of ones. We can then solve
for
:
and compute
This can be generalised for multiple values of
, one per row.
-by-
matrix in which each row is the
barycentric coordinates of corresponding row of
P
. If the
simplex is degenerate a warning is issued and the function returns
NULL
.
Based on the Octave function by David Bateman.
David Sterratt
## Define simplex in 2D (i.e. a triangle) X <- rbind(c(0, 0), c(0, 1), c(1, 0)) ## Cartesian coordinates of points P <- rbind(c(0.5, 0.5), c(0.1, 0.8)) ## Plot triangle and points trimesh(rbind(1:3), X) text(X[,1], X[,2], 1:3) # Label vertices points(P) cart2bary(X, P)
## Define simplex in 2D (i.e. a triangle) X <- rbind(c(0, 0), c(0, 1), c(1, 0)) ## Cartesian coordinates of points P <- rbind(c(0.5, 0.5), c(0.1, 0.8)) ## Plot triangle and points trimesh(rbind(1:3), X) text(X[,1], X[,2], 1:3) # Label vertices points(P) cart2bary(X, P)
The inputs x
, y
(, and z
) must be the same shape, or
scalar. If called with a single matrix argument then each row of C
represents the Cartesian coordinate (x
, y
(, z
)).
cart2pol(x, y = NULL, z = NULL)
cart2pol(x, y = NULL, z = NULL)
x |
x-coordinates or matrix with three columns |
y |
y-coordinates (optional, if |
z |
z-coordinates (optional, if |
A matrix P
where each row represents one
polar/(cylindrical) coordinate (theta
, r
, (,
z
)).
Kai Habel
David Sterratt
If called with a single matrix argument then each row of c
represents the Cartesian coordinate (x
, y
, z
).
cart2sph(x, y = NULL, z = NULL)
cart2sph(x, y = NULL, z = NULL)
x |
x-coordinates or matrix with three columns |
y |
y-coordinates (optional, if |
z |
z-coordinates (optional, if |
Matrix with columns:
theta |
the angle relative to the positive x-axis |
phi |
the angle relative to the xy-plane |
r |
the distance to the origin |
Kai Habel
David Sterratt
Returns information about the smallest convex complex of a set of
input points in -dimensional space (the convex hull of the
points). By default, indices to points forming the facets of the
hull are returned; optionally normals to the facets and the
generalised surface area and volume can be returned. This function
interfaces the Qhull library.
convhulln( p, options = "Tv", output.options = NULL, return.non.triangulated.facets = FALSE )
convhulln( p, options = "Tv", output.options = NULL, return.non.triangulated.facets = FALSE )
p |
An |
options |
String containing extra options for the underlying Qhull command; see details below and Qhull documentation at ../doc/qhull/html/qconvex.html#synopsis. |
output.options |
String containing Qhull options to generate
extra output. Currently |
return.non.triangulated.facets |
logical defining whether the
output facets should be triangulated; |
By default (return.non.triangulated.facets
is
FALSE
), return an -by-
matrix in which each
row contains the indices of the points in
p
forming an
-dimensional facet. e.g In 3 dimensions, there are 3
indices in each row describing the vertices of 2-dimensional
triangles.
If return.non.triangulated.facets
is TRUE
then the
number of columns equals the maximum number of vertices in a
facet, and each row defines a polygon corresponding to a facet
of the convex hull with its vertices followed by NA
s
until the end of the row.
If the output.options
or options
argument contains
FA
or n
, return a list with class convhulln
comprising the named elements:
p
The points passed to convnhulln
hull
The convex hull, represented as a matrix indexing p
, as
described above
area
If FA
is specified, the generalised area of
the hull. This is the surface area of a 3D hull or the length of
the perimeter of a 2D hull.
See ../doc/qhull/html/qh-optf.html#FA.
vol
If FA
is specified, the generalised volume of
the hull. This is volume of a 3D hull or the area of a 2D hull.
See ../doc/qhull/html/qh-optf.html#FA.
normals
If n
is specified, this is a matrix
hyperplane normals with offsets. See ../doc/qhull/html/qh-opto.html#n.
This function was originally a port of the Octave convhulln function written by Kai Habel.
See further notes in delaunayn
.
Raoul Grasman, Robert B. Gramacy, Pavlo Mozharovskyi and David Sterratt [email protected]
Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T., “The Quickhull algorithm for convex hulls,” ACM Trans. on Mathematical Software, Dec 1996.
intersectn
, delaunayn
,
surf.tri
, convex.hull
## Points in a sphere ps <- matrix(rnorm(3000), ncol=3) ps <- sqrt(3)*ps/drop(sqrt((ps^2) %*% rep(1, 3))) ts.surf <- t(convhulln(ps)) # see the qhull documentations for the options ## Not run: rgl::triangles3d(ps[ts.surf,1],ps[ts.surf,2],ps[ts.surf,3],col="blue",alpha=.2) for(i in 1:(8*360)) rgl::view3d(i/8) ## End(Not run) ## Square pq <- rbox(0, C=0.5, D=2) # Return indices only convhulln(pq) # Return convhulln object with normals, generalised area and volume ch <- convhulln(pq, output.options=TRUE) plot(ch) ## Cube pc <- rbox(0, C=0.5, D=3) # Return indices of triangles on surface convhulln(pc) # Return indices of squares on surface convhulln(pc, return.non.triangulated.facets=TRUE)
## Points in a sphere ps <- matrix(rnorm(3000), ncol=3) ps <- sqrt(3)*ps/drop(sqrt((ps^2) %*% rep(1, 3))) ts.surf <- t(convhulln(ps)) # see the qhull documentations for the options ## Not run: rgl::triangles3d(ps[ts.surf,1],ps[ts.surf,2],ps[ts.surf,3],col="blue",alpha=.2) for(i in 1:(8*360)) rgl::view3d(i/8) ## End(Not run) ## Square pq <- rbox(0, C=0.5, D=2) # Return indices only convhulln(pq) # Return convhulln object with normals, generalised area and volume ch <- convhulln(pq, output.options=TRUE) plot(ch) ## Cube pc <- rbox(0, C=0.5, D=3) # Return indices of triangles on surface convhulln(pc) # Return indices of squares on surface convhulln(pc, return.non.triangulated.facets=TRUE)
The Delaunay triangulation is a tessellation of the convex hull of
the points such that no -sphere defined by the
-
triangles contains any other points from the set.
delaunayn(p, options = NULL, output.options = NULL, full = FALSE)
delaunayn(p, options = NULL, output.options = NULL, full = FALSE)
p |
An |
options |
String containing extra control options for the underlying Qhull command; see the Qhull documentation (../doc/qhull/html/qdelaun.html) for the available options. The If |
output.options |
String containing Qhull options to control
output. Currently |
full |
Deprecated and will be removed in a future release.
Adds options |
If output.options
is NULL
(the default),
return the Delaunay triangulation as a matrix with rows
and
columns in which each row contains a set of
indices to the input points
p
. Thus each row describes a
simplex of dimension , e.g. a triangle in 2D or a
tetrahedron in 3D.
If the output.options
argument is TRUE
or is a
string containing Fn
or Fa
, return a list with
class delaunayn
comprising the named elements:
tri
The Delaunay triangulation described above
areas
If TRUE
or if Fa
is specified, an
-dimensional vector containing the generalised area of
each simplex (e.g. in 2D the areas of triangles; in 3D the volumes
of tetrahedra). See ../doc/qhull/html/qh-optf.html#Fa.
neighbours
If TRUE
or if Fn
is specified,
a list of neighbours of each simplex. Note that a negative number
corresponds to "facet" (="edge" in 2D or "face" in 3D) that has no
neighbour, as will be the case for some simplices on the boundary
of the triangulation.
See ../doc/qhull/html/qh-optf.html#Fn
This function interfaces the Qhull library and is a port from Octave (https://octave.org/) to R. Qhull computes convex hulls, Delaunay triangulations, halfspace intersections about a point, Voronoi diagrams, furthest-site Delaunay triangulations, and furthest-site Voronoi diagrams. It runs in 2D, 3D, 4D, and higher dimensions. It implements the Quickhull algorithm for computing the convex hull. Qhull handles round-off errors from floating point arithmetic. It computes volumes, surface areas, and approximations to the convex hull. See the Qhull documentation included in this distribution (the doc directory ../doc/qhull/index.html).
Qhull does not support constrained Delaunay triangulations, triangulation of non-convex surfaces, mesh generation of non-convex objects, or medium-sized inputs in 9D and higher. A rudimentary algorithm for mesh generation in non-convex regions using Delaunay triangulation is implemented in distmesh2d (currently only 2D).
Raoul Grasman and Robert B. Gramacy; based on the corresponding Octave sources of Kai Habel.
Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T., “The Quickhull algorithm for convex hulls,” ACM Trans. on Mathematical Software, Dec 1996.
tri.mesh
, convhulln
,
surf.tri
, distmesh2d
# example delaunayn d <- c(-1,1) pc <- as.matrix(rbind(expand.grid(d,d,d),0)) tc <- delaunayn(pc) # example tetramesh ## Not run: rgl::view3d(60) rgl::light3d(120,60) tetramesh(tc,pc, alpha=0.9) ## End(Not run) tc1 <- delaunayn(pc, output.options="Fa") ## sum of generalised areas is total volume of cube sum(tc1$areas)
# example delaunayn d <- c(-1,1) pc <- as.matrix(rbind(expand.grid(d,d,d),0)) tc <- delaunayn(pc) # example tetramesh ## Not run: rgl::view3d(60) rgl::light3d(120,60) tetramesh(tc,pc, alpha=0.9) ## End(Not run) tc1 <- delaunayn(pc, output.options="Fa") ## sum of generalised areas is total volume of cube sum(tc1$areas)
An unstructured simplex requires a choice of mesh points (vertex nodes) and
a triangulation. This is a simple and short algorithm that improves the
quality of a mesh by relocating the mesh points according to a relaxation
scheme of forces in a truss structure. The topology of the truss is reset
using Delaunay triangulation. A (sufficiently smooth) user supplied signed
distance function (fd
) indicates if a given node is inside or
outside the region. Points outside the region are projected back to the
boundary.
distmesh2d( fd, fh, h0, bbox, p = NULL, pfix = array(0, dim = c(0, 2)), ..., dptol = 0.001, ttol = 0.1, Fscale = 1.2, deltat = 0.2, geps = 0.001 * h0, deps = sqrt(.Machine$double.eps) * h0, maxiter = 1000, plot = TRUE )
distmesh2d( fd, fh, h0, bbox, p = NULL, pfix = array(0, dim = c(0, 2)), ..., dptol = 0.001, ttol = 0.1, Fscale = 1.2, deltat = 0.2, geps = 0.001 * h0, deps = sqrt(.Machine$double.eps) * h0, maxiter = 1000, plot = TRUE )
fd |
Vectorized signed distance function, for example
|
fh |
Vectorized function, for example
|
h0 |
Initial distance between mesh nodes. (Ignored of
|
bbox |
Bounding box |
p |
An |
pfix |
|
... |
parameters to be passed to |
dptol |
Algorithm stops when all node movements are smaller
than |
ttol |
Controls how far the points can move (relatively)
before a retriangulation with |
Fscale |
“Internal pressure” in the edges. |
deltat |
Size of the time step in Euler's method. |
geps |
Tolerance in the geometry evaluations. |
deps |
Stepsize |
maxiter |
Maximum iterations. |
plot |
logical. If |
This is an implementation of original Matlab software of Per-Olof Persson.
Excerpt (modified) from the reference below:
‘The algorithm is based on a mechanical analogy between a triangular
mesh and a 2D truss structure. In the physical model, the edges of the
Delaunay triangles of a set of points correspond to bars of a truss. Each
bar has a force-displacement relationship
depending on its current length
and its unextended length
.’
‘External forces on the structure come at the boundaries, on which
external forces have normal orientations. These external forces are just
large enough to prevent nodes from moving outside the boundary. The
position of the nodes are the unknowns, and are found by solving for a
static force equilibrium. The hope is that (when fh = function(p)
return(rep(1,nrow(p)))
), the lengths of all the bars at equilibrium will
be nearly equal, giving a well-shaped triangular mesh.’
See the references below for all details. Also, see the comments in the source file.
n
-by-2
matrix with node positions.
Implement in C/Fortran
Implement an n
D version as provided in the Matlab
package
Translate other functions of the Matlab package
Raoul Grasman
http://persson.berkeley.edu/distmesh/
P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM Review, Volume 46 (2), pp. 329-345, June 2004
tri.mesh
, delaunayn
,
mesh.dcircle
, mesh.drectangle
,
mesh.diff
, mesh.union
,
mesh.intersect
# examples distmesh2d fd <- function(p, ...) sqrt((p^2)%*%c(1,1)) - 1 # also predefined as `mesh.dcircle' fh <- function(p,...) rep(1,nrow(p)) bbox <- matrix(c(-1,1,-1,1),2,2) p <- distmesh2d(fd,fh,0.2,bbox, maxiter=100) # this may take a while: # press Esc to get result of current iteration # example with non-convex region fd <- function(p, ...) mesh.diff(p , mesh.drectangle, mesh.dcircle, radius=.3) # fd defines difference of square and circle p <- distmesh2d(fd,fh,0.05,bbox,radius=0.3,maxiter=4) p <- distmesh2d(fd,fh,0.05,bbox,radius=0.3, maxiter=10) # continue on previous mesh
# examples distmesh2d fd <- function(p, ...) sqrt((p^2)%*%c(1,1)) - 1 # also predefined as `mesh.dcircle' fh <- function(p,...) rep(1,nrow(p)) bbox <- matrix(c(-1,1,-1,1),2,2) p <- distmesh2d(fd,fh,0.2,bbox, maxiter=100) # this may take a while: # press Esc to get result of current iteration # example with non-convex region fd <- function(p, ...) mesh.diff(p , mesh.drectangle, mesh.dcircle, radius=.3) # fd defines difference of square and circle p <- distmesh2d(fd,fh,0.05,bbox,radius=0.3,maxiter=4) p <- distmesh2d(fd,fh,0.05,bbox,radius=0.3, maxiter=10) # continue on previous mesh
An unstructured simplex requires a choice of mesh points (vertex nodes) and
a triangulation. This is a simple and short algorithm that improves the
quality of a mesh by relocating the mesh points according to a relaxation
scheme of forces in a truss structure. The topology of the truss is reset
using Delaunay triangulation. A (sufficiently smooth) user supplied signed
distance function (fd
) indicates if a given node is inside or
outside the region. Points outside the region are projected back to the
boundary.
distmeshnd( fdist, fh, h, box, pfix = array(dim = c(0, ncol(box))), ..., ptol = 0.001, ttol = 0.1, deltat = 0.1, geps = 0.1 * h, deps = sqrt(.Machine$double.eps) * h )
distmeshnd( fdist, fh, h, box, pfix = array(dim = c(0, ncol(box))), ..., ptol = 0.001, ttol = 0.1, deltat = 0.1, geps = 0.1 * h, deps = sqrt(.Machine$double.eps) * h )
fdist |
Vectorized signed distance function, for example
|
fh |
Vectorized function, for example |
h |
Initial distance between mesh nodes. |
box |
|
pfix |
|
... |
parameters that are passed to |
ptol |
Algorithm stops when all node movements are smaller than
|
ttol |
Controls how far the points can move (relatively) before a
retriangulation with |
deltat |
Size of the time step in Euler's method. |
geps |
Tolerance in the geometry evaluations. |
deps |
Stepsize |
This is an implementation of original Matlab software of Per-Olof Persson.
Excerpt (modified) from the reference below:
‘The algorithm is based on a mechanical analogy between a triangular
mesh and a n-D truss structure. In the physical model, the edges of the
Delaunay triangles of a set of points correspond to bars of a truss. Each
bar has a force-displacement relationship
depending on its current length
and its unextended length
.’
‘External forces on the structure come at the boundaries, on which
external forces have normal orientations. These external forces are just
large enough to prevent nodes from moving outside the boundary. The
position of the nodes are the unknowns, and are found by solving for a
static force equilibrium. The hope is that (when fh = function(p)
return(rep(1,nrow(p)))
), the lengths of all the bars at equilibrium will
be nearly equal, giving a well-shaped triangular mesh.’
See the references below for all details. Also, see the comments in the
source file of distmesh2d
.
m
-by-n
matrix with node positions.
Implement in C/Fortran
Translate other functions of the Matlab package
Raoul Grasman; translated from original Matlab sources of Per-Olof Persson.
http://persson.berkeley.edu/distmesh/
P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM Review, Volume 46 (2), pp. 329-345, June 2004
distmesh2d
, tri.mesh
,
delaunayn
, mesh.dsphere
,
mesh.hunif
,mesh.diff
,
mesh.union
, mesh.intersect
## Not run: # examples distmeshnd require(rgl) fd = function(p, ...) sqrt((p^2)%*%c(1,1,1)) - 1 # also predefined as `mesh.dsphere' fh = function(p,...) rep(1,nrow(p)) # also predefined as `mesh.hunif' bbox = matrix(c(-1,1),2,3) p = distmeshnd(fd,fh,0.2,bbox, maxiter=100) # this may take a while: # press Esc to get result of current iteration ## End(Not run)
## Not run: # examples distmeshnd require(rgl) fd = function(p, ...) sqrt((p^2)%*%c(1,1,1)) - 1 # also predefined as `mesh.dsphere' fh = function(p,...) rep(1,nrow(p)) # also predefined as `mesh.hunif' bbox = matrix(c(-1,1),2,3) p = distmeshnd(fd,fh,0.2,bbox, maxiter=100) # this may take a while: # press Esc to get result of current iteration ## End(Not run)
If x
and y
are matrices, calculate the dot-product
along the first non-singleton dimension. If the optional argument
d
is given, calculate the dot-product along this
dimension.
dot(x, y, d = NULL)
dot(x, y, d = NULL)
x |
Matrix of vectors |
y |
Matrix of vectors |
d |
Dimension along which to calculate the dot product |
Vector with length of d
th dimension
David Sterratt
entry.value
retrieves or sets the values in an array a
at the
positions indicated by the rows of a matrix idx
.
entry.value(a, idx)
entry.value(a, idx)
a |
An array. |
idx |
Numerical matrix with the same number of columns as the number
of dimensions of |
value |
An array of length |
entry.value(a,idx)
returns a vector of values at the
indicated cells. entry.value(a,idx) <- val
changes the indicated
cells of a
to val
.
Raoul Grasman
a = array(1:(4^4),c(4,4,4,4)) entry.value(a,cbind(1:4,1:4,1:4,1:4)) entry.value(a,cbind(1:4,1:4,1:4,1:4)) <- 0 entry.value(a, as.matrix(expand.grid(1:4,1:4,1:4,1:4))) # same as `c(a[1:4,1:4,1:4,1:4])' which is same as `c(a)'
a = array(1:(4^4),c(4,4,4,4)) entry.value(a,cbind(1:4,1:4,1:4,1:4)) entry.value(a,cbind(1:4,1:4,1:4,1:4)) <- 0 entry.value(a, as.matrix(expand.grid(1:4,1:4,1:4,1:4))) # same as `c(a[1:4,1:4,1:4,1:4])' which is same as `c(a)'
Computes the external product
of the 3D vectors in x and y.
extprod3d(x, y, drop = TRUE)
extprod3d(x, y, drop = TRUE)
x |
|
y |
|
drop |
logical. If |
If n
is greater than 1 or drop
is
FALSE
, n
-by-3 matrix; if n
is 1 and
drop
is TRUE
, a vector of length 3.
Raoul Grasman
Find point that lies somewhere in interesction of two convex
hulls. If such a point does not exist, return NA
. The
feasible point is found using a linear program similar to the one
suggested at ../doc/qhull/html/qhalf.html#notes
feasible.point(ch1, ch2, tol = 0)
feasible.point(ch1, ch2, tol = 0)
ch1 |
First convex hull with normals |
ch2 |
Second convex hull with normals |
tol |
The point must be at least this far within the facets of both convex hulls |
Compute halfspace intersection about a point
halfspacen(p, fp, options = "Tv")
halfspacen(p, fp, options = "Tv")
p |
An |
fp |
A “feasible” point that is within the space contained within all the halfspaces. |
options |
String containing extra options, separated by spaces, for the underlying Qhull command; see Qhull documentation at ../doc/qhull/html/qhalf.html. |
A -column matrix containing the intersection
points of the hyperplanes ../doc/qhull/html/qhalf.html.
halfspacen
was introduced in geometry 0.4.0, and is
still under development. It is worth checking results for
unexpected behaviour.
David Sterratt
Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T., “The Quickhull algorithm for convex hulls,” ACM Trans. on Mathematical Software, Dec 1996.
p <- rbox(0, C=0.5) # Generate points on a unit cube centered around the origin ch <- convhulln(p, "n") # Generate convex hull, including normals to facets, with "n" option # Intersections of half planes # These points should be the same as the orginal points pn <- halfspacen(ch$normals, c(0, 0, 0))
p <- rbox(0, C=0.5) # Generate points on a unit cube centered around the origin ch <- convhulln(p, "n") # Generate convex hull, including normals to facets, with "n" option # Intersections of half planes # These points should be the same as the orginal points pn <- halfspacen(ch$normals, c(0, 0, 0))
Tests if a set of points lies within a convex hull, returning a
boolean vector in which each element is TRUE
if the
corresponding point lies within the hull and FALSE
if it
lies outwith the hull or on one of its facets.
inhulln(ch, p)
inhulln(ch, p)
ch |
Convex hull produced using |
p |
An |
A boolean vector with elements
inhulln
was introduced in geometry 0.4.0, and is
still under development. It is worth checking results for
unexpected behaviour.
David Sterratt
convhulln
, point.in.polygon
in sp
p <- cbind(c(-1, -1, 1), c(-1, 1, -1)) ch <- convhulln(p) ## First point should be in the hull; last two outside inhulln(ch, rbind(c(-0.5, -0.5), c( 1 , 1), c(10 , 0))) ## Test hypercube p <- rbox(D=4, B=1) ch <- convhulln(p) tp <- cbind(seq(-1.9, 1.9, by=0.2), 0, 0, 0) pin <- inhulln(ch, tp) ## Points on x-axis should be in box only betw,een -1 and 1 pin == (tp[,1] < 1 & tp[,1] > -1)
p <- cbind(c(-1, -1, 1), c(-1, 1, -1)) ch <- convhulln(p) ## First point should be in the hull; last two outside inhulln(ch, rbind(c(-0.5, -0.5), c( 1 , 1), c(10 , 0))) ## Test hypercube p <- rbox(D=4, B=1) ch <- convhulln(p) tp <- cbind(seq(-1.9, 1.9, by=0.2), 0, 0, 0) pin <- inhulln(ch, tp) ## Points on x-axis should be in box only betw,een -1 and 1 pin == (tp[,1] < 1 & tp[,1] > -1)
Compute convex hull of intersection of two sets of points
intersectn( ps1, ps2, tol = 0, return.chs = TRUE, options = "Tv", fp = NULL, autoscale = FALSE )
intersectn( ps1, ps2, tol = 0, return.chs = TRUE, options = "Tv", fp = NULL, autoscale = FALSE )
ps1 |
First set of points |
ps2 |
Second set of points |
tol |
Tolerance used to determine if a feasible point lies within the convex hulls of both points and to round off the points generated by the halfspace intersection, which sometimes produces points very close together. |
return.chs |
If |
options |
Options passed to |
fp |
Coordinates of feasible point, i.e. a point known to lie
in the hulls of |
autoscale |
Experimental in v0.4.2 Automatically scale the points to lie in a sensible numeric range. May help to correct some numerical issues. |
List containing named elements: ch1
, the convex
hull of the first set of points, with volumes, areas and normals
(see convhulln
; ch2
, the convex hull of the
first set of points, with volumes, areas and normals; ps
,
the intersection points of convex hulls ch1
and
ch2
; and ch
, the convex hull of the intersection
points, with volumes, areas and normals.
intersectn
was introduced in geometry 0.4.0, and is
still under development. It is worth checking results for
unexpected behaviour.
David Sterratt
convhulln
, halfspacen
,
inhulln
, feasible.point
# Two overlapping boxes ps1 <- rbox(0, C=0.5) ps2 <- rbox(0, C=0.5) + 0.5 out <- intersectn(ps1, ps2) message("Volume of 1st convex hull: ", out$ch1$vol) message("Volume of 2nd convex hull: ", out$ch2$vol) message("Volume of intersection convex hull: ", out$ch$vol)
# Two overlapping boxes ps1 <- rbox(0, C=0.5) ps2 <- rbox(0, C=0.5) + 0.5 out <- intersectn(ps1, ps2) message("Volume of 1st convex hull: ", out$ch1$vol) message("Volume of 2nd convex hull: ", out$ch2$vol) message("Volume of intersection convex hull: ", out$ch$vol)
Compute maximum or minimum of each row, or sort each row of a matrix, or a set of (equal length) vectors.
matmax(...)
matmax(...)
... |
A numeric matrix or a set of numeric vectors (that are
column-wise bind together into a matrix with |
matmin
and matmax
return a vector of length
nrow(cbind(...))
. matsort
returns a matrix of dimension
dim(cbind(...))
with in each row of cbind(...)
sorted.
matsort(x)
is a lot faster than, e.g., t(apply(x,1,sort))
,
if x
is tall (i.e., nrow(x)
>>ncol(x)
and
ncol(x)
<30. If ncol(x)
>30 then matsort
simply calls
't(apply(x,1,sort))
'. matorder
returns a permutation which
rearranges its first argument into ascending order, breaking ties by
further arguments.
Raoul Grasman
example(Unique)
example(Unique)
Signed distance from points p
to boundary of circle to
allow easy definition of regions in distmesh2d
.
mesh.dcircle(p, radius = 1, ...)
mesh.dcircle(p, radius = 1, ...)
p |
A matrix with 2 columns (3 in |
radius |
radius of circle |
... |
additional arguments (not used) |
A vector of length nrow(p)
containing the signed
distances to the circle
Raoul Grasman; translated from original Matlab sources of Per-Olof Persson.
http://persson.berkeley.edu/distmesh/
P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM Review, Volume 46 (2), pp. 329-345, June 2004
distmesh2d
, mesh.drectangle
,
mesh.diff
, mesh.intersect
,
mesh.union
example(distmesh2d)
example(distmesh2d)
Compute the signed distances from points p
to a region
defined by the difference, union or intersection of regions
specified by the functions regionA
and regionB
.
regionA
and regionB
must accept a matrix p
with 2 columns as their first argument, and must return a vector
of length nrow(p)
containing the signed distances of the
supplied points in p
to their respective regions.
mesh.diff(p, regionA, regionB, ...)
mesh.diff(p, regionA, regionB, ...)
p |
A matrix with 2 columns (3 in |
regionA |
vectorized function describing region A in the union / intersection / difference |
regionB |
vectorized function describing region B in the union / intersection / difference |
... |
additional arguments passed to |
A vector of length nrow(p)
containing the signed
distances to the boundary of the region.
Raoul Grasman; translated from original Matlab sources of Per-Olof Persson.
distmesh2d
, mesh.dcircle
,
mesh.drectangle
mesh.dsphere
Signed distance from points p
to boundary of rectangle to
allow easy definition of regions in distmesh2d
.
mesh.drectangle(p, x1 = -1/2, y1 = -1/2, x2 = 1/2, y2 = 1/2, ...)
mesh.drectangle(p, x1 = -1/2, y1 = -1/2, x2 = 1/2, y2 = 1/2, ...)
p |
A matrix with 2 columns, each row representing a point in the plane. |
x1 |
lower left corner of rectangle |
y1 |
lower left corner of rectangle |
x2 |
upper right corner of rectangle |
y2 |
upper right corner of rectangle |
... |
additional arguments (not used) |
a vector of length nrow(p)
containing the signed distances
Raoul Grasman; translated from original Matlab sources of Per-Olof Persson.
http://persson.berkeley.edu/distmesh/
P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM Review, Volume 46 (2), pp. 329-345, June 2004
distmesh2d
, mesh.drectangle
,
mesh.diff
, mesh.intersect
,
mesh.union
example(distmesh2d)
example(distmesh2d)
Signed distance from points p
to boundary of sphere to
allow easy definition of regions in distmeshnd
.
mesh.dsphere(p, radius = 1, ...)
mesh.dsphere(p, radius = 1, ...)
p |
A matrix with 2 columns (3 in |
radius |
radius of sphere |
... |
additional arguments (not used) |
A vector of length nrow(p)
containing the signed
distances to the sphere
Raoul Grasman; translated from original Matlab sources of Per-Olof Persson.
http://persson.berkeley.edu/distmesh/
P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM Review, Volume 46 (2), pp. 329-345, June 2004
example(distmeshnd)
example(distmeshnd)
Uniform desired edge length function of position to allow easy
definition of regions when passed as the fh
argument of
distmesh2d
or distmeshnd
.
mesh.hunif(p, ...)
mesh.hunif(p, ...)
p |
A |
... |
additional arguments (not used) |
Vector of ones of length n
.
Raoul Grasman; translated from original Matlab sources of Per-Olof Persson.
distmesh2d
and distmeshnd
.
The inputs theta
, r
, (and z
) must be the same shape, or
scalar. If called with a single matrix argument then each row of P
represents the polar/(cylindrical) coordinate (theta
, r
(, z
)).
pol2cart(theta, r = NULL, z = NULL)
pol2cart(theta, r = NULL, z = NULL)
theta |
describes the angle relative to the positive x-axis. |
r |
is the distance to the z-axis (0, 0, z). |
z |
(optional) is the z-coordinate |
a matrix C
where each row represents one Cartesian
coordinate (x
, y
(, z
)).
Kai Habel
David Sterratt
Determines area of a polygon by triangle method. The variables
x
and y
define the vertex pairs, and must therefore
have the same shape. They can be either vectors or arrays. If
they are arrays then the columns of x
and y
are
treated separately and an area returned for each.
polyarea(x, y, d = 1)
polyarea(x, y, d = 1)
x |
X coordinates of vertices. |
y |
Y coordinates of vertices. |
d |
Dimension of array to work along. |
If the optional dim
argument is given, then polyarea
works along this dimension of the arrays x
and y
.
Area(s) of polygon(s).
David Sterratt based on the octave sources by David M. Doolin
x <- c(1, 1, 3, 3, 1) y <- c(1, 3, 3, 1, 1) polyarea(x, y) polyarea(cbind(x, x), cbind(y, y)) ## c(4, 4) polyarea(cbind(x, x), cbind(y, y), 1) ## c(4, 4) polyarea(rbind(x, x), rbind(y, y), 2) ## c(4, 4)
x <- c(1, 1, 3, 3, 1) y <- c(1, 3, 3, 1, 1) polyarea(x, y) polyarea(cbind(x, x), cbind(y, y)) ## c(4, 4) polyarea(cbind(x, x), cbind(y, y), 1) ## c(4, 4) polyarea(rbind(x, x), rbind(y, y), 2) ## c(4, 4)
Default is corners of a hypercube.
rbox(n = 3000, D = 3, B = 0.5, C = NA)
rbox(n = 3000, D = 3, B = 0.5, C = NA)
n |
number of random points in hypercube |
D |
number of dimensions of hypercube |
B |
bounding box coordinate - faces will be |
C |
add a unit hypercube to the output - faces will be |
Matrix of points
David Sterratt
The inputs theta
, phi
, and r
must be the same
shape, or scalar. If called with a single matrix argument then
each row of S
represents the spherical coordinate
(theta
, phi
, r
).
sph2cart(theta, phi = NULL, r = NULL)
sph2cart(theta, phi = NULL, r = NULL)
theta |
describes the angle relative to the positive x-axis. |
phi |
is the angle relative to the xy-plane. |
r |
is the distance to the origin If only a single return argument is requested then return a matrix
|
Kai Habel
David Sterratt
Find surface triangles from tetrahedral mesh typically obtained
with delaunayn
.
surf.tri(p, t)
surf.tri(p, t)
p |
An |
t |
Matrix with 4 columns, interpreted as output of
|
surf.tri
and convhulln
serve a similar purpose in 3D,
but surf.tri
also works for non-convex meshes obtained e.g. with
distmeshnd
. It also does not produce currently unavoidable
diagnostic output on the console as convhulln
does at the Rterm
console–i.e., surf.tri
is silent.
An m
-by-3
index matrix of which each row defines a
triangle. The indices refer to the rows in p
.
surf.tri
was based on Matlab code for mesh of Per-Olof Persson
(http://persson.berkeley.edu/distmesh/).
Raoul Grasman
tri.mesh
, convhulln
,
surf.tri
, distmesh2d
## Not run: # more extensive example of surf.tri # url's of publically available data: data1.url = "http://neuroimage.usc.edu/USCPhantom/mesh_data.bin" data2.url = "http://neuroimage.usc.edu/USCPhantom/CT_PCS_trans.bin" meshdata = R.matlab::readMat(url(data1.url)) elec = R.matlab::readMat(url(data2.url))$eeg.ct2pcs/1000 brain = meshdata$mesh.brain[,c(1,3,2)] scalp = meshdata$mesh.scalp[,c(1,3,2)] skull = meshdata$mesh.skull[,c(1,3,2)] tbr = t(surf.tri(brain, delaunayn(brain))) tsk = t(surf.tri(skull, delaunayn(skull))) tsc = t(surf.tri(scalp, delaunayn(scalp))) rgl::triangles3d(brain[tbr,1], brain[tbr,2], brain[tbr,3],col="gray") rgl::triangles3d(skull[tsk,1], skull[tsk,2], skull[tsk,3],col="white", alpha=0.3) rgl::triangles3d(scalp[tsc,1], scalp[tsc,2], scalp[tsc,3],col="#a53900", alpha=0.6) rgl::view3d(-40,30,.4,zoom=.03) lx = c(-.025,.025); ly = -c(.02,.02); rgl::spheres3d(elec[,1],elec[,3],elec[,2],radius=.0025,col='gray') rgl::spheres3d( lx, ly,.11,radius=.015,col="white") rgl::spheres3d( lx, ly,.116,radius=.015*.7,col="brown") rgl::spheres3d( lx, ly,.124,radius=.015*.25,col="black") ## End(Not run)
## Not run: # more extensive example of surf.tri # url's of publically available data: data1.url = "http://neuroimage.usc.edu/USCPhantom/mesh_data.bin" data2.url = "http://neuroimage.usc.edu/USCPhantom/CT_PCS_trans.bin" meshdata = R.matlab::readMat(url(data1.url)) elec = R.matlab::readMat(url(data2.url))$eeg.ct2pcs/1000 brain = meshdata$mesh.brain[,c(1,3,2)] scalp = meshdata$mesh.scalp[,c(1,3,2)] skull = meshdata$mesh.skull[,c(1,3,2)] tbr = t(surf.tri(brain, delaunayn(brain))) tsk = t(surf.tri(skull, delaunayn(skull))) tsc = t(surf.tri(scalp, delaunayn(scalp))) rgl::triangles3d(brain[tbr,1], brain[tbr,2], brain[tbr,3],col="gray") rgl::triangles3d(skull[tsk,1], skull[tsk,2], skull[tsk,3],col="white", alpha=0.3) rgl::triangles3d(scalp[tsc,1], scalp[tsc,2], scalp[tsc,3],col="#a53900", alpha=0.6) rgl::view3d(-40,30,.4,zoom=.03) lx = c(-.025,.025); ly = -c(.02,.02); rgl::spheres3d(elec[,1],elec[,3],elec[,2],radius=.0025,col='gray') rgl::spheres3d( lx, ly,.11,radius=.015,col="white") rgl::spheres3d( lx, ly,.116,radius=.015*.7,col="brown") rgl::spheres3d( lx, ly,.124,radius=.015*.25,col="black") ## End(Not run)
tetramesh(T, X, col)
uses the rgl package to
display the tetrahedrons defined in the m-by-4 matrix T as mesh.
Each row of T
specifies a tetrahedron by giving the 4
indices of its points in X
.
tetramesh(T, X, col = grDevices::heat.colors(nrow(T)), clear = TRUE, ...)
tetramesh(T, X, col = grDevices::heat.colors(nrow(T)), clear = TRUE, ...)
T |
T is a |
X |
X is an n-by-2/n-by-3 matrix. The rows of X represent |
col |
The tetrahedron colour. See rgl documentation for details. |
clear |
Should the current rendering device be cleared? |
... |
Parameters to the rendering device. See the rgl package. |
Raoul Grasman
trimesh
, rgl
, delaunayn
,
convhulln
, surf.tri
## Not run: # example delaunayn d = c(-1,1) pc = as.matrix(rbind(expand.grid(d,d,d),0)) tc = delaunayn(pc) # example tetramesh clr = rep(1,3) %o% (1:nrow(tc)+1) rgl::view3d(60,fov=20) rgl::light3d(270,60) tetramesh(tc,pc,alpha=0.7,col=clr) ## End(Not run)
## Not run: # example delaunayn d = c(-1,1) pc = as.matrix(rbind(expand.grid(d,d,d),0)) tc = delaunayn(pc) # example tetramesh clr = rep(1,3) %o% (1:nrow(tc)+1) rgl::view3d(60,fov=20) rgl::light3d(270,60) tetramesh(tc,pc,alpha=0.7,col=clr) ## End(Not run)
trimesh(T, p)
displays the triangles defined in the m-by-3
matrix T
and points p
as a mesh. Each row of
T
specifies a triangle by giving the 3 indices of its
points in X
.
trimesh(T, p, p2, add = FALSE, axis = FALSE, boxed = FALSE, ...)
trimesh(T, p, p2, add = FALSE, axis = FALSE, boxed = FALSE, ...)
T |
T is a |
p |
A vector or a matrix. |
p2 |
if |
add |
Add to existing plot in current active device? |
axis |
Draw axes? |
boxed |
Plot box? |
... |
Parameters to the rendering device. See the rgl package. |
Raoul Grasman
tetramesh
, rgl
,
delaunayn
, convhulln
,
surf.tri
#example trimesh p = cbind(x=rnorm(30), y=rnorm(30)) tt = delaunayn(p) trimesh(tt,p)
#example trimesh p = cbind(x=rnorm(30), y=rnorm(30)) tt = delaunayn(p) trimesh(tt,p)
For t <- delaunay(cbind(x, y))
, where (x, y)
is a 2D set of
points, tsearch(x, y, t, xi, yi)
finds the index in t
containing the points (xi, yi)
. For points outside the convex hull
the index is NA
.
tsearch(x, y, t, xi, yi, bary = FALSE, method = "quadtree")
tsearch(x, y, t, xi, yi, bary = FALSE, method = "quadtree")
x |
X-coordinates of triangulation points |
y |
Y-coordinates of triangulation points |
t |
Triangulation, e.g. produced by |
xi |
X-coordinates of points to test |
yi |
Y-coordinates of points to test |
bary |
If |
method |
One of |
If bary
is FALSE
, the index in t
containing the points
(xi, yi)
. For points outside the convex hull the index is NA
.
If bary
is TRUE
, a list containing:
the index in t
containing the points (xi, yi)
a 3-column matrix containing the barycentric coordinates with
respect to the enclosing triangle of each point (xi, yi)
.
The original Octave function is Copyright (C) 2007-2012 David Bateman
Jean-Romain Roussel (Quadtree algorithm), David Sterratt (Octave-based implementation)
For t = delaunayn(x)
, where x
is a set of points in
dimensions,
tsearchn(x, t, xi)
finds the index in t
containing the points xi
. For points outside the convex hull,
idx
is NA
. tsearchn
also returns the barycentric
coordinates p
of the enclosing triangles.
tsearchn(x, t, xi, ...)
tsearchn(x, t, xi, ...)
x |
An |
t |
A matrix with |
xi |
An |
... |
Additional arguments |
If x
is NA
and the t
is a
delaunayn
object produced by
delaunayn
with the full
option, then use the
Qhull library to perform the search. Please note that this is
experimental in geometry version 0.4.0 and is only partly tested
for 3D hulls, and does not yet work for hulls of 4 dimensions and
above.
A list containing:
idx
An -long vector containing the indices
of the row of
t
in which each point in xi
is found.
p
An -by-
matrix containing the
barycentric coordinates with respect to the enclosing simplex
of each point in
xi
.
Based on the Octave function Copyright (C) 2007-2012 David Bateman.
David Sterratt
‘Unique’ returns a vector, data frame or array like 'x' but with duplicate elements removed.
Unique(X, rows.are.sets = FALSE)
Unique(X, rows.are.sets = FALSE)
X |
Numerical matrix. |
rows.are.sets |
If ‘ |
Matrix of the same number of columns as x
, with the unique
rows in x
sorted according to the columns of x
. If
rows.are.sets = TRUE
the rows are also sorted.
‘Unique
’ is (under circumstances) much quicker than the
more generic base function ‘unique
’.
Raoul Grasman
# `Unique' is faster than `unique' x = matrix(sample(1:(4*8),4*8),ncol=4) y = x[sample(1:nrow(x),3000,TRUE), ] gc(); system.time(unique(y)) gc(); system.time(Unique(y)) # z = Unique(y) x[matorder(x),] z[matorder(z),]
# `Unique' is faster than `unique' x = matrix(sample(1:(4*8),4*8),ncol=4) y = x[sample(1:nrow(x),3000,TRUE), ] gc(); system.time(unique(y)) gc(); system.time(Unique(y)) # z = Unique(y) x[matorder(x),] z[matorder(z),]