Mercurial > hg > octave-nkf

comparison scripts/geometry/voronoi.m @ 6852:a34d59fc7a91

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-08-31 20:51:58 by dbateman]
author dbateman
date Fri, 31 Aug 2007 20:51:58 +0000
parents 8618f29520c6
children 93c65f2a5668
comparison
equal deleted inserted replaced
6851:05f6f120a65f 6852:a34d59fc7a91
105 105
106 if (lx != ly) 106 if (lx != ly)
107 error ("voronoi: arguments must be vectors of same length"); 107 error ("voronoi: arguments must be vectors of same length");
108 endif 108 endif
109 109
110 [p, c, infi] = __voronoi__ ([x(:),y(:)], opts{:}); 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{:});
111 127
112 idx = find (!infi); 128 idx = find (!infi);
113 129
114 ll = length (idx); 130 ll = length (idx);
115 k = 0;r = 1; 131 k = 0;r = 1;
116 132
117 for i = 1:ll 133 for i = 1 : ll
118 k += length (c{idx(i)}); 134 k += length (c{idx(i)});
119 endfor 135 endfor
120 136
121 edges = zeros (2, k); 137 edges = zeros (2, k);
122 138
123 for i=1:ll 139 for i = 1 : ll
124 fac = c {idx(i)}; 140 fac = c {idx(i)};
125 lf = length(fac); 141 lf = length (fac);
126 fac = [fac, fac(1)]; 142 fac = [fac, fac(1)];
127 fst = fac(1:length(fac) - 1); 143 fst = fac (1 : length(fac) - 1);
128 sec = fac(2:length(fac)); 144 sec = fac(2 : length(fac));
129 edges(:,r:r+lf-1) = [fst; sec]; 145 edges (:, r : r + lf - 1) = [fst; sec];
130 r += lf; 146 r += lf;
131 endfor 147 endfor
132 148
133 ## Identify the unique edges of the Voronoi diagram 149 ## Identify the unique edges of the Voronoi diagram
134 edges = sortrows (sort (edges).').'; 150 edges = sortrows (sort (edges).').';
135 edges = edges (:,[(edges(1,1:end-1) != edges(1,2:end) | ... 151 edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ...
136 edges(2,1:end-1) != edges(2,2:end)), true]); 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) = [];
137 161
138 ## Get points of the diagram 162 ## Get points of the diagram
139 vx = reshape(p(edges, 1), size(edges)); 163 vx = reshape (p (edges, 1), size(edges));
140 vy = reshape(p(edges, 2), size(edges)); 164 vy = reshape (p (edges, 2), size(edges));
141 165
142 if (nargout < 2) 166 if (nargout < 2)
143 lim = [min(x(:)), max(x(:)), min(y(:)), max(y(:))]; 167 lim = [xmin, xmax, ymin, ymax];
144 h = plot (handl, vx, vy, linespec{:}, x, y, '+'); 168 h = plot (handl, vx, vy, linespec{:}, x, y, '+');
145 axis(lim+0.1*[[-1,1]*(lim(2)-lim(1)),[-1,1]*(lim(4)-lim(3))]); 169 axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ...
170 [-1, 1] * (lim (4) - lim (3))]);
146 if (nargout == 1) 171 if (nargout == 1)
147 vxx = h; 172 vxx = h;
148 endif 173 endif
149 elseif (nargout == 2) 174 elseif (nargout == 2)
150 vvx = vx; 175 vvx = vx;