1
@c Copyright (C) 2007-2013 John W. Eaton and David Bateman
3
@c This file is part of Octave.
5
@c Octave is free software; you can redistribute it and/or modify it
6
@c under the terms of the GNU General Public License as published by the
7
@c Free Software Foundation; either version 3 of the License, or (at
8
@c your option) any later version.
10
@c Octave is distributed in the hope that it will be useful, but WITHOUT
11
@c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12
@c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15
@c You should have received a copy of the GNU General Public License
16
@c along with Octave; see the file COPYING. If not, see
17
@c <http://www.gnu.org/licenses/>.
22
Much of the geometry code in Octave is based on the Qhull
23
library@footnote{Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T.,
24
@cite{The Quickhull Algorithm for Convex Hulls}, ACM Trans. on Mathematical
25
Software, 22(4):469--483, Dec 1996, @url{http://www.qhull.org}}.
26
Some of the documentation for Qhull, particularly for the options that
27
can be passed to @code{delaunay}, @code{voronoi} and @code{convhull},
28
etc., is relevant to Octave users.
31
* Delaunay Triangulation::
34
* Interpolation on Scattered Data::
37
@node Delaunay Triangulation
38
@section Delaunay Triangulation
40
The Delaunay triangulation is constructed from a set of
41
circum-circles. These circum-circles are chosen so that there are at
42
least three of the points in the set to triangulation on the
43
circumference of the circum-circle. None of the points in the set of
44
points falls within any of the circum-circles.
46
In general there are only three points on the circumference of any
47
circum-circle. However, in some cases, and in particular for the
48
case of a regular grid, 4 or more points can be on a single
49
circum-circle. In this case the Delaunay triangulation is not unique.
53
The 3- and N-dimensional extension of the Delaunay triangulation are
54
given by @code{delaunay3} and @code{delaunayn} respectively.
55
@code{delaunay3} returns a set of tetrahedra that satisfy the
56
Delaunay circum-circle criteria. Similarly, @code{delaunayn} returns the
57
N-dimensional simplex satisfying the Delaunay circum-circle criteria.
58
The N-dimensional extension of a triangulation is called a tessellation.
64
An example of a Delaunay triangulation of a set of points is
72
X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
73
Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
75
plot (X, Y, "b", x, y, "r*");
81
The result of which can be seen in @ref{fig:delaunay}.
83
@float Figure,fig:delaunay
84
@center @image{delaunay,4in}
85
@caption{Delaunay triangulation of a random set of points}
90
* Plotting the Triangulation::
91
* Identifying Points in Triangulation::
94
@node Plotting the Triangulation
95
@subsection Plotting the Triangulation
97
Octave has the functions @code{triplot}, @code{trimesh}, and @code{trisurf}
98
to plot the Delaunay triangulation of a 2-dimensional set of points.
99
@code{tetramesh} will plot the triangulation of a 3-dimensional set of points.
107
@DOCSTRING(tetramesh)
109
The difference between @code{triplot}, and @code{trimesh} or @code{triplot},
110
is that the former only plots the 2-dimensional triangulation itself, whereas
111
the second two plot the value of a function @code{f (@var{x}, @var{y})}. An
112
example of the use of the @code{triplot} function is
119
tri = delaunay (x, y);
125
which plots the Delaunay triangulation of a set of random points in
128
The output of the above can be seen in @ref{fig:triplot}.
130
@float Figure,fig:triplot
131
@center @image{triplot,4in}
132
@caption{Delaunay triangulation of a random set of points}
136
@node Identifying Points in Triangulation
137
@subsection Identifying Points in Triangulation
139
It is often necessary to identify whether a particular point in the
140
N-dimensional space is within the Delaunay tessellation of a set of
141
points in this N-dimensional space, and if so which N-simplex contains
142
the point and which point in the tessellation is closest to the desired
143
point. The functions @code{tsearch} and @code{dsearch} perform this
144
function in a triangulation, and @code{tsearchn} and @code{dsearchn} in
145
an N-dimensional tessellation.
147
To identify whether a particular point represented by a vector @var{p}
148
falls within one of the simplices of an N-simplex, we can write the
149
Cartesian coordinates of the point in a parametric form with respect to
150
the N-simplex. This parametric form is called the Barycentric
151
Coordinates of the point. If the points defining the N-simplex are given
152
by @code{@var{N} + 1} vectors @var{t}(@var{i},:), then the Barycentric
153
coordinates defining the point @var{p} are given by
156
@var{p} = sum (@var{beta}(1:@var{N}+1) * @var{t}(1:@var{N}+1),:)
160
where there are @code{@var{N} + 1} values @code{@var{beta}(@var{i})}
161
that together as a vector represent the Barycentric coordinates of the
162
point @var{p}. To ensure a unique solution for the values of
163
@code{@var{beta}(@var{i})} an additional criteria of
166
sum (@var{beta}(1:@var{N}+1)) == 1
170
is imposed, and we can therefore write the above as
174
@var{p} - @var{t}(end, :) = @var{beta}(1:end-1) * (@var{t}(1:end-1, :)
175
- ones (@var{N}, 1) * @var{t}(end, :)
180
Solving for @var{beta} we can then write
184
@var{beta}(1:end-1) = (@var{p} - @var{t}(end, :)) / (@var{t}(1:end-1, :)
185
- ones (@var{N}, 1) * @var{t}(end, :))
186
@var{beta}(end) = sum (@var{beta}(1:end-1))
191
which gives the formula for the conversion of the Cartesian coordinates
192
of the point @var{p} to the Barycentric coordinates @var{beta}. An
193
important property of the Barycentric coordinates is that for all points
197
0 <= @var{beta}(@var{i}) <= 1
201
Therefore, the test in @code{tsearch} and @code{tsearchn} essentially
202
only needs to express each point in terms of the Barycentric coordinates
203
of each of the simplices of the N-simplex and test the values of
204
@var{beta}. This is exactly the implementation used in
205
@code{tsearchn}. @code{tsearch} is optimized for 2-dimensions and the
206
Barycentric coordinates are not explicitly formed.
212
An example of the use of @code{tsearch} can be seen with the simple
217
@var{x} = [-1; -1; 1; 1];
218
@var{y} = [-1; 1; -1; 1];
219
@var{tri} = [1, 2, 3; 2, 3, 1];
224
consisting of two triangles defined by @var{tri}. We can then identify
225
which triangle a point falls in like
229
tsearch (@var{x}, @var{y}, @var{tri}, -0.5, -0.5)
231
tsearch (@var{x}, @var{y}, @var{tri}, 0.5, 0.5)
237
and we can confirm that a point doesn't lie within one of the triangles like
241
tsearch (@var{x}, @var{y}, @var{tri}, 2, 2)
246
The @code{dsearch} and @code{dsearchn} find the closest point in a
247
tessellation to the desired point. The desired point does not
248
necessarily have to be in the tessellation, and even if it the returned
249
point of the tessellation does not have to be one of the vertexes of the
250
N-simplex within which the desired point is found.
256
An example of the use of @code{dsearch}, using the above values of
257
@var{x}, @var{y} and @var{tri} is
261
dsearch (@var{x}, @var{y}, @var{tri}, -2, -2)
266
If you wish the points that are outside the tessellation to be flagged,
267
then @code{dsearchn} can be used as
271
dsearchn ([@var{x}, @var{y}], @var{tri}, [-2, -2], NaN)
273
dsearchn ([@var{x}, @var{y}], @var{tri}, [-0.5, -0.5], NaN)
279
where the point outside the tessellation are then flagged with @code{NaN}.
281
@node Voronoi Diagrams
282
@section Voronoi Diagrams
284
A Voronoi diagram or Voronoi tessellation of a set of points @var{s} in
285
an N-dimensional space, is the tessellation of the N-dimensional space
286
such that all points in @code{@var{v}(@var{p})}, a partitions of the
287
tessellation where @var{p} is a member of @var{s}, are closer to @var{p}
288
than any other point in @var{s}. The Voronoi diagram is related to the
289
Delaunay triangulation of a set of points, in that the vertexes of the
290
Voronoi tessellation are the centers of the circum-circles of the
291
simplices of the Delaunay tessellation.
297
An example of the use of @code{voronoi} is
304
tri = delaunay (x, y);
305
[vx, vy] = voronoi (x, y, tri);
306
triplot (tri, x, y, "b");
314
The result of which can be seen in @ref{fig:voronoi}. Note that the
315
circum-circle of one of the triangles has been added to this figure, to
316
make the relationship between the Delaunay tessellation and the Voronoi
319
@float Figure,fig:voronoi
320
@center @image{voronoi,4in}
321
@caption{Delaunay triangulation and Voronoi diagram of a random set of points}
325
Additional information about the size of the facets of a Voronoi
326
diagram, and which points of a set of points is in a polygon can be had
327
with the @code{polyarea} and @code{inpolygon} functions respectively.
331
An example of the use of @code{polyarea} might be
338
[c, f] = voronoin ([x, y]);
339
af = zeros (size (f));
340
for i = 1 : length (f)
341
af(i) = polyarea (c (f @{i, :@}, 1), c (f @{i, :@}, 2));
346
Facets of the Voronoi diagram with a vertex at infinity have infinity
347
area. A simplified version of @code{polyarea} for rectangles is
348
available with @code{rectint}
352
@DOCSTRING(inpolygon)
354
An example of the use of @code{inpolygon} might be
361
vx = cos (pi * [-1 : 0.1: 1]);
362
vy = sin (pi * [-1 : 0.1 : 1]);
363
in = inpolygon (x, y, vx, vy);
364
plot (vx, vy, x(in), y(in), "r+", x(!in), y(!in), "bo");
365
axis ([-2, 2, -2, 2]);
371
The result of which can be seen in @ref{fig:inpolygon}.
373
@float Figure,fig:inpolygon
374
@center @image{inpolygon,4in}
375
@caption{Demonstration of the @code{inpolygon} function to determine the
376
points inside a polygon}
383
The convex hull of a set of points is the minimum convex envelope
384
containing all of the points. Octave has the functions @code{convhull}
385
and @code{convhulln} to calculate the convex hull of 2-dimensional and
386
N-dimensional sets of points.
390
@DOCSTRING(convhulln)
392
An example of the use of @code{convhull} is
399
plot (x(k), y(k), "r-", x, y, "b+");
400
axis ([-3.05, 3.05, -0.05, 1.05]);
406
The output of the above can be seen in @ref{fig:convhull}.
408
@float Figure,fig:convhull
409
@center @image{convhull,4in}
410
@caption{The convex hull of a simple set of points}
414
@node Interpolation on Scattered Data
415
@section Interpolation on Scattered Data
417
An important use of the Delaunay tessellation is that it can be used to
418
interpolate from scattered data to an arbitrary set of points. To do
419
this the N-simplex of the known set of points is calculated with
420
@code{delaunay}, @code{delaunay3} or @code{delaunayn}. Then the
421
simplices in to which the desired points are found are
422
identified. Finally the vertices of the simplices are used to
423
interpolate to the desired points. The functions that perform this
424
interpolation are @code{griddata}, @code{griddata3} and
429
@DOCSTRING(griddata3)
431
@DOCSTRING(griddatan)
433
An example of the use of the @code{griddata} function is
438
x = 2*rand (1000,1) - 1;
439
y = 2*rand (size (x)) - 1;
440
z = sin (2*(x.^2+y.^2));
441
[xx,yy] = meshgrid (linspace (-1,1,32));
442
griddata (x,y,z,xx,yy);
447
that interpolates from a random scattering of points, to a uniform
450
The output of the above can be seen in @ref{fig:griddata}.
452
@float Figure,fig:griddata
453
@center @image{griddata,4in}
454
@caption{Interpolation from a scattered data to a regular grid}
added
removed
doc/interpreter/geometry.txi
