Mathematics 
Tessellation and Interpolation of Scattered Data in Higher Dimensions
Many applications in science, engineering, statistics, and mathematics require structures like convex hulls, Voronoi diagrams, and Delaunay tessellations. Using Qhull [1], MATLAB functions enable you to geometrically analyze data sets in any dimension.
Function 
Description 

nD convex hull. 

nD Delaunay tessellation. 

nD nearest point search. 

nD data gridding and hypersurface fitting. 

nD closest simplex search. 

nD Voronoi diagrams. 
This section demonstrates these geometric analysis techniques:
Convex Hulls
The convex hull of a data set in ndimensional space is defined as the smallest convex region that contains the data set.
Computing a Convex Hull. The convhulln
function returns the indices of the points in a data set that comprise the facets of the convex hull for the set. For example, suppose X
is an 8by3 matrix that consists of the 8 vertices of a cube. The convex hull of X
then consists of 12 facets.
d = [1 1]; [x,y,z] = meshgrid(d,d,d); X = [x(:),y(:),z(:)]; % 8 corner points of a cube C = convhulln(X) C = 3 1 5 1 2 5 2 1 3 7 3 5 8 7 5 7 8 3 6 8 5 2 6 5 6 2 8 8 4 3 4 2 3 2 4 8
Because the data is threedimensional, the facets that make up the convex hull are triangles. The 12 rows of C
represent 12 triangles. The elements of C are indices of points in X
. For example, the first row, 3 1 5, means that the first triangle has X(3,:)
, X(1,:)
, and X(5,:)
as its vertices.
For threedimensional convex hulls, you can use trisurf
to plot the output. However, using patch
to plot the output gives you more control over the color of the facets. Note that you cannot plot convhulln
output for n > 3
.
This code plots the convex hull by drawing the triangles as threedimensional patches.
figure, hold on d = [1 2 3 1 % Index into C column. for i = 1:size(C,1) % Draw each triangle. j= C(i,d); % Get the ith C to make a patch. h(i)=patch(X(j,1),X(j,2),X(j,3),i,'FaceAlpha',0.9); end % 'FaceAlpha' is used to make it transparent. hold off view(3), axis equal, axis off camorbit(90,5); % To view it from another angle title('Convex hull of a cube')
Delaunay Tessellations
A Delaunay tessellation is a set of simplices such that no data points are contained in any simplex's circumsphere. In twodimensional space, a simplex is a triangle. In threedimensional space, a simplex is a tetrahedron.
Computing a Delaunay Tessellation. The delaunayn
function returns the indices of the points in a data set that comprise the simplices of an ndimensional Delaunay tessellation of the data set.
This example uses the same X
as in the convex hull example, i.e. the 8 corner points of a cube, with the addition of a center point.
d = [1 1]; [x,y,z] = meshgrid(d,d,d); X = [x(:),y(:),z(:)]; % 8 corner points of a cube X(9,:) = [0 0 0]; % Add center to the vertex list. T = delaunayn(X) % Generate Delaunay tessellation. T = 9 1 5 6 3 9 1 5 2 9 1 6 2 3 9 4 2 3 9 1 7 9 5 6 7 3 9 5 8 7 9 6 8 2 9 6 8 2 9 4 8 3 9 4 8 7 3 9
The 12 rows of T
represent the 12 simplices, in this case irregular tetrahedrons, that partition the cube. Each row represents one tetrahedron, and the row elements are indices of points in X
.
For threedimensional tessellations, you can use tetramesh
to plot the output. However, using patch
to plot the output gives you more control over the color of the facets. Note that you cannot plot delaunayn
output for n > 3
.
This code plots the tessellation T
by drawing the tetrahedrons using threedimensional patches.
figure, hold on d = [1 1 1 2; 2 2 3 3; 3 4 4 4]; % Index into T for i = 1:size(T,1) % Draw each tetrahedron. y = T(i,d); % Get the ith T to make a patch. x1 = reshape(X(y,1),3,4); x2 = reshape(X(y,2),3,4); x3 = reshape(X(y,3),3,4); h(i)=patch(x1,x2,x3,(1:4)*i,'FaceAlpha',0.9); end hold off view(3), axis equal axis off camorbit(65,120) % To view it from another angle title('Delaunay tessellation of a cube with a center point')
You can use cameramenu
to rotate the figure in any direction.
Given m data points in ndimensional space, a Voronoi diagram is the partition of ndimensional space into m polyhedral regions, one region for each data point. Such a region is called a Voronoi cell. A Voronoi cell satisfies the condition that it contains all points that are closer to its data point than any other data point in the set.
Computing a Voronoi Diagram. The voronoin
function returns two outputs:
V
is an mbyn matrix of m points in nspace. Each row of V
represents a Voronoi vertex.
C
is a cell array of vectors. Each vector in the cell array C
represents a Voronoi cell. The vector contains indices of the points in V
that are the vertices of the Voronoi cell. Each Voronoi cell may have a different number of points.
Because a Voronoi cell can be unbounded, the first row of V
is a point at infinity. Then any unbounded Voronoi cell in C
includes the point at infinity, i.e., the first point in V
.
This example uses the same X
as in the Delaunay example, i.e., the 8 corner points of a cube and its center. Random noise is added to make the cube less regular. The resulting Voronoi diagram has 9 Voronoi cells.
d = [1 1]; [x,y,z] = meshgrid(d,d,d); X = [x(:),y(:),z(:)]; % 8 corner points of a cube X(9,:) = [0 0 0]; % Add center to the vertex list. X = X+0.01*rand(size(X)); % Make the cube less regular. [V,C] = voronoin(X); V = Inf Inf Inf 0.0055 1.5054 0.0004 0.0037 0.0101 1.4990 0.0052 0.0087 1.4990 0.0030 1.5054 0.0030 0.0072 0.0072 1.4971 1.7912 0.0000 0.0044 1.4886 0.0011 0.0036 1.4886 0.0002 0.0045 0.0101 0.0044 1.4971 1.5115 0.0074 0.0033 1.5115 0.0081 0.0040 0.0104 1.4846 0.0007 0.0026 1.4846 0.0071 C = [1x8 double] [1x6 double] [1x4 double] [1x6 double] [1x6 double] [1x6 double] [1x6 double] [1x6 double] [1x12 double]
In this example, V
is a 13by3 matrix, the 13 rows are the coordinates of the 13 Voronoi vertices. The first row of V
is a point at infinity. C
is a 9by1 cell array, where each cell in the array contains an index vector into V
corresponding to one of the 9 Voronoi cells. For example, the 9th cell of the Voronoi diagram is
If any index in a cell of the cell array is 1
, then the corresponding Voronoi cell contains the first point in V
, a point at infinity. This means the Voronoi cell is unbounded.
To view a bounded Voronoi cell, i.e., one that does not contain a point at infinity, use the convhulln
function to compute the vertices of the facets that make up the Voronoi cell. Then use patch
and other plot functions to generate the figure. For example, this code plots the Voronoi cell defined by the 9th cell in C
.
X = V(C{9},:); % View 9th Voronoi cell. K = convhulln(X); figure hold on d = [1 2 3 1]; % Index into K for i = 1:size(K,1) j = K(i,d); h(i)=patch(X(j,1),X(j,2),X(j,3),i,'FaceAlpha',0.9); end hold off view(3) axis equal title('One cell of a Voronoi diagram')
Interpolating NDimensional Data
Use the griddatan
function to interpolate multidimensional data, particularly scattered data. griddatan
uses the delaunayn
function to tessellate the data, and then interpolates based on the tessellation.
Suppose you want to visualize a function that you have evaluated at a set of n
scattered points. In this example, X
is an n
by3 matrix of points, each row containing the (x
,y
,z
) coordinates for one of the points. The vector v
contains the n
function values at these points. The function for this example is the squared distance from the origin, v = x.^2 + y.^2 + z.^2
.
Start by generating n
= 5000 points at random in threedimensional space, and computing the value of a function on those points.
The next step is to use interpolation to compute function values over a grid. Use meshgrid
to create the grid, and griddatan
to do the interpolation.
delta = 0.05; d = 1:delta:1; [x0,y0,z0] = meshgrid(d,d,d); X0 = [x0(:), y0(:), z0(:)]; v0 = griddatan(X,v,X0); v0 = reshape(v0, size(x0));
Then use isosurface and related functions to visualize the surface that consists of the (x
,y
,z
) values for which the function takes a constant value. You could pick any value, but the example uses the value 0.6. Since the function is the squared distance from the origin, the surface at a constant value is a sphere.
p = patch(isosurface(x0,y0,z0,v0,0.6)); isonormals(x0,y0,z0,v0,p); set(p,'FaceColor','red','EdgeColor','none'); view(3); camlight; lighting phong axis equal title('Interpolated sphere from scattered data')
Triangulation and Interpolation of Scattered Data  Selected Bibliography 