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 `convhulln` n-D convex hull. `delaunayn` n-D Delaunay tessellation. `dsearchn` n-D nearest point search. `griddatan` n-D data gridding and hypersurface fitting. `tsearchn` n-D closest simplex search. `voronoin` 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.

Voronoi Diagrams

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

• ```C{9} = 2 3 4 5 6 7 8 9 10 11 12 13
```

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.

• ```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.

 Triangulation and Interpolation of Scattered Data Selected Bibliography