|
7018
|
1 @c Copyright (C) 2007 John W. Eaton and David Bateman |
|
|
2 @c |
|
|
3 @c This file is part of Octave. |
|
|
4 @c |
|
|
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. |
|
|
9 @c |
|
|
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 |
|
|
13 @c for more details. |
|
|
14 @c |
|
|
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/>. |
|
6558
|
18 |
|
|
19 @node Geometry |
|
|
20 @chapter Geometry |
|
|
21 |
|
6832
|
22 Much of geometry code in Octave is based on the QHull @footnote{Barber, |
|
|
23 C.B., Dobkin, D.P., and Huhdanpaa, H.T., "The Quickhull algorithm for |
|
|
24 convex hulls," ACM Trans. on Mathematical Software, 22(4):469-483, Dec |
|
|
25 1996, @url{http://www.qhull.org}}. Some of the documentation for Qhull, |
|
|
26 particularly for the options that can be passed to @code{delaunay}, |
|
|
27 @code{voronoi} and @code{convhull}, etc, is relevant to Octave users. |
|
|
28 |
|
6823
|
29 @menu |
|
|
30 * Delaunay Triangulation:: |
|
|
31 * Voronoi Diagrams:: |
|
|
32 * Convex Hull:: |
|
|
33 * Interpolation on Scattered Data:: |
|
|
34 @end menu |
|
|
35 |
|
|
36 @node Delaunay Triangulation |
|
|
37 @section Delaunay Triangulation |
|
|
38 |
|
6832
|
39 The Delaunay triangulation is constructed from a set of |
|
|
40 circum-circles. These circum-circles are chosen so that there are at |
|
|
41 least three of the points in the set to triangulation on the |
|
|
42 circumference of the circum-circle. None of the points in the set of |
|
|
43 points falls within any of the circum-circles. |
|
|
44 |
|
|
45 In general there are only three points on the circumference of any |
|
|
46 circum-circle. However, in the some cases, and in particular for the |
|
|
47 case of a regular grid, 4 or more points can be on a single |
|
|
48 circum-circle. In this case the Delaunay triangulation is not unique. |
|
|
49 |
|
6823
|
50 @DOCSTRING(delaunay) |
|
|
51 |
|
|
52 The 3- and N-dimensional extension of the Delaunay triangulation are |
|
|
53 given by @code{delaunay3} and @code{delaunayn} respectively. |
|
|
54 @code{delaunay3} returns a set of tetrahedra that satisfy the |
|
|
55 Delaunay circum-circle criteria. Similarly, @code{delaunayn} returns the |
|
|
56 N-dimensional simplex satisfying the Delaunay circum-circle criteria. |
|
7007
|
57 The N-dimensional extension of a triangulation is called a tessellation. |
|
6823
|
58 |
|
|
59 @DOCSTRING(delaunay3) |
|
|
60 |
|
|
61 @DOCSTRING(delaunayn) |
|
|
62 |
|
6832
|
63 An example of a Delaunay triangulation of a set of points is |
|
|
64 |
|
|
65 @example |
|
|
66 @group |
|
|
67 rand ("state", 2); |
|
|
68 x = rand (10, 1); |
|
|
69 y = rand (10, 1); |
|
|
70 T = delaunay (x, y); |
|
|
71 X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; |
|
|
72 Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; |
|
|
73 axis ([0, 1, 0, 1]); |
|
|
74 plot(X, Y, "b", x, y, "r*"); |
|
|
75 @end group |
|
|
76 @end example |
|
|
77 |
|
6855
|
78 @ifset HAVE_QHULL |
|
6832
|
79 @ifnotinfo |
|
|
80 @noindent |
|
|
81 The result of which can be seen in @ref{fig:delaunay}. |
|
|
82 |
|
|
83 @float Figure,fig:delaunay |
|
|
84 @image{delaunay,8cm} |
|
|
85 @caption{Delaunay triangulation of a random set of points} |
|
|
86 @end float |
|
|
87 @end ifnotinfo |
|
6855
|
88 @end ifset |
|
6832
|
89 |
|
6823
|
90 @menu |
|
6832
|
91 * Plotting the Triangulation:: |
|
6823
|
92 * Identifying points in Triangulation:: |
|
|
93 @end menu |
|
|
94 |
|
6832
|
95 @node Plotting the Triangulation |
|
|
96 @subsection Plotting the Triangulation |
|
|
97 |
|
|
98 Octave has the functions @code{triplot} and @code{trimesh} to plot the |
|
|
99 Delaunay triangulation of a 2-dimensional set of points. |
|
|
100 |
|
|
101 @DOCSTRING(triplot) |
|
|
102 |
|
|
103 @DOCSTRING(trimesh) |
|
|
104 |
|
|
105 The difference between @code{triplot} and @code{trimesh} is that the |
|
|
106 former only plots the 2-dimensional triangulation itself, whereas the |
|
|
107 second plots the value of some function @code{f (@var{x}, @var{y})}. |
|
|
108 An example of the use of the @code{triplot} function is |
|
|
109 |
|
|
110 @example |
|
|
111 @group |
|
|
112 rand ("state", 2) |
|
|
113 x = rand (20, 1); |
|
|
114 y = rand (20, 1); |
|
|
115 tri = delaunay (x, y); |
|
|
116 triplot (tri, x, y); |
|
|
117 @end group |
|
|
118 @end example |
|
|
119 |
|
|
120 that plot the Delaunay triangulation of a set of random points in |
|
|
121 2-dimensions. |
|
|
122 @ifnotinfo |
|
|
123 The output of the above can be seen in @ref{fig:triplot}. |
|
|
124 |
|
|
125 @float Figure,fig:triplot |
|
|
126 @image{triplot,8cm} |
|
|
127 @caption{Delaunay triangulation of a random set of points} |
|
|
128 @end float |
|
|
129 @end ifnotinfo |
|
|
130 |
|
6823
|
131 @node Identifying points in Triangulation |
|
|
132 @subsection Identifying points in Triangulation |
|
|
133 |
|
|
134 It is often necessary to identify whether a particular point in the |
|
7007
|
135 N-dimensional space is within the Delaunay tessellation of a set of |
|
6823
|
136 points in this N-dimensional space, and if so which N-Simplex contains |
|
7007
|
137 the point and which point in the tessellation is closest to the desired |
|
6823
|
138 point. The functions @code{tsearch} and @code{dsearch} perform this |
|
|
139 function in a triangulation, and @code{tsearchn} and @code{dsearchn} in |
|
7007
|
140 an N-dimensional tessellation. |
|
6823
|
141 |
|
|
142 To identify whether a particular point represented by a vector @var{p} |
|
|
143 falls within one of the simplices of an N-Simplex, we can write the |
|
|
144 Cartesian coordinates of the point in a parametric form with respect to |
|
|
145 the N-Simplex. This parametric form is called the Barycentric |
|
|
146 Coordinates of the point. If the points defining the N-Simplex are given |
|
|
147 by @code{@var{N} + 1} vectors @var{t}(@var{i},:), then the Barycentric |
|
|
148 coordinates defining the point @var{p} is given by |
|
|
149 |
|
|
150 @example |
|
|
151 @var{p} = sum (@var{beta}(1:@var{N}+1) * @var{t}(1:@var{N}+1),:) |
|
|
152 @end example |
|
|
153 |
|
|
154 @noindent |
|
|
155 where there are @code{@var{N} + 1} values @code{@var{beta}(@var{i})} |
|
|
156 that together as a vector represent the Barycentric coordinates of the |
|
|
157 point @var{p}. To ensure a unique solution for the values of |
|
|
158 @code{@var{beta}(@var{i})} an additional criteria of |
|
|
159 |
|
|
160 @example |
|
|
161 sum (@var{beta}(1:@var{N}+1)) == 1 |
|
|
162 @end example |
|
|
163 |
|
|
164 @noindent |
|
|
165 is imposed, and we can therefore write the above as |
|
|
166 |
|
|
167 @example |
|
|
168 @var{p} - @var{t}(end, :) = @var{beta}(1:end-1) * (@var{t}(1:end-1, :) |
|
|
169 - ones(@var{N}, 1) * @var{t}(end, :) |
|
|
170 @end example |
|
|
171 |
|
|
172 @noindent |
|
|
173 Solving for @var{beta} we can then write |
|
|
174 |
|
|
175 @example |
|
|
176 @var{beta}(1:end-1) = (@var{p} - @var{t}(end, :)) / (@var{t}(1:end-1, :) |
|
|
177 - ones(@var{N}, 1) * @var{t}(end, :)) |
|
|
178 @var{beta}(end) = sum(@var{beta}(1:end-1)) |
|
|
179 @end example |
|
|
180 |
|
|
181 @noindent |
|
|
182 which gives the formula for the conversion of the Cartesian coordinates |
|
|
183 of the point @var{p} to the Barycentric coordinates @var{beta}. An |
|
|
184 important property of the Barycentric coordinates is that for all points |
|
|
185 in the N-Simplex |
|
|
186 |
|
|
187 @example |
|
|
188 0 <= @var{beta}(@var{i}) <= 1 |
|
|
189 @end example |
|
|
190 |
|
|
191 @noindent |
|
|
192 Therefore, the test in @code{tsearch} and @code{tsearchn} essentially |
|
|
193 only needs to express each point in terms of the Barycentric coordinates |
|
|
194 of each of the simplices of the N-Simplex and test the values of |
|
|
195 @var{beta}. This is exactly the implementation used in |
|
|
196 @code{tsearchn}. @code{tsearch} is optimized for 2-dimensions and the |
|
|
197 Barycentric coordinates are not explicitly formed. |
|
|
198 |
|
|
199 @DOCSTRING(tsearch) |
|
|
200 |
|
|
201 @DOCSTRING(tsearchn) |
|
|
202 |
|
|
203 An example of the use of @code{tsearch} can be seen with the simple |
|
|
204 triangulation |
|
|
205 |
|
|
206 @example |
|
|
207 @group |
|
|
208 @var{x} = [-1; -1; 1; 1]; |
|
|
209 @var{y} = [-1; 1; -1; 1]; |
|
|
210 @var{tri} = [1, 2, 3; 2, 3, 1]; |
|
|
211 @end group |
|
|
212 @end example |
|
|
213 |
|
|
214 @noindent |
|
|
215 consisting of two triangles defined by @var{tri}. We can then identify |
|
|
216 which triangle a point falls in like |
|
|
217 |
|
|
218 @example |
|
|
219 @group |
|
|
220 tsearch (@var{x}, @var{y}, @var{tri}, -0.5, -0.5) |
|
|
221 @result{} 1 |
|
|
222 tsearch (@var{x}, @var{y}, @var{tri}, 0.5, 0.5) |
|
|
223 @result{} 2 |
|
|
224 @end group |
|
|
225 @end example |
|
|
226 |
|
|
227 @noindent |
|
|
228 and we can confirm that a point doesn't lie within one of the triangles like |
|
|
229 |
|
|
230 @example |
|
|
231 @group |
|
|
232 tsearch (@var{x}, @var{y}, @var{tri}, 2, 2) |
|
|
233 @result{} NaN |
|
|
234 @end group |
|
|
235 @end example |
|
|
236 |
|
|
237 The @code{dsearch} and @code{dsearchn} find the closest point in a |
|
|
238 tessellation to the desired point. The desired point does not |
|
|
239 necessarily have to be in the tessellation, and even if it the returned |
|
6832
|
240 point of the tessellation does not have to be one of the vertexes of the |
|
6823
|
241 N-simplex within which the desired point is found. |
|
|
242 |
|
|
243 @DOCSTRING(dsearch) |
|
|
244 |
|
|
245 @DOCSTRING(dsearchn) |
|
|
246 |
|
|
247 An example of the use of @code{dsearch}, using the above values of |
|
|
248 @var{x}, @var{y} and @var{tri} is |
|
|
249 |
|
|
250 @example |
|
|
251 @group |
|
|
252 dsearch (@var{x}, @var{y}, @var{tri}, -2, -2) |
|
|
253 @result{} 1 |
|
|
254 @end group |
|
|
255 @end example |
|
|
256 |
|
|
257 If you wish the points that are outside the tessellation to be flagged, |
|
|
258 then @code{dsearchn} can be used as |
|
|
259 |
|
|
260 @example |
|
|
261 @group |
|
|
262 dsearchn ([@var{x}, @var{y}], @var{tri}, [-2, -2], NaN) |
|
|
263 @result{} NaN |
|
|
264 dsearchn ([@var{x}, @var{y}], @var{tri}, [-0.5, -0.5], NaN) |
|
|
265 @result{} 1 |
|
|
266 @end group |
|
|
267 @end example |
|
|
268 |
|
|
269 @noindent |
|
|
270 where the point outside the tessellation are then flagged with @code{NaN}. |
|
|
271 |
|
|
272 @node Voronoi Diagrams |
|
|
273 @section Voronoi Diagrams |
|
|
274 |
|
|
275 A Voronoi diagram or Voronoi tessellation of a set of points @var{s} in |
|
|
276 an N-dimensional space, is the tessellation of the N-dimensional space |
|
|
277 such that all points in @code{@var{v}(@var{p})}, a partitions of the |
|
|
278 tessellation where @var{p} is a member of @var{s}, are closer to @var{p} |
|
|
279 than any other point in @var{s}. The Voronoi diagram is related to the |
|
6832
|
280 Delaunay triangulation of a set of points, in that the vertexes of the |
|
|
281 Voronoi tessellation are the center's of the circum-circles of the |
|
|
282 simplicies of the Delaunay tessellation. |
|
6823
|
283 |
|
|
284 @DOCSTRING(voronoi) |
|
|
285 |
|
|
286 @DOCSTRING(voronoin) |
|
|
287 |
|
6832
|
288 An example of the use of @code{voronoi} is |
|
|
289 |
|
|
290 @example |
|
|
291 @group |
|
|
292 rand("state",9); |
|
|
293 x = rand(10,1); |
|
|
294 y = rand(10,1); |
|
|
295 tri = delaunay (x, y); |
|
|
296 [vx, vy] = voronoi (x, y, tri); |
|
|
297 triplot (tri, x, y, "b"); |
|
|
298 hold on; |
|
|
299 plot (vx, vy, "r"); |
|
|
300 @end group |
|
|
301 @end example |
|
|
302 |
|
6855
|
303 @ifset HAVE_QHULL |
|
6832
|
304 @ifnotinfo |
|
|
305 @noindent |
|
|
306 The result of which can be seen in @ref{fig:voronoi}. Note that the |
|
|
307 circum-circle of one of the triangles has been added to this figure, to |
|
|
308 make the relationship between the Delaunay tessellation and the Voronoi |
|
|
309 diagram clearer. |
|
|
310 |
|
|
311 @float Figure,fig:voronoi |
|
|
312 @image{voronoi,8cm} |
|
|
313 @caption{Delaunay triangulation and Voronoi diagram of a random set of points} |
|
|
314 @end float |
|
|
315 @end ifnotinfo |
|
6855
|
316 @end ifset |
|
6832
|
317 |
|
6847
|
318 Additional information about the size of the facets of a Voronoi |
|
|
319 diagram, and which points of a set of points is in a polygon can be had |
|
|
320 with the @code{polyarea} and @code{inpolygon} functions respectively. |
|
6832
|
321 |
|
|
322 @DOCSTRING(polyarea) |
|
|
323 |
|
|
324 An example of the use of @code{polyarea} might be |
|
|
325 |
|
|
326 @example |
|
|
327 @group |
|
|
328 rand ("state", 2); |
|
|
329 x = rand (10, 1); |
|
|
330 y = rand (10, 1); |
|
|
331 [c, f] = voronoin ([x, y]); |
|
|
332 af = zeros (size(f)); |
|
|
333 for i = 1 : length (f) |
|
|
334 af(i) = polyarea (c (f @{i, :@}, 1), c (f @{i, :@}, 2)); |
|
|
335 endfor |
|
|
336 @end group |
|
|
337 @end example |
|
|
338 |
|
|
339 Facets of the Voronoi diagram with a vertex at infinity have infinity area. |
|
|
340 |
|
6847
|
341 @DOCSTRING(inpolygon) |
|
|
342 |
|
|
343 An example of the use of @code{inpolygon} might be |
|
|
344 |
|
|
345 @example |
|
|
346 @group |
|
|
347 randn ("state", 2); |
|
|
348 x = randn (100, 1); |
|
|
349 y = randn (100, 1); |
|
|
350 vx = cos (pi * [-1 : 0.1: 1]); |
|
|
351 vy = sin (pi * [-1 : 0.1 : 1]); |
|
|
352 in = inpolygon (x, y, vx, vy); |
|
|
353 plot(vx, vy, x(in), y(in), "r+", x(!in), y(!in), "bo"); |
|
|
354 axis ([-2, 2, -2, 2]); |
|
|
355 @end group |
|
|
356 @end example |
|
|
357 |
|
|
358 @ifnotinfo |
|
|
359 @noindent |
|
|
360 The result of which can be seen in @ref{fig:inpolygon}. |
|
|
361 |
|
|
362 @float Figure,fig:inpolygon |
|
|
363 @image{inpolygon,8cm} |
|
|
364 @caption{Demonstration of the @code{inpolygon} function to determine the |
|
|
365 points inside a polygon} |
|
|
366 @end float |
|
|
367 @end ifnotinfo |
|
|
368 |
|
6823
|
369 @node Convex Hull |
|
|
370 @section Convex Hull |
|
|
371 |
|
7001
|
372 The convex hull of a set of points is the minimum convex envelope |
|
6832
|
373 containing all of the points. Octave has the functions @code{convhull} |
|
7007
|
374 and @code{convhulln} to calculate the convex hull of 2-dimensional and |
|
6832
|
375 N-dimensional sets of points. |
|
|
376 |
|
6823
|
377 @DOCSTRING(convhull) |
|
|
378 |
|
|
379 @DOCSTRING(convhulln) |
|
|
380 |
|
6832
|
381 An example of the use of @code{convhull} is |
|
6823
|
382 |
|
6832
|
383 @example |
|
|
384 @group |
|
|
385 x = -3:0.05:3; |
|
|
386 y = abs (sin (x)); |
|
|
387 k = convhull (x, y); |
|
|
388 plot (x(k), y(k), "r-", x, y, "b+"); |
|
|
389 axis ([-3.05, 3.05, -0.05, 1.05]); |
|
|
390 @end group |
|
|
391 @end example |
|
6823
|
392 |
|
6855
|
393 @ifset HAVE_QHULL |
|
6832
|
394 @ifnotinfo |
|
|
395 @noindent |
|
|
396 The output of the above can be seen in @ref{fig:convhull}. |
|
6823
|
397 |
|
6832
|
398 @float Figure,fig:convhull |
|
|
399 @image{convhull,8cm} |
|
|
400 @caption{The convex hull of a simple set of points} |
|
|
401 @end float |
|
|
402 @end ifnotinfo |
|
6855
|
403 @end ifset |
|
6823
|
404 |
|
|
405 @node Interpolation on Scattered Data |
|
|
406 @section Interpolation on Scattered Data |
|
|
407 |
|
6832
|
408 An important use of the Delaunay tessellation is that it can be used to |
|
|
409 interpolate from scattered data to an arbitrary set of points. To do |
|
|
410 this the N-simplex of the known set of points is calculated with |
|
|
411 @code{delaunay}, @code{delaunay3} or @code{delaunayn}. Then the |
|
|
412 simplicies in to which the desired points are found are |
|
|
413 identified. Finally the vertices of the simplicies are used to |
|
|
414 interpolate to the desired points. The functions that perform this |
|
|
415 interpolation are @code{griddata}, @code{griddata3} and |
|
|
416 @code{griddatan}. |
|
|
417 |
|
6823
|
418 @DOCSTRING(griddata) |
|
|
419 |
|
|
420 @DOCSTRING(griddata3) |
|
|
421 |
|
|
422 @DOCSTRING(griddatan) |
|
6832
|
423 |
|
|
424 An example of the use of the @code{griddata} function is |
|
|
425 |
|
|
426 @example |
|
|
427 @group |
|
|
428 rand("state",1); |
|
|
429 x=2*rand(1000,1)-1; |
|
|
430 y=2*rand(size(x))-1; |
|
|
431 z=sin(2*(x.^2+y.^2)); |
|
|
432 [xx,yy]=meshgrid(linspace(-1,1,32)); |
|
|
433 griddata(x,y,z,xx,yy); |
|
|
434 @end group |
|
|
435 @end example |
|
|
436 |
|
6855
|
437 @ifset HAVE_QHULL |
|
6832
|
438 @noindent |
|
|
439 that interpolates from a random scattering of points, to a uniform |
|
|
440 grid. |
|
|
441 @ifnotinfo |
|
|
442 The output of the above can be seen in @ref{fig:griddata}. |
|
|
443 |
|
|
444 @float Figure,fig:griddata |
|
|
445 @image{griddata,8cm} |
|
|
446 @caption{Interpolation from a scattered data to a regular grid} |
|
|
447 @end float |
|
|
448 @end ifnotinfo |
|
6855
|
449 @end ifset |
