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