|
6823
|
1 ## Copyright (C) 2000 Kai Habel |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
8 ## any later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
|
19 |
|
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {} voronoi (@var{x}, @var{y}) |
|
|
22 ## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle") |
|
|
23 ## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle", @var{options}) |
|
|
24 ## @deftypefnx {Function File} {[@var{vx}, @var{vy}] =} voronoi (@dots{}) |
|
|
25 ## plots voronoi diagram of points @code{(@var{x}, @var{y})}. |
|
|
26 ## The voronoi facets with points at infinity are not drawn. |
|
|
27 ## [@var{vx}, @var{vy}] = voronoi(...) returns the vertices instead plotting the |
|
|
28 ## diagram. plot (@var{vx}, @var{vy}) shows the voronoi diagram. |
|
|
29 ## |
|
|
30 ## A fourth optional argument, which must be a string, contains extra options |
|
|
31 ## passed to the underlying qhull command. See the documentation for the |
|
|
32 ## Qhull library for details. |
|
|
33 ## |
|
|
34 ## @example |
|
|
35 ## @group |
|
6826
|
36 ## x = rand (10, 1); |
|
|
37 ## y = rand (size (x)); |
|
|
38 ## h = convhull (x, y); |
|
|
39 ## [vx, vy] = voronoi (x, y); |
|
|
40 ## plot (vx, vy, "-b", x, y, "o", x(h), y(h), "-g") |
|
|
41 ## legend ("", "points", "hull"); |
|
6823
|
42 ## @end group |
|
|
43 ## @end example |
|
|
44 ## |
|
|
45 ## @seealso{voronoin, delaunay, convhull} |
|
|
46 ## @end deftypefn |
|
|
47 |
|
|
48 ## Author: Kai Habel <kai.habel@gmx.de> |
|
|
49 ## First Release: 20/08/2000 |
|
|
50 |
|
|
51 ## 2002-01-04 Paul Kienzle <pkienzle@users.sf.net> |
|
|
52 ## * limit the default graph to the input points rather than the whole diagram |
|
|
53 ## * provide example |
|
|
54 ## * use unique(x,"rows") rather than __unique_rows__ |
|
|
55 |
|
|
56 ## 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net> |
|
|
57 ## Added optional fourth argument to pass options to the underlying |
|
|
58 ## qhull command |
|
|
59 |
|
|
60 function [vvx, vvy] = voronoi (varargin) |
|
|
61 |
|
|
62 if (nargin < 1) |
|
|
63 print_usage (); |
|
|
64 endif |
|
|
65 |
|
|
66 narg = 1; |
|
|
67 if (isscalar (varargin{1}) && ishandle (varargin{1})) |
|
|
68 handl = varargin{1}; |
|
|
69 narg++; |
|
6826
|
70 if (! strcmp (get (handl, "type"), "axes")) |
|
6823
|
71 error ("voronoi: expecting first argument to be an axes object"); |
|
|
72 endif |
|
|
73 else |
|
|
74 if (nargout < 2) |
|
6826
|
75 handl = gca (); |
|
6823
|
76 endif |
|
|
77 endif |
|
|
78 |
|
|
79 if (nargin < 1 + narg || nargin > 3 + narg) |
|
|
80 print_usage (); |
|
|
81 endif |
|
|
82 |
|
|
83 x = varargin{narg++}; |
|
|
84 y = varargin{narg++}; |
|
|
85 |
|
|
86 opts = {}; |
|
|
87 if (narg <= nargin) |
|
|
88 if (iscell (varargin{narg})) |
|
|
89 opts = varargin(narg++); |
|
|
90 elseif (ismatrix (varargin{narg})) |
|
|
91 ## Accept but ignore the triangulation |
|
|
92 narg++; |
|
|
93 endif |
|
|
94 endif |
|
|
95 |
|
|
96 linespec = {"b"}; |
|
|
97 if (narg <= nargin) |
|
|
98 if (ischar (varargin{narg})) |
|
|
99 linespec = varargin(narg); |
|
|
100 endif |
|
|
101 endif |
|
|
102 |
|
|
103 lx = length (x); |
|
|
104 ly = length (y); |
|
|
105 |
|
|
106 if (lx != ly) |
|
|
107 error ("voronoi: arguments must be vectors of same length"); |
|
|
108 endif |
|
|
109 |
|
6852
|
110 ## Add big box to approximate rays to infinity |
|
|
111 xmax = max (x(:)); |
|
|
112 xmin = min (x(:)); |
|
|
113 ymax = max (y(:)); |
|
|
114 ymin = min (y(:)); |
|
|
115 xdelta = xmax - xmin; |
|
|
116 ydelta = ymax - ymin; |
|
|
117 scale = 10e4; |
|
|
118 |
|
|
119 xbox = [xmin - scale * xdelta; xmin - scale * xdelta; ... |
|
|
120 xmax + scale * xdelta; xmax + scale * xdelta]; |
|
|
121 ybox = [xmin - scale * xdelta; xmax + scale * xdelta; ... |
|
|
122 xmax + scale * xdelta; xmin - scale * xdelta]; |
|
|
123 x = [x(:) ; xbox(:)]; |
|
|
124 y = [y(:) ; ybox(:)]; |
|
|
125 |
|
|
126 [p, c, infi] = __voronoi__ ([x(:), y(:)], opts{:}); |
|
6823
|
127 |
|
|
128 idx = find (!infi); |
|
|
129 |
|
|
130 ll = length (idx); |
|
|
131 k = 0;r = 1; |
|
|
132 |
|
6852
|
133 for i = 1 : ll |
|
6823
|
134 k += length (c{idx(i)}); |
|
|
135 endfor |
|
|
136 |
|
|
137 edges = zeros (2, k); |
|
|
138 |
|
6852
|
139 for i = 1 : ll |
|
6823
|
140 fac = c {idx(i)}; |
|
6852
|
141 lf = length (fac); |
|
6823
|
142 fac = [fac, fac(1)]; |
|
6852
|
143 fst = fac (1 : length(fac) - 1); |
|
|
144 sec = fac(2 : length(fac)); |
|
|
145 edges (:, r : r + lf - 1) = [fst; sec]; |
|
6823
|
146 r += lf; |
|
|
147 endfor |
|
|
148 |
|
|
149 ## Identify the unique edges of the Voronoi diagram |
|
|
150 edges = sortrows (sort (edges).').'; |
|
6852
|
151 edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ... |
|
|
152 edges(2, 1 :end - 1) != edges(2, 2 : end)), true]); |
|
|
153 |
|
|
154 ## Eliminate the edges of the diagram representing the box |
|
|
155 poutside = (1 : rows(p)) ... |
|
|
156 (p (:, 1) < xmin - xdelta | p (:, 1) > xmax + xdelta | ... |
|
|
157 p (:, 2) < ymin - ydelta | p (:, 2) > ymax + ydelta); |
|
|
158 edgeoutside = ismember (edges (1, :), poutside) & ... |
|
|
159 ismember (edges (2, :), poutside); |
|
|
160 edges (:, edgeoutside) = []; |
|
6823
|
161 |
|
|
162 ## Get points of the diagram |
|
6852
|
163 vx = reshape (p (edges, 1), size(edges)); |
|
|
164 vy = reshape (p (edges, 2), size(edges)); |
|
6823
|
165 |
|
|
166 if (nargout < 2) |
|
6852
|
167 lim = [xmin, xmax, ymin, ymax]; |
|
6823
|
168 h = plot (handl, vx, vy, linespec{:}, x, y, '+'); |
|
6852
|
169 axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ... |
|
|
170 [-1, 1] * (lim (4) - lim (3))]); |
|
6823
|
171 if (nargout == 1) |
|
|
172 vxx = h; |
|
|
173 endif |
|
|
174 elseif (nargout == 2) |
|
|
175 vvx = vx; |
|
|
176 vvy = vy; |
|
|
177 else |
|
|
178 error ("voronoi: only two or zero output arguments supported") |
|
|
179 endif |
|
|
180 |
|
|
181 endfunction |
