smoothMesh

SMOOTHMESH Smooth mesh by replacing each vertex by the average of its neighbors.

   V2 = smoothMesh(V, F)
   [V2, F2] = smoothMesh(V, F)
   Performs smoothing of the values given in V, by using adjacency
   information given in F. 
   V is a numeric array representing either vertex coordinate, or value
   field associated to each vertex. F is an array of faces, given either
   as a NF-by-3 or NF-by-4 numeric array, or as a cell array. 
   Artifact adjacencies are added if faces have more than 4 vertices.

   ... = smoothMesh(V, F, NITER)
   Repeat the smoothing procedure NITER times. This is equivalent to
   calling the smoothMesh function NITER times.


   Example
     [v f] = torusMesh([50 50 50 30 10 30 45]);
     v = v + randn(size(v));
     [v2 f] = smoothMesh(v, f, 3);
     figure; drawMesh(v2, f);
     l = light; lighting gouraud

   See also
     meshes3d, meshAdjacencyMatrix, triangulateFaces, drawMesh

0001 function varargout = smoothMesh(varargin)
0002 %SMOOTHMESH Smooth mesh by replacing each vertex by the average of its neighbors.
0003 %
0004 %   V2 = smoothMesh(V, F)
0005 %   [V2, F2] = smoothMesh(V, F)
0006 %   Performs smoothing of the values given in V, by using adjacency
0007 %   information given in F.
0008 %   V is a numeric array representing either vertex coordinate, or value
0009 %   field associated to each vertex. F is an array of faces, given either
0010 %   as a NF-by-3 or NF-by-4 numeric array, or as a cell array.
0011 %   Artifact adjacencies are added if faces have more than 4 vertices.
0012 %
0013 %   ... = smoothMesh(V, F, NITER)
0014 %   Repeat the smoothing procedure NITER times. This is equivalent to
0015 %   calling the smoothMesh function NITER times.
0016 %
0017 %
0018 %   Example
0019 %     [v f] = torusMesh([50 50 50 30 10 30 45]);
0020 %     v = v + randn(size(v));
0021 %     [v2 f] = smoothMesh(v, f, 3);
0022 %     figure; drawMesh(v2, f);
0023 %     l = light; lighting gouraud
0024 %
0025 %   See also
0026 %     meshes3d, meshAdjacencyMatrix, triangulateFaces, drawMesh
0027 %
0028 
0029 % ------
0030 % Author: David Legland
0031 % e-mail: david.legland@inra.fr
0032 % Created: 2013-04-29,    using Matlab 7.9.0.529 (R2009b)
0033 % Copyright 2013 INRA - Cepia Software Platform.
0034 
0035 var1 = varargin{1};
0036 if isstruct(var1)
0037     vertices = var1.vertices;
0038     faces = var1.faces;
0039     varargin(1) = [];
0040 else
0041     vertices = varargin{1};
0042     faces = varargin{2};
0043     varargin(1:2) = [];
0044 end
0045 
0046 % determine number of iterations
0047 nIter = 1;
0048 if ~isempty(varargin)
0049     nIter = varargin{1};
0050 end
0051 
0052 % compute adjacency matrix,
0053 % result is a Nv-by-Nv matrix with zeros on the diagonal
0054 adj = meshAdjacencyMatrix(faces);
0055 
0056 % ensure the size of the matrix is Nv-by-Nv
0057 % (this can not be the case if some vertices are not referenced)
0058 nv = size(vertices, 1);
0059 if size(adj, 1) < nv
0060     adj(nv, nv) = 0;
0061 end
0062 
0063 % Add "self adjacencies"
0064 adj = adj + speye(nv);
0065 
0066 % weight each vertex by the number of its neighbors
0067 w = spdiags(full(sum(adj, 2).^(-1)), 0, nv, nv);
0068 adj = w * adj;
0069 
0070 % do averaging to smooth the field
0071 v2 = vertices;
0072 for k = 1:nIter
0073     v2 = adj * v2;
0074 end
0075 
0076 varargout = formatMeshOutput(nargout, v2, faces);
0077 
0078 %% Old version
0079 % % Compute vertex adjacencies
0080 % edges = computeMeshEdges(faces);
0081 % v2 = zeros(size(vertices));
0082 %
0083 % % apply several smoothing
0084 % for iter = 1:nIter
0085 %
0086 %     % replace the coords of each vertex by the average coordinate in the
0087 %     % neighborhood
0088 %     for i = 1:size(vertices, 1)
0089 %         edgeInds = sum(edges == i, 2) > 0;
0090 %         neighInds = unique(edges(edgeInds, :));
0091 %         v2(i, :) = mean(vertices(neighInds, :));
0092 %     end
0093 %
0094 %     % update for next iteration
0095 %     vertices = v2;
0096 % end