octave_packages/geometry-1.5.0/polygons2d/medialAxisConvex.m

   1 ## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr>
   2 ## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   3 ## All rights reserved.
   4 ##
   5 ## Redistribution and use in source and binary forms, with or without
   6 ## modification, are permitted provided that the following conditions are met:
   7 ##
   8 ##     1 Redistributions of source code must retain the above copyright notice,
   9 ##       this list of conditions and the following disclaimer.
  10 ##     2 Redistributions in binary form must reproduce the above copyright
  11 ##       notice, this list of conditions and the following disclaimer in the
  12 ##       documentation and/or other materials provided with the distribution.
  13 ##
  14 ## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
  15 ## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16 ## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17 ## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
  18 ## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19 ## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  20 ## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  21 ## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  22 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  23 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  24 ##
  25 ## The views and conclusions contained in the software and documentation are
  26 ## those of the authors and should not be interpreted as representing official
  27 ## policies, either expressed or implied, of the copyright holders.
  28
  29 ## -*- texinfo -*-
  30 ## @deftypefn {Function File} {[@var{n} @var{e}] = } medialAxisConvex (@var{polygon})
  31 ## Compute medial axis of a convex polygon
  32 ##
  33 ##   @var{polygon} is given as a set of points [x1 y1;x2 y2 ...], returns
  34 ##   the medial axis of the polygon as a graph.
  35 ##   @var{n} is a set of nodes. The first elements of @var{n} are the vertices of the
  36 ##   original polygon.
  37 ##   @var{e} is a set of edges, containing indices of source and target nodes.
  38 ##   Edges are sorted according to order of creation. Index of first vertex
  39 ##   is lower than index of last vertex, i.e. edges always point to newly
  40 ##   created nodes.
  41 ##
  42 ##   Notes:
  43 ##   - Is not fully implemented, need more development (usually crashes for
  44 ##       polygons with more than 6-7 points...)
  45 ##   - Works only for convex polygons.
  46 ##   - Complexity is not optimal: this algorithm is O(n*log n), but linear
  47 ##   algorithms exist.
  48 ##
  49 ## @seealso{polygons2d, bisector}
  50 ## @end deftypefn
  51
  52 function [nodes, edges] = medialAxisConvex(points)
  53
  54   # eventually remove the last point if it is the same as the first one
  55   if points(1,:) == points(end, :)
  56       nodes = points(1:end-1, :);
  57   else
  58       nodes = points;
  59   end
  60
  61   # special case of triangles:
  62   # compute directly the gravity center, and simplify computation.
  63   if size(nodes, 1)==3
  64       nodes = [nodes; mean(nodes, 1)];
  65       edges = [1 4;2 4;3 4];
  66       return
  67   end
  68
  69   # number of nodes, and also of initial rays
  70   N = size(nodes, 1);
  71
  72   # create ray of each vertex
  73   rays = zeros(N, 4);
  74   rays(1, 1:4) = bisector(nodes([2 1 N], :));
  75   rays(N, 1:4) = bisector(nodes([1 N N-1], :));
  76   for i=2:N-1
  77       rays(i, 1:4) = bisector(nodes([i+1, i, i-1], :));
  78   end
  79
  80   # add indices of edges producing rays (indices of first vertex, second
  81   # vertex is obtained by adding one modulo N).
  82   rayEdges = [[N (1:N-1)]' (1:N)'];
  83
  84   pint = intersectLines(rays, rays([2:N 1], :));
  85   #ti   = linePosition(pint, rays);
  86   #ti   = min(linePosition(pint, rays), linePosition(pint, rays([2:N 1], :)));
  87   ti = distancePointLine(pint, ...
  88       createLine(points([N (1:N-1)]', :), points((1:N)', :)));
  89
  90   # create list of events.
  91   # terms are : R1 R2 X Y t0
  92   # R1 and R2 are indices of involved rays
  93   # X and Y is coordinate of intersection point
  94   # t0 is position of point on rays
  95   events = sortrows([ (1:N)' [2:N 1]' pint ti], 5);
  96
  97   # initialize edges
  98   edges = zeros(0, 2);
  99
 100
 101   # -------------------
 102   # process each event until there is no more
 103
 104   # start after index of last vertex, and process N-3 intermediate rays
 105   for i=N+1:2*N-3
 106       # add new node at the rays intersection
 107       nodes(i,:) = events(1, 3:4);
 108
 109       # add new couple of edges
 110       edges = [edges; events(1,1) i; events(1,2) i];
 111
 112       # find the two edges creating the new emanating ray
 113       n1 = rayEdges(events(1, 1), 1);
 114       n2 = rayEdges(events(1, 2), 2);
 115
 116       # create the new ray
 117       line1 = createLine(nodes(n1, :), nodes(mod(n1,N)+1, :));
 118       line2 = createLine(nodes(mod(n2,N)+1, :), nodes(n2, :));
 119       ray0 = bisector(line1, line2);
 120
 121       # set its origin to emanating point
 122       ray0(1:2) = nodes(i, :);
 123
 124       # add the new ray to the list
 125       rays = [rays; ray0];
 126       rayEdges(size(rayEdges, 1)+1, 1:2) = [n1 n2];
 127
 128       # find the two neighbour rays
 129       ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2)==0;
 130       ir = unique(events(ind, 1:2));
 131       ir = ir(~ismember(ir, events(1,1:2)));
 132
 133       # create new intersections
 134       pint = intersectLines(ray0, rays(ir, :));
 135       #ti   = min(linePosition(pint, ray0), linePosition(pint, rays(ir, :))) + events(1,5);
 136       ti = distancePointLine(pint, line1);
 137
 138       # remove all events involving old intersected rays
 139       ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2)==0;
 140       events = events(ind, :);
 141
 142       # add the newly formed events
 143       events = [events; ir(1) i pint(1,:) ti(1); ir(2) i pint(2,:) ti(2)];
 144
 145       # and sort them according to 'position' parameter
 146       events = sortrows(events, 5);
 147   end
 148
 149   # centroid computation for last 3 rays
 150   nodes = [nodes; mean(events(:, 3:4))];
 151   edges = [edges; [unique(events(:,1:2)) ones(3, 1)*(2*N-2)]];
 152
 153 endfunction