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
---------
author : David Legland
INRA - TPV URPOI - BIA IMASTE
created the 07/07/2005.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 % INRA - TPV URPOI - BIA IMASTE 0027 % created the 07/07/2005. 0028 % 0029 0030 % HISTORY 0031 % 18/04/2007: fix some typos, note the function to be unimplemented 0032 0033 % TODO: is not fully implemented, need to finish it 0034 0035 % eventually remove the last point if it is the same as the first one 0036 if points(1,:) == points(end, :) 0037 nodes = points(1:end-1, :); 0038 else 0039 nodes = points; 0040 end 0041 0042 % special case of triangles: 0043 % compute directly the gravity center, and simplify computation. 0044 if size(nodes, 1)==3 0045 nodes = [nodes; mean(nodes, 1)]; 0046 edges = [1 4;2 4;3 4]; 0047 return 0048 end 0049 0050 % number of nodes, and also of initial rays 0051 N = size(nodes, 1); 0052 0053 % create ray of each vertex 0054 rays = zeros(N, 4); 0055 rays(1, 1:4) = bisector(nodes([2 1 N], :)); 0056 rays(N, 1:4) = bisector(nodes([1 N N-1], :)); 0057 for i=2:N-1 0058 rays(i, 1:4) = bisector(nodes([i+1, i, i-1], :)); 0059 end 0060 0061 % add indices of edges producing rays (indices of first vertex, second 0062 % vertex is obtained by adding one modulo N). 0063 rayEdges = [[N (1:N-1)]' (1:N)']; 0064 0065 pint = intersectLines(rays, rays([2:N 1], :)); 0066 0067 % compute the distance between each intersection point and the closest 0068 % edge. This distance is used as marker to propagate processing front. 0069 ti = zeros(N, 1); 0070 for i = 1:N 0071 line = createLine(points(mod(i-2, N)+1, :), points(i, :)); 0072 ti(i) = abs(distancePointLine(pint(i,:), line)); 0073 end 0074 0075 % create list of events. 0076 % terms are : R1 R2 X Y t0 0077 % R1 and R2 are indices of involved rays 0078 % X and Y is coordinate of intersection point 0079 % t0 is position of point on rays 0080 events = sortrows([ (1:N)' [2:N 1]' pint ti], 5); 0081 0082 % initialize edges 0083 edges = zeros(0, 2); 0084 0085 0086 % ------------------- 0087 % process each event until there is no more 0088 0089 % start after index of last vertex, and process N-3 intermediate rays 0090 for i = N+1:2*N-3 0091 % add new node at the rays intersection 0092 nodes(i,:) = events(1, 3:4); 0093 0094 % add new couple of edges 0095 edges = [edges; events(1,1) i; events(1,2) i]; %#ok<AGROW> 0096 0097 % find the two edges creating the new emanating ray 0098 n1 = rayEdges(events(1, 1), 1); 0099 n2 = rayEdges(events(1, 2), 2); 0100 0101 % create the new ray 0102 line1 = createLine(nodes(n1, :), nodes(mod(n1,N)+1, :)); 0103 line2 = createLine(nodes(mod(n2,N)+1, :), nodes(n2, :)); 0104 ray0 = bisector(line1, line2); 0105 0106 % set its origin to emanating point 0107 ray0(1:2) = nodes(i, :); 0108 0109 % add the new ray to the list 0110 rays = [rays; ray0]; %#ok<AGROW> 0111 rayEdges(size(rayEdges, 1)+1, 1:2) = [n1 n2]; 0112 0113 % find the two neighbour rays 0114 ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2) ~= 0; 0115 ir = unique(events(ind, 1:2)); 0116 ir = ir(~ismember(ir, events(1,1:2))); 0117 0118 % create new intersections 0119 pint = intersectLines(ray0, rays(ir, :)); 0120 ti = abs(distancePointLine(pint, line1)); 0121 0122 % remove all events involving old intersected rays 0123 ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2) == 0; 0124 events = events(ind, :); 0125 0126 % add the newly formed events 0127 events = [events; ir(1) i pint(1,:) ti(1); ir(2) i pint(2,:) ti(2)]; %#ok<AGROW> 0128 0129 % and sort them according to 'position' parameter 0130 events = sortrows(events, 5); 0131 end 0132 0133 % centroid computation for last 3 rays 0134 nodes = [nodes; mean(events(:, 3:4))]; 0135 edges = [edges; [unique(events(:,1:2)) ones(3, 1)*(2*N-2)]];