Skeletonization of a polygon with a dense distribution of vertices.
[V, ADJ] = polygonSkeleton(POLY)
POLY is given as a N-by-2 array of polygon vertex coordinates. The
result is given a Nv-by-2 array of skeleton vertex coordinates, and an
adjacency list, as a NV-by-1 cell array. Each cell contains indices of
vertices adjacent to the vertex indexed by the cell.
[V, ADJ, RAD] = polygonSkeleton(POLY)
Also returns the radius of each vertex, corresponding to the distance
between the vertex and the closest point of the original contour
polygon.
SKEL = polygonSkeleton(POLY)
Concatenates the results in a struct SKEL with following fields:
* vertices the Nv-by-2 array of skeleton vertex coordinates
* adjList the adjacency list of each vertex, as a Nv-by-1 cell array.
* radius the Nv-by-1 array of radius for each vertex
Example
% start from a binary shape
img = imread('circles.png');
img = imFillHoles(img);
figure; imshow(img); hold on;
% compute a smooth contour
cntList = imContours(img);
cnts = smoothPolygon(cntList{1}, 5);
drawPolygon(cnts, 'g');
% compute skeleton
[vertices, adjList] = polygonSkeleton(poly);
% convert adjacency list to an edge array
edges = adjacencyListToEdges(adjList);
% draw the skeleton graph
drawGraphEdges(vertices, edges);
See also
graphs, adjacencyListToEdges0001 function varargout = polygonSkeleton(poly, varargin) 0002 % Skeletonization of a polygon with a dense distribution of vertices. 0003 % 0004 % [V, ADJ] = polygonSkeleton(POLY) 0005 % POLY is given as a N-by-2 array of polygon vertex coordinates. The 0006 % result is given a Nv-by-2 array of skeleton vertex coordinates, and an 0007 % adjacency list, as a NV-by-1 cell array. Each cell contains indices of 0008 % vertices adjacent to the vertex indexed by the cell. 0009 % 0010 % [V, ADJ, RAD] = polygonSkeleton(POLY) 0011 % Also returns the radius of each vertex, corresponding to the distance 0012 % between the vertex and the closest point of the original contour 0013 % polygon. 0014 % 0015 % SKEL = polygonSkeleton(POLY) 0016 % Concatenates the results in a struct SKEL with following fields: 0017 % * vertices the Nv-by-2 array of skeleton vertex coordinates 0018 % * adjList the adjacency list of each vertex, as a Nv-by-1 cell array. 0019 % * radius the Nv-by-1 array of radius for each vertex 0020 % 0021 % Example 0022 % % start from a binary shape 0023 % img = imread('circles.png'); 0024 % img = imFillHoles(img); 0025 % figure; imshow(img); hold on; 0026 % % compute a smooth contour 0027 % cntList = imContours(img); 0028 % cnts = smoothPolygon(cntList{1}, 5); 0029 % drawPolygon(cnts, 'g'); 0030 % % compute skeleton 0031 % [vertices, adjList] = polygonSkeleton(poly); 0032 % % convert adjacency list to an edge array 0033 % edges = adjacencyListToEdges(adjList); 0034 % % draw the skeleton graph 0035 % drawGraphEdges(vertices, edges); 0036 % 0037 % See also 0038 % graphs, adjacencyListToEdges 0039 0040 % ------ 0041 % Author: David Legland 0042 % e-mail: david.legland@inrae.fr 0043 % INRAE - BIA Research Unit - BIBS Platform (Nantes) 0044 % Created: 2020-05-29, using Matlab 9.8.0.1323502 (R2020a) 0045 % Copyright 2020 INRAE. 0046 0047 0048 %% Voronoi Diagram computation 0049 0050 % Compute Voronoi Diagram, using polygon vertices as germs. 0051 [V, C] = voronoin(poly); 0052 0053 % compute number of elements of each array 0054 nVertices = size(V, 1); 0055 nCells = size(C, 1); 0056 0057 % Detection of the vertices located inside the contour 0058 insideFlag = inpolygon(V(:,1), V(:,2), poly(:,1), poly(:,2)); 0059 innerVertices = V(insideFlag, :); 0060 0061 % indices of inner vertices 0062 indsInside = find(insideFlag); 0063 nInnerVertices = length(indsInside); 0064 0065 % compute map between voronoi vertex indices and skeleton vertex indices. 0066 vertexIndexMap = zeros(nVertices, 1); 0067 for iVertex = 1:length(indsInside) 0068 vertexIndexMap(indsInside(iVertex)) = iVertex; 0069 end 0070 0071 0072 %% Compute the topology of the skeleton 0073 % 0074 % Compute the topology as a list of adjacent vertex indices for each vertex 0075 % inside the polygon. 0076 % Need to convert between voronoi indices and skeleton indices. 0077 0078 % allocate adjacncy list 0079 adjList = cell(nInnerVertices, 1); 0080 vertexGermInds = zeros(nInnerVertices, 1); 0081 0082 % iterate on voronoi cells to compute skeleton by linking adjacent vertices 0083 % (avoiding first cell which is located at infinity by convention) 0084 for iGerm = 2:nCells 0085 % vertices of current cell 0086 cellVertices = C{iGerm}; 0087 nCellVertices = length(cellVertices); 0088 0089 % iterate on vertices of current cell 0090 for k = 1:nCellVertices 0091 % index of current voronoi vertex 0092 iVertex = cellVertices(k); 0093 0094 % process only vertices within the contour 0095 if insideFlag(iVertex) == 0 0096 continue; 0097 end 0098 0099 % convert voronoi vertex index to skeleton vertex index 0100 indV1 = vertexIndexMap(iVertex); 0101 0102 % update the reference germ associated to current skeleton vertex 0103 vertexGermInds(indV1) = iGerm; 0104 0105 % compute indices of adjacent vertices (in cell) 0106 iPrev = cellVertices(mod(k - 2, nCellVertices) + 1); 0107 iNext = cellVertices(mod(k, nCellVertices) + 1); 0108 0109 % keep only the neighbors within the polygon 0110 if insideFlag(iPrev) == 1 0111 adjList{indV1} = [adjList{indV1} vertexIndexMap(iPrev)]; 0112 end 0113 if insideFlag(iNext) == 1 0114 adjList{indV1} = [adjList{indV1} vertexIndexMap(iNext)]; 0115 end 0116 end 0117 end 0118 0119 % cleanup to avoid duplicate indices 0120 for iVertex = 1:nInnerVertices 0121 adjList{iVertex} = unique(adjList{iVertex}); 0122 end 0123 0124 0125 %% Compute radius list 0126 0127 % for each voronoi vertex inside the polygon, compute the distance to 0128 % original polygon. 0129 % Find indices of germs associated to each vertex. 0130 % By construction, each vertex is the circumcenter of three germs. 0131 radiusList = zeros(nInnerVertices, 1); 0132 for iVertex = 1:nInnerVertices 0133 radiusList(iVertex) = norm(poly(vertexGermInds(iVertex),:) - innerVertices(iVertex,:)); 0134 end 0135 0136 0137 %% Format output 0138 0139 if nargout <= 1 0140 % format output to a struct 0141 varargout{1} = struct('vertices', {innerVertices}, 'adjList', {adjList}, 'radius', {radiusList}); 0142 elseif nargout == 2 0143 varargout = {innerVertices, adjList}; 0144 elseif nargout == 3 0145 varargout = {innerVertices, adjList, radiusList}; 0146 end