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)]];