File: meshVoronoiDiagram.m

package info (click to toggle)
octave-matgeom 1.2.4-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 3,584 kB
  • sloc: objc: 469; makefile: 10
file content (131 lines) | stat: -rw-r--r-- 5,268 bytes parent folder | download
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