| 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 |
|
n-D convex hull. |
|
n-D Delaunay tessellation. |
|
n-D nearest point search. |
|
n-D data gridding and hypersurface fitting. |
|
n-D closest simplex search. |
|
n-D Voronoi diagrams. |
This section demonstrates these geometric analysis techniques:
Convex Hulls
The convex hull of a data set in n-dimensional 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 8-by-3 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 three-dimensional, 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 three-dimensional 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 three-dimensional 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 two-dimensional space, a simplex is a triangle. In three-dimensional 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 n-dimensional 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 three-dimensional 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 three-dimensional 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 n-dimensional space, a Voronoi diagram is the partition of n-dimensional 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 m-by-n matrix of m points in n-space. 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 13-by-3 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 9-by-1 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 N-Dimensional 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-by-3 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 three-dimensional 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 | |
