1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
|
## Copyright (C) 2024 David Legland
## All rights reserved.
##
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are met:
##
## 1 Redistributions of source code must retain the above copyright notice,
## this list of conditions and the following disclaimer.
## 2 Redistributions in binary form must reproduce the above copyright
## notice, this list of conditions and the following disclaimer in the
## documentation and/or other materials provided with the distribution.
##
## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
##
## The views and conclusions contained in the software and documentation are
## those of the authors and should not be interpreted as representing official
## policies, either expressed or implied, of the copyright holders.
function [regions, distances] = meshVoronoiDiagram(vertices, faces, germInds)
%MESHVORONOIDIAGRAM Voronoi Diagram on the surface of a polygonal mesh.
%
% REGIONS = meshVoronoiDiagram(V, F, GERM_INDS)
% Computes the region associated to each vertex of the input mesh (V,F),
% i.e. the index of the germ closest to the vertex.
% V is a NV-by-3 array of vertx coordinates, and F is a NF-by-3 array
% containing vertex indices of each face.
% IGERMS is an array of vertex indices.
% REGIONS is a column vector with as many rows as the number of vertices,
% containing for each vertex the index of the closest germ.
% The regions are computed by propagating distances along edges.
%
% [REGIONS, DISTANCES] = meshVoronoiDiagram(V, F, GERM_INDS)
% Also returns the (geodesic) distance from each vertex to the closest
% germ.
%
%
% Example
% % Create a triangular mesh based on an icosahedron shape
% [v, f] = createIcosahedron; v = v - mean(v);
% [v, f] = subdivideMesh(v, f, 10); v = normalizeVector3d(v);
% figure; hold on; axis equal; view(3);
% drawMesh(v, f, 'faceColor', [.7 .7 .7]);
% % generate germs within the mesh (identify with indices)
% nGerms = 10;
% inds = randperm(size(v, 1), nGerms);
% drawPoint3d(v(inds,:), 'bo');
% % Compute Voronoi Diagram
% [regions, distances] = meshVoronoiDiagram(v, f, inds);
% % Display regions
% figure; hold on; axis equal; view(3);
% cmap = jet(nGerms);
% patch('Vertices', v, 'Faces', f, 'FaceVertexCData', cmap(regions, :), 'FaceColor', 'interp', 'LineStyle', 'none');
% % Display distance maps
% figure; hold on; axis equal; view(3);
% patch('Vertices', v, 'Faces', f, 'FaceVertexCData', distances, 'FaceColor', 'interp', 'LineStyle', 'none');
%
% See also
% meshes3d, drawMesh
% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2020-04-16, using Matlab 9.7.0.1247435 (R2019b) Update 2
% Copyright 2020-2023 INRAE - BIA Research Unit - BIBS Platform (Nantes)
% choose initial germ indices if input is a scalar
nVertices = size(vertices, 1);
if isscalar(germInds)
germInds = randperm(nVertices, germInds);
end
nGerms = length(germInds);
% init info for vertices
distances = inf * ones(nVertices, 1);
regions = zeros(nVertices, 1);
% initialize vertex data for each germ
for iGerm = 1:nGerms
ind = germInds(iGerm);
distances(ind) = 0;
regions(ind) = iGerm;
end
% compute the adjacencey matrix, as a Nv-by-Nv sparse matrix of 0 and 1.
adj = meshAdjacencyMatrix(faces);
% create the queue of vertex to process, initialized with germ indices.
vertexQueue = germInds;
processed = false(nVertices, 1);
% process vertices in the queue, using distance as priority
while ~isempty(vertexQueue)
% find the vertex with smallest distance
[~, ind] = min(distances(vertexQueue));
iVertex = vertexQueue(ind);
vertexQueue(ind) = [];
% info for current vertex
vertex = vertices(iVertex, :);
dist0 = distances(iVertex);
% neighbors of current vertex
neighbors = find(adj(iVertex,:));
% keep only vertices not yet processed
neighbors = neighbors(~processed(neighbors));
for iNeigh = 1:length(neighbors)
indNeigh = neighbors(iNeigh);
dist = dist0 + distancePoints3d(vertex, vertices(indNeigh, :));
if dist < distances(indNeigh)
distances(indNeigh) = dist;
regions(indNeigh) = regions(iVertex);
vertexQueue = [vertexQueue indNeigh]; %#ok<AGROW>
end
end
% mark current vertex as processed
processed(iVertex) = true;
end
|