| MATLAB® | |
| On this page… |
|---|
Functions for Analysis of Closest-Point Problems and Geometric Analysis
Function | Description |
|---|---|
Convex hull | |
Delaunay triangulation | |
3-D Delaunay tessellation | |
Nearest point search of Delaunay triangulation | |
True for points inside polygonal region | |
Area of polygon | |
Area of intersection for two or more rectangles | |
Closest triangle search | |
Voronoi diagram |
Note Examples in this section use the MATLAB® seamount data set. Seamounts are underwater mountains. They are valuable sources of information about marine geology. The seamount data set represents the surface, in 1984, of the seamount designated LR148.8W located at 48.2°S, 148.8°W on the Louisville Ridge in the South Pacific. The seamount data set provides longitude (x), latitude (y) and depth-in-feet (z) data for 294 points on the seamount LR148.8W. |
The convhull function returns the indices of the points in a data set that comprise the convex hull for the set. Use the plot function to plot the output of convhull.
This example loads the seamount data and plots the longitudinal (x) and latitudinal (y) data as a scatter plot. It then generates the convex hull and uses plot to plot the convex hull:
load seamount plot(x,y,'.','markersize',10) k = convhull(x,y); hold on, plot(x(k),y(k),'-r'), hold off grid on
Given a set of coplanar data points, Delaunay triangulation is a set of lines connecting each point to its natural neighbors. The delaunay function returns a Delaunay triangulation as a set of triangles having the property that, for each triangle, the unique circle circumscribed about the triangle contains no data points.
You can use triplot to print the resulting triangles in two-dimensional space. You can also add data for a third dimension to the output of delaunay and plot the result as a surface with trisurf, or as a mesh with trimesh.
To try delaunay, load the seamount data set and view the longitude (x) and latitude (y) data as a scatter plot:
load seamount
plot(x,y,'.','markersize',12)
xlabel('Longitude'), ylabel('Latitude')
grid on
Apply Delaunay triangulation and use triplot to overplot the resulting triangles on the scatter plot:
tri = delaunay(x,y); hold on, triplot(tri,x,y), hold off
Add the depth data (z) from seamount, to the Delaunay triangulation, and use trimesh to produce a mesh in three-dimensional space. Similarly, you can use trisurf to produce a surface:
figure
hidden on
trimesh(tri,x,y,z)
grid on
xlabel('Longitude'); ylabel('Latitude'); zlabel('Depth in Feet')
This code uses meshgrid, griddata, and contour to produce a contour plot of the seamount data:
figure
[xi,yi] = meshgrid(210.8:.01:211.8,-48.5:.01:-47.9);
zi = griddata(x,y,z,xi,yi,'cubic');
[c,h] = contour(xi,yi,zi,'b-');
clabel(c,h)
xlabel('Longitude'), ylabel('Latitude')
The arguments for meshgrid encompass the largest and smallest x and y values in the original seamount data. To obtain these values, use min(x), max(x), min(y), and max(y).
You can search through the Delaunay triangulation data with two functions:
dsearch finds the indices of the (x,y) points in a Delaunay triangulation closest to the points you specify. This code searches for the point closest to (211.32, -48.35) in the triangulation of the seamount data.
xi = 211.32; yi = -48.35;
p = dsearch(x,y,tri,xi,yi);
[x(p), y(p)]
ans =
211.3400 -48.3700
tsearch finds the indices into the delaunay output that specify the enclosing triangles of the points you specify. This example uses the index of the enclosing triangle for the point (211.32, -48.35) to obtain the coordinates of the vertices of the triangle:
xi = 211.32; yi = -48.35; t = tsearch(x,y,tri,xi,yi); r = tri(t,:); A = [x(r) y(r)] A = 211.2800 -48.3200 211.3400 -48.3700 211.3000 -48.3000
Voronoi diagrams are a closest-point plotting technique related to Delaunay triangulation.
For each point in a set of coplanar points, you can draw a polygon that encloses all the intermediate points that are closer to that point than to any other point in the set. Such a polygon is called a Voronoi polygon, and the set of all Voronoi polygons for a given point set is called a Voronoi diagram.
The voronoi function can plot the cells of the Voronoi diagram, or return the vertices of the edges of the diagram. This example loads the seamount data, then uses the voronoi function to produce the Voronoi diagram for the longitudinal (x) and latitudinal (y) dimensions. Note that voronoi plots only the bounded cells of the Voronoi diagram:
load seamount
voronoi(x,y)
grid on
xlabel('Longitude'), ylabel('Latitude')
Note See the voronoi function for an example that uses the vertices of the edges to plot a Voronoi diagram. |
Many applications in science, engineering, statistics, and mathematics require structures like convex hulls, Voronoi diagrams, and Delaunay tessellations.
Functions for Multidimensional Geometrical Analysis
Function | Description |
|---|---|
N-dimensional convex hull | |
N-dimensional Delaunay tessellation | |
N-dimensional nearest point search | |
N-dimensional data gridding and hypersurface fitting | |
N-dimensional closest simplex search | |
N-dimensional Voronoi diagrams |
This section demonstrates these geometric analysis techniques:
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 =
4 2 1
3 4 1
7 3 1
5 7 1
7 4 3
4 7 8
2 6 1
6 5 1
4 6 2
6 4 8
6 7 5
7 6 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, 4 2 1, means that the first triangle has X(4,:), X(2,:), and X(1,:) 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' adds transparency
hold off
view(3), axis equal, axis off
camorbit(90,-5); % To view it from another angle
title('Convex hull of a cube')
A Delaunay tessellation is a set of simplices with the property that, for each simplex, the unique sphere circumscribed about the simplex contains no data points. 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 Computing a Convex Hull, 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 =
4 3 9 1
4 9 2 1
7 9 3 1
7 5 9 1
7 9 4 3
7 8 4 9
6 2 9 1
6 9 5 1
6 4 9 2
6 4 8 9
6 9 7 5
6 8 7 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 Computing a Delaunay Tessellation, 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
rand('twister',5489); % Initialize the generator
X = X+0.01*rand(size(X)); % Make the cube less regular
[V,C] = voronoin(X)
V =
Inf Inf Inf
0.0029 -1.4858 0.0079
0.1702 -0.0719 193.1848
0.0018 0.0089 1.5064
0.0067 0.0040 1.5064
0.0060 2.9397 0.0073
0.0033 1.5095 0.0119
0.0117 1.5095 0.0035
-1.4873 0.0055 0.0107
-1.6607 0.0051 0.0100
-1.4873 0.0050 0.0101
0.0054 -1.4858 0.0054
4.6864 -0.0017 0.0080
1.5105 0.0107 0.0022
1.5105 0.0025 0.0104
0.0108 0.0120 -1.4850
0.0032 0.0070 -3.1326
0.0043 0.0056 -1.4849
C =
[1x7 double]
[1x10 double]
[1x8 double]
[1x7 double]
[1x8 double]
[1x7 double]
[1x7 double]
[1x10 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 ninth cell of the Voronoi diagram is
C{9} = 2 4 5 7 8 9 11 12 14 15 16 18
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')
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:
n = 5000; X = 2*rand(n,3)-1; v = sum(X.^2,2);
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')
Note A smaller delta produces a smoother sphere, but increases the compute time. |
| Gridded Data | Function Functions | |
| © 1984-2008- The MathWorks, Inc. - Site Help - Patents - Trademarks - Privacy Policy - Preventing Piracy - RSS |