1## Copyright (C) 2021 David Legland 2## All rights reserved. 3## 4## Redistribution and use in source and binary forms, with or without 5## modification, are permitted provided that the following conditions are met: 6## 7## 1 Redistributions of source code must retain the above copyright notice, 8## this list of conditions and the following disclaimer. 9## 2 Redistributions in binary form must reproduce the above copyright 10## notice, this list of conditions and the following disclaimer in the 11## documentation and/or other materials provided with the distribution. 12## 13## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' 14## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 15## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 16## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR 17## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 18## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 19## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 20## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 21## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 22## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 23## 24## The views and conclusions contained in the software and documentation are 25## those of the authors and should not be interpreted as representing official 26## policies, either expressed or implied, of the copyright holders. 27 28function varargout = polygonSkeleton(poly, varargin) 29% Skeletonization of a polygon with a dense distribution of vertices. 30% 31% [V, ADJ] = polygonSkeleton(POLY) 32% POLY is given as a N-by-2 array of polygon vertex coordinates. The 33% result is given a Nv-by-2 array of skeleton vertex coordinates, and an 34% adjacency list, as a NV-by-1 cell array. Each cell contains indices of 35% vertices adjacent to the vertex indexed by the cell. 36% 37% [V, ADJ, RAD] = polygonSkeleton(POLY) 38% Also returns the radius of each vertex, corresponding to the distance 39% between the vertex and the closest point of the original contour 40% polygon. 41% 42% SKEL = polygonSkeleton(POLY) 43% Concatenates the results in a struct SKEL with following fields: 44% * vertices the Nv-by-2 array of skeleton vertex coordinates 45% * adjList the adjacency list of each vertex, as a Nv-by-1 cell array. 46% * radius the Nv-by-1 array of radius for each vertex 47% 48% Example 49% % start from a binary shape 50% img = imread('circles.png'); 51% img = imFillHoles(img); 52% figure; imshow(img); hold on; 53% % compute a smooth contour 54% cntList = imContours(img); 55% cnts = smoothPolygon(cntList{1}, 5); 56% drawPolygon(cnts, 'g'); 57% % compute skeleton 58% [vertices, adjList] = polygonSkeleton(poly); 59% % convert adjacency list to an edge array 60% edges = adjacencyListToEdges(adjList); 61% % draw the skeleton graph 62% drawGraphEdges(vertices, edges); 63% 64% See also 65% graphs, adjacencyListToEdges 66 67% ------ 68% Author: David Legland 69% e-mail: david.legland@inrae.fr 70% INRAE - BIA Research Unit - BIBS Platform (Nantes) 71% Created: 2020-05-29, using Matlab 9.8.0.1323502 (R2020a) 72% Copyright 2020 INRAE. 73 74 75%% Voronoi Diagram computation 76 77% Compute Voronoi Diagram, using polygon vertices as germs. 78[V, C] = voronoin(poly); 79 80% compute number of elements of each array 81nVertices = size(V, 1); 82nCells = size(C, 1); 83 84% Detection of the vertices located inside the contour 85insideFlag = inpolygon(V(:,1), V(:,2), poly(:,1), poly(:,2)); 86innerVertices = V(insideFlag, :); 87 88% indices of inner vertices 89indsInside = find(insideFlag); 90nInnerVertices = length(indsInside); 91 92% compute map between voronoi vertex indices and skeleton vertex indices. 93vertexIndexMap = zeros(nVertices, 1); 94for iVertex = 1:length(indsInside) 95 vertexIndexMap(indsInside(iVertex)) = iVertex; 96end 97 98 99%% Compute the topology of the skeleton 100% 101% Compute the topology as a list of adjacent vertex indices for each vertex 102% inside the polygon. 103% Need to convert between voronoi indices and skeleton indices. 104 105% allocate adjacncy list 106adjList = cell(nInnerVertices, 1); 107vertexGermInds = zeros(nInnerVertices, 1); 108 109% iterate on voronoi cells to compute skeleton by linking adjacent vertices 110% (avoiding first cell which is located at infinity by convention) 111for iGerm = 2:nCells 112 % vertices of current cell 113 cellVertices = C{iGerm}; 114 nCellVertices = length(cellVertices); 115 116 % iterate on vertices of current cell 117 for k = 1:nCellVertices 118 % index of current voronoi vertex 119 iVertex = cellVertices(k); 120 121 % process only vertices within the contour 122 if insideFlag(iVertex) == 0 123 continue; 124 end 125 126 % convert voronoi vertex index to skeleton vertex index 127 indV1 = vertexIndexMap(iVertex); 128 129 % update the reference germ associated to current skeleton vertex 130 vertexGermInds(indV1) = iGerm; 131 132 % compute indices of adjacent vertices (in cell) 133 iPrev = cellVertices(mod(k - 2, nCellVertices) + 1); 134 iNext = cellVertices(mod(k, nCellVertices) + 1); 135 136 % keep only the neighbors within the polygon 137 if insideFlag(iPrev) == 1 138 adjList{indV1} = [adjList{indV1} vertexIndexMap(iPrev)]; 139 end 140 if insideFlag(iNext) == 1 141 adjList{indV1} = [adjList{indV1} vertexIndexMap(iNext)]; 142 end 143 end 144end 145 146% cleanup to avoid duplicate indices 147for iVertex = 1:nInnerVertices 148 adjList{iVertex} = unique(adjList{iVertex}); 149end 150 151 152%% Compute radius list 153 154% for each voronoi vertex inside the polygon, compute the distance to 155% original polygon. 156% Find indices of germs associated to each vertex. 157% By construction, each vertex is the circumcenter of three germs. 158radiusList = zeros(nInnerVertices, 1); 159for iVertex = 1:nInnerVertices 160 radiusList(iVertex) = norm(poly(vertexGermInds(iVertex),:) - innerVertices(iVertex,:)); 161end 162 163 164%% Format output 165 166if nargout <= 1 167 % format output to a struct 168 varargout{1} = struct('vertices', {innerVertices}, 'adjList', {adjList}, 'radius', {radiusList}); 169elseif nargout == 2 170 varargout = {innerVertices, adjList}; 171elseif nargout == 3 172 varargout = {innerVertices, adjList, radiusList}; 173end 174