Mercurial > hg > octave-avbm
changeset 6852:a34d59fc7a91
[project @ 2007-08-31 20:51:58 by dbateman]
author | dbateman |
---|---|
date | Fri, 31 Aug 2007 20:51:58 +0000 |
parents | 05f6f120a65f |
children | 02106856b8f6 |
files | scripts/ChangeLog scripts/geometry/voronoi.m |
diffstat | 2 files changed, 43 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2007-08-31 David Bateman <dbateman@free.fr> + + * geometry/voronoi.m: Add large box around data to get a good + approximation of the rays to infinity. + 2007-08-31 Michael goffioul <michael.goffioul@gmail.com> * plot/axes.m: Allow parent to be specified when creating axes
--- a/scripts/geometry/voronoi.m +++ b/scripts/geometry/voronoi.m @@ -107,42 +107,67 @@ error ("voronoi: arguments must be vectors of same length"); endif - [p, c, infi] = __voronoi__ ([x(:),y(:)], opts{:}); + ## Add big box to approximate rays to infinity + xmax = max (x(:)); + xmin = min (x(:)); + ymax = max (y(:)); + ymin = min (y(:)); + xdelta = xmax - xmin; + ydelta = ymax - ymin; + scale = 10e4; + + xbox = [xmin - scale * xdelta; xmin - scale * xdelta; ... + xmax + scale * xdelta; xmax + scale * xdelta]; + ybox = [xmin - scale * xdelta; xmax + scale * xdelta; ... + xmax + scale * xdelta; xmin - scale * xdelta]; + x = [x(:) ; xbox(:)]; + y = [y(:) ; ybox(:)]; + + [p, c, infi] = __voronoi__ ([x(:), y(:)], opts{:}); idx = find (!infi); ll = length (idx); k = 0;r = 1; - for i = 1:ll + for i = 1 : ll k += length (c{idx(i)}); endfor edges = zeros (2, k); - for i=1:ll + for i = 1 : ll fac = c {idx(i)}; - lf = length(fac); + lf = length (fac); fac = [fac, fac(1)]; - fst = fac(1:length(fac) - 1); - sec = fac(2:length(fac)); - edges(:,r:r+lf-1) = [fst; sec]; + fst = fac (1 : length(fac) - 1); + sec = fac(2 : length(fac)); + edges (:, r : r + lf - 1) = [fst; sec]; r += lf; endfor ## Identify the unique edges of the Voronoi diagram edges = sortrows (sort (edges).').'; - edges = edges (:,[(edges(1,1:end-1) != edges(1,2:end) | ... - edges(2,1:end-1) != edges(2,2:end)), true]); + edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ... + edges(2, 1 :end - 1) != edges(2, 2 : end)), true]); + + ## Eliminate the edges of the diagram representing the box + poutside = (1 : rows(p)) ... + (p (:, 1) < xmin - xdelta | p (:, 1) > xmax + xdelta | ... + p (:, 2) < ymin - ydelta | p (:, 2) > ymax + ydelta); + edgeoutside = ismember (edges (1, :), poutside) & ... + ismember (edges (2, :), poutside); + edges (:, edgeoutside) = []; ## Get points of the diagram - vx = reshape(p(edges, 1), size(edges)); - vy = reshape(p(edges, 2), size(edges)); + vx = reshape (p (edges, 1), size(edges)); + vy = reshape (p (edges, 2), size(edges)); if (nargout < 2) - lim = [min(x(:)), max(x(:)), min(y(:)), max(y(:))]; + lim = [xmin, xmax, ymin, ymax]; h = plot (handl, vx, vy, linespec{:}, x, y, '+'); - axis(lim+0.1*[[-1,1]*(lim(2)-lim(1)),[-1,1]*(lim(4)-lim(3))]); + axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ... + [-1, 1] * (lim (4) - lim (3))]); if (nargout == 1) vxx = h; endif