Mercurial > hg > octave-nkf

comparison scripts/geometry/voronoi.m @ 8440:e792c736b1ac

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix for dense grids and speed up for voronoi function
author David Bateman <dbateman@free.fr>
date Mon, 29 Dec 2008 23:41:18 +0100
parents fa78cb8d8a5c
children e07e93c04080
comparison
equal deleted inserted replaced
8439:a6b4d8fdbea1 8440:e792c736b1ac
104 104
105 if (lx != ly) 105 if (lx != ly)
106 error ("voronoi: arguments must be vectors of same length"); 106 error ("voronoi: arguments must be vectors of same length");
107 endif 107 endif
108 108
109 ## Add big box to approximate rays to infinity 109 ## Add box to approximate rays to infinity. For Voronoi diagrams the
110 ## box can (and should) be close to the points themselves. To make the
111 ## job of finding the exterior edges it should be at least two times the
112 ## delta below however
110 xmax = max (x(:)); 113 xmax = max (x(:));
111 xmin = min (x(:)); 114 xmin = min (x(:));
112 ymax = max (y(:)); 115 ymax = max (y(:));
113 ymin = min (y(:)); 116 ymin = min (y(:));
114 xdelta = xmax - xmin; 117 xdelta = xmax - xmin;
115 ydelta = ymax - ymin; 118 ydelta = ymax - ymin;
116 scale = 10e4; 119 scale = 2;
117 120
118 xbox = [xmin - scale * xdelta; xmin - scale * xdelta; ... 121 xbox = [xmin - scale * xdelta; xmin - scale * xdelta; ...
119 xmax + scale * xdelta; xmax + scale * xdelta]; 122 xmax + scale * xdelta; xmax + scale * xdelta];
120 ybox = [xmin - scale * xdelta; xmax + scale * xdelta; ... 123 ybox = [xmin - scale * xdelta; xmax + scale * xdelta; ...
121 xmax + scale * xdelta; xmin - scale * xdelta]; 124 xmax + scale * xdelta; xmin - scale * xdelta];
122 x = [x(:) ; xbox(:)];
123 y = [y(:) ; ybox(:)];
124 125
125 [p, c, infi] = __voronoi__ ([x(:), y(:)], opts{:}); 126 [p, c, infi] = __voronoi__ ([[x(:) ; xbox(:)], [y(:); ybox(:)]], opts{:});
126 127
127 idx = find (!infi); 128 idx = find (!infi);
128
129 ll = length (idx); 129 ll = length (idx);
130 k = 0;r = 1; 130 c = c(idx).';
131 131 k = sum (cellfun ('length', c));
132 for i = 1 : ll 132 edges = cell2mat(cellfun (@(x) [x ; [x(end), x(1:end-1)]], c,
133 k += length (c{idx(i)}); 133 "UniformOutput", false));
134 endfor
135
136 edges = zeros (2, k);
137
138 for i = 1 : ll
139 fac = c{idx(i)};
140 lf = length (fac);
141 fac = [fac, fac(1)];
142 fst = fac (1 : length(fac) - 1);
143 sec = fac(2 : length(fac));
144 edges (:, r : r + lf - 1) = [fst; sec];
145 r += lf;
146 endfor
147 134
148 ## Identify the unique edges of the Voronoi diagram 135 ## Identify the unique edges of the Voronoi diagram
149 edges = sortrows (sort (edges).').'; 136 edges = sortrows (sort (edges).').';
150 edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ... 137 edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ...
151 edges(2, 1 :end - 1) != edges(2, 2 : end)), true]); 138 edges(2, 1 :end - 1) != edges(2, 2 : end)), true]);