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.
5 ## Redistribution and use in source and binary forms, with or without
6 ## modification, are permitted provided that the following conditions are met:
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.
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.
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.
30 ## @deftypefn {Function File} {[@var{n} @var{e}] = } medialAxisConvex (@var{polygon})
31 ## Compute medial axis of a convex polygon
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
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
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
49 ## @seealso{polygons2d, bisector}
52 function [nodes, edges] = medialAxisConvex(points)
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, :);
61 # special case of triangles:
62 # compute directly the gravity center, and simplify computation.
64 nodes = [nodes; mean(nodes, 1)];
65 edges = [1 4;2 4;3 4];
69 # number of nodes, and also of initial rays
72 # create ray of each vertex
74 rays(1, 1:4) = bisector(nodes([2 1 N], :));
75 rays(N, 1:4) = bisector(nodes([1 N N-1], :));
77 rays(i, 1:4) = bisector(nodes([i+1, i, i-1], :));
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)'];
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)', :)));
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);
101 # -------------------
102 # process each event until there is no more
104 # start after index of last vertex, and process N-3 intermediate rays
106 # add new node at the rays intersection
107 nodes(i,:) = events(1, 3:4);
109 # add new couple of edges
110 edges = [edges; events(1,1) i; events(1,2) i];
112 # find the two edges creating the new emanating ray
113 n1 = rayEdges(events(1, 1), 1);
114 n2 = rayEdges(events(1, 2), 2);
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);
121 # set its origin to emanating point
122 ray0(1:2) = nodes(i, :);
124 # add the new ray to the list
126 rayEdges(size(rayEdges, 1)+1, 1:2) = [n1 n2];
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)));
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);
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, :);
142 # add the newly formed events
143 events = [events; ir(1) i pint(1,:) ti(1); ir(2) i pint(2,:) ti(2)];
145 # and sort them according to 'position' parameter
146 events = sortrows(events, 5);
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)]];