medialAxisConvex

MEDIALAXISCONVEX Compute medial axis of a convex polygon.

   [N, E] = medialAxisConvex(POLYGON);
   where POLYGON is given as a set of points [x1 y1;x2 y2 ...], returns
   the medial axis of the polygon as a graph.
   N is a set of nodes. The first elements of N are the vertices of the
   original polygon.
   E is a set of edges, containing indices of source and target nodes.
   Edges are sorted according to order of creation. Index of first vertex
   is lower than index of last vertex, i.e. edges always point to newly
   created nodes.

   Notes:
   - Is not fully implemented, need more development (usually crashes for
       polygons with more than 6-7 points...)
   - Works only for convex polygons.
   - Complexity is not optimal: this algorithm is O(n*log n), but linear
   algorithms exist.

   See also
   polygons2d, bisector

0001 function [nodes, edges] = medialAxisConvex(points)
0002 %MEDIALAXISCONVEX Compute medial axis of a convex polygon.
0003 %
0004 %   [N, E] = medialAxisConvex(POLYGON);
0005 %   where POLYGON is given as a set of points [x1 y1;x2 y2 ...], returns
0006 %   the medial axis of the polygon as a graph.
0007 %   N is a set of nodes. The first elements of N are the vertices of the
0008 %   original polygon.
0009 %   E is a set of edges, containing indices of source and target nodes.
0010 %   Edges are sorted according to order of creation. Index of first vertex
0011 %   is lower than index of last vertex, i.e. edges always point to newly
0012 %   created nodes.
0013 %
0014 %   Notes:
0015 %   - Is not fully implemented, need more development (usually crashes for
0016 %       polygons with more than 6-7 points...)
0017 %   - Works only for convex polygons.
0018 %   - Complexity is not optimal: this algorithm is O(n*log n), but linear
0019 %   algorithms exist.
0020 %
0021 %   See also
0022 %   polygons2d, bisector
0023 
0024 % ------
0025 % Author: David Legland
0026 % E-mail: david.legland@inrae.fr
0027 % Created: 2005-07-07
0028 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0029 
0030 % eventually remove the last point if it is the same as the first one
0031 if points(1,:) == points(end, :)
0032     nodes = points(1:end-1, :);
0033 else
0034     nodes = points;
0035 end
0036 
0037 % special case of triangles:
0038 % compute directly the gravity center, and simplify computation.
0039 if size(nodes, 1)==3
0040     nodes = [nodes; mean(nodes, 1)];
0041     edges = [1 4;2 4;3 4];
0042     return
0043 end
0044 
0045 % number of nodes, and also of initial rays
0046 N = size(nodes, 1);
0047 
0048 % create ray of each vertex
0049 rays = zeros(N, 4);
0050 rays(1, 1:4) = bisector(nodes([2 1 N], :));
0051 rays(N, 1:4) = bisector(nodes([1 N N-1], :));
0052 for i=2:N-1
0053     rays(i, 1:4) = bisector(nodes([i+1, i, i-1], :));
0054 end
0055 
0056 % add indices of edges producing rays (indices of first vertex, second
0057 % vertex is obtained by adding one modulo N).
0058 rayEdges = [[N (1:N-1)]' (1:N)'];
0059 
0060 pint = intersectLines(rays, rays([2:N 1], :));
0061 
0062 % compute the distance between each intersection point and the closest
0063 % edge. This distance is used as marker to propagate processing front.
0064 ti = zeros(N, 1);
0065 for i = 1:N
0066     line = createLine(points(mod(i-2, N)+1, :), points(i, :));
0067     ti(i) = abs(distancePointLine(pint(i,:), line));
0068 end
0069 
0070 % create list of events.
0071 % terms are : R1 R2 X Y t0
0072 % R1 and R2 are indices of involved rays
0073 % X and Y is coordinate of intersection point
0074 % t0 is position of point on rays
0075 events = sortrows([ (1:N)' [2:N 1]' pint ti], 5);
0076 
0077 % initialize edges
0078 edges = zeros(0, 2);
0079 
0080 %% process each event until there is no more
0081 
0082 % start after index of last vertex, and process N-3 intermediate rays
0083 for i = N+1:2*N-3
0084     % add new node at the rays intersection
0085     nodes(i,:) = events(1, 3:4);
0086     
0087     % add new couple of edges
0088     edges = [edges; events(1,1) i; events(1,2) i]; %#ok<AGROW>
0089             
0090     % find the two edges creating the new emanating ray
0091     n1 = rayEdges(events(1, 1), 1);
0092     n2 = rayEdges(events(1, 2), 2);    
0093     
0094     % create the new ray
0095     line1 = createLine(nodes(n1, :), nodes(mod(n1,N)+1, :));
0096     line2 = createLine(nodes(mod(n2,N)+1, :), nodes(n2, :));
0097     ray0 = bisector(line1, line2);
0098     
0099     % set its origin to emanating point
0100     ray0(1:2) = nodes(i, :);
0101 
0102     % add the new ray to the list
0103     rays = [rays; ray0]; %#ok<AGROW>
0104     rayEdges(size(rayEdges, 1)+1, 1:2) = [n1 n2];
0105     
0106     % find the two neighbour rays
0107     ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2) ~= 0;
0108     ir = unique(events(ind, 1:2));
0109     ir = ir(~ismember(ir, events(1,1:2)));
0110     
0111     % create new intersections
0112     pint = intersectLines(ray0, rays(ir, :));
0113     ti = abs(distancePointLine(pint, line1));
0114     
0115     % remove all events involving old intersected rays
0116     ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2) == 0;
0117     events = events(ind, :);
0118     
0119     % add the newly formed events
0120     events = [events; ir(1) i pint(1,:) ti(1); ir(2) i pint(2,:) ti(2)]; %#ok<AGROW>
0121 
0122     % and sort them according to 'position' parameter
0123     events = sortrows(events, 5);
0124 end
0125 
0126 % centroid computation for last 3 rays
0127 nodes = [nodes; mean(events(:, 3:4))];
0128 edges = [edges; [unique(events(:,1:2)) ones(3, 1)*(2*N-2)]];