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

0001 function [regions, distances] = meshVoronoiDiagram(vertices, faces, germInds)
0002 % Voronoi Diagram on the surface of a polygonal mesh.
0003 %
0004 %   REGIONS = meshVoronoiDiagram(V, F, GERM_INDS)
0005 %   Computes the region associated to each vertex of the input mesh (V,F),
0006 %   i.e. the index of the germ closest to the vertex.
0007 %   V is a NV-by-3 array of vertx coordinates, and F is a NF-by-3 array
0008 %   containing vertex indices of each face.
0009 %   IGERMS is an array of vertex indices.
0010 %   REGIONS is a column vector with as many rows as the number of vertices,
0011 %   containing for each vertex the index of the closest germ.
0012 %   The regions are computed by propagating distances along edges.
0013 %
0014 %   [REGIONS, DISTANCES] = meshVoronoiDiagram(V, F, GERM_INDS)
0015 %   Also returns the (geodesic) distance from each vertex to the closest
0016 %   germ.
0017 %
0018 %
0019 %   Example
0020 %     % Create a triangular mesh based on an icosahedron shape
0021 %     [v, f] = createIcosahedron; v = v - mean(v);
0022 %     [v, f] = subdivideMesh(v, f, 10); v = normalizeVector3d(v);
0023 %     figure; hold on; axis equal; view(3);
0024 %     drawMesh(v, f, 'faceColor', [.7 .7 .7]);
0025 %     % generate germs within the mesh (identify with indices)
0026 %     nGerms = 10;
0027 %     inds = randperm(size(v, 1), nGerms);
0028 %     drawPoint3d(v(inds,:), 'bo');
0029 %     % Compute Voronoi Diagram
0030 %     [regions, distances] = meshVoronoiDiagram(v, f, inds);
0031 %     % Display regions
0032 %     figure; hold on; axis equal; view(3);
0033 %     cmap = jet(nGerms);
0034 %     patch('Vertices', v, 'Faces', f, 'FaceVertexCData', cmap(regions, :), 'FaceColor', 'interp', 'LineStyle', 'none');
0035 %     % Display distance maps
0036 %     figure; hold on; axis equal; view(3);
0037 %     patch('Vertices', v, 'Faces', f, 'FaceVertexCData', distances, 'FaceColor', 'interp', 'LineStyle', 'none');
0038 %
0039 %   See also
0040 %     meshes3d, drawMesh
0041  
0042 % ------
0043 % Author: David Legland
0044 % e-mail: david.legland@inrae.fr
0045 % INRAE - BIA Research Unit - BIBS Platform (Nantes)
0046 % Created: 2020-04-16,    using Matlab 9.7.0.1247435 (R2019b) Update 2
0047 % Copyright 2020 INRAE.
0048 
0049 % choose initial germ indices if input is a scalar
0050 nVertices = size(vertices, 1);
0051 if isscalar(germInds)
0052     germInds = randperm(nVertices, germInds);
0053 end
0054 nGerms = length(germInds);
0055 
0056 % init info for vertices
0057 distances = inf * ones(nVertices, 1);
0058 regions = zeros(nVertices, 1);
0059 
0060 % initialize vertex data for each germ
0061 for iGerm = 1:nGerms
0062     ind = germInds(iGerm);
0063     distances(ind) = 0;
0064     regions(ind) = iGerm;
0065 end
0066 
0067 % compute the adjacencey matrix, as a Nv-by-Nv sparse matrix of 0 and 1.
0068 adj = meshAdjacencyMatrix(faces);
0069 
0070 
0071 % create the queue of vertex to process, initialized with germ indices.
0072 vertexQueue = germInds;
0073 processed = false(nVertices, 1);
0074 
0075 % process vertices in the queue, using distance as priority
0076 while ~isempty(vertexQueue)
0077     % find the vertex with smallest distance
0078     [~, ind] = min(distances(vertexQueue));
0079     iVertex = vertexQueue(ind);
0080     vertexQueue(ind) = [];
0081     
0082     % info for current vertex
0083     vertex = vertices(iVertex, :);
0084     dist0 = distances(iVertex);
0085     
0086     % neighbors of current vertex
0087     neighbors = find(adj(iVertex,:));
0088     
0089     % keep only vertices not yet processed
0090     neighbors = neighbors(~processed(neighbors));
0091     
0092     for iNeigh = 1:length(neighbors)
0093         indNeigh = neighbors(iNeigh);
0094         dist = dist0 + distancePoints3d(vertex, vertices(indNeigh, :));
0095         
0096         if dist < distances(indNeigh)
0097             distances(indNeigh) = dist;
0098             regions(indNeigh) = regions(iVertex);
0099             vertexQueue = [vertexQueue indNeigh]; %#ok<AGROW>
0100         end
0101     end
0102     
0103     % mark current vertex as processed
0104     processed(iVertex) = true;
0105 end