MERGECOPLANARFACES Merge coplanar faces of a polyhedral mesh.
[NODES, FACES] = mergeCoplanarFaces(NODES, FACES)
[NODES, EDGES, FACES] = mergeCoplanarFaces(NODES, EDGES, FACES)
NODES is a set of 3D points (as a nNodes-by-3 array),
and FACES is one of:
- a nFaces-by-X array containing vertex indices of each face, with each
face having the same number of vertices,
- a nFaces-by-1 cell array, each cell containing indices of a face.
The function groups faces which are coplanar and contiguous, resulting
in a "lighter" mesh. This can be useful for visualizing binary 3D
images for example.
FACES = mergeCoplanarFaces(..., PRECISION)
Adjust the threshold for deciding if two faces are coplanar or
parallel. Default value is 1e-5.
Example
[v, e, f] = createCube;
figure; drawMesh(v, f); view(3); axis equal;
[v2, f2] = mergeCoplanarFaces(v, f);
figure; drawMesh(v2, f2);
view(3); axis equal; view(3);
See also
meshes3d, drawMesh, minConvexHull, triangulateFaces0001 function varargout = mergeCoplanarFaces(nodes, varargin) 0002 %MERGECOPLANARFACES Merge coplanar faces of a polyhedral mesh. 0003 % 0004 % [NODES, FACES] = mergeCoplanarFaces(NODES, FACES) 0005 % [NODES, EDGES, FACES] = mergeCoplanarFaces(NODES, EDGES, FACES) 0006 % NODES is a set of 3D points (as a nNodes-by-3 array), 0007 % and FACES is one of: 0008 % - a nFaces-by-X array containing vertex indices of each face, with each 0009 % face having the same number of vertices, 0010 % - a nFaces-by-1 cell array, each cell containing indices of a face. 0011 % The function groups faces which are coplanar and contiguous, resulting 0012 % in a "lighter" mesh. This can be useful for visualizing binary 3D 0013 % images for example. 0014 % 0015 % FACES = mergeCoplanarFaces(..., PRECISION) 0016 % Adjust the threshold for deciding if two faces are coplanar or 0017 % parallel. Default value is 1e-5. 0018 % 0019 % Example 0020 % [v, e, f] = createCube; 0021 % figure; drawMesh(v, f); view(3); axis equal; 0022 % [v2, f2] = mergeCoplanarFaces(v, f); 0023 % figure; drawMesh(v2, f2); 0024 % view(3); axis equal; view(3); 0025 % 0026 % See also 0027 % meshes3d, drawMesh, minConvexHull, triangulateFaces 0028 % 0029 0030 % ------ 0031 % Author: David Legland 0032 % e-mail: david.legland@inra.fr 0033 % Created: 2006-07-05 0034 % Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas). 0035 0036 % 20/07/2006 add tolerance for coplanarity test 0037 % 21/08/2006 fix small bug due to difference of methods to test 0038 % coplanaritity, sometimes resulting in 3 points of a face not coplanar ! 0039 % Also add control on precision 0040 % 14/08/2007 rename minConvexHull->meshReduce, and extend to non convex 0041 % shapes 0042 % 2011-01-14 code clean up 0043 % 2013-02-22 rename from meshReduce to mergeCoplanarFaces 0044 0045 0046 %% Process input arguments 0047 0048 % set up precision 0049 acc = 1e-5; 0050 if ~isempty(varargin) 0051 var = varargin{end}; 0052 if length(var) == 1 0053 acc = var; 0054 varargin(end) = []; 0055 end 0056 end 0057 0058 % extract faces and edges 0059 if length(varargin) == 1 0060 faces = varargin{1}; 0061 else 0062 faces = varargin{2}; 0063 end 0064 0065 0066 %% Initialisations 0067 0068 % number of faces 0069 nNodes = size(nodes, 1); 0070 nFaces = size(faces, 1); 0071 0072 % compute number of vertices of each face 0073 Fn = ones(nFaces, 1) * size(faces, 2); 0074 0075 % compute normal of each faces 0076 normals = meshFaceNormals(nodes, faces); 0077 0078 % initialize empty faces and edges 0079 faces2 = cell(0, 1); 0080 edges2 = zeros(0, 2); 0081 0082 % Processing flag for each face 0083 % 1: face to process, 0: already processed 0084 % in the beginning, every triangle face need to be processed 0085 flag = ones(nFaces, 1); 0086 0087 0088 %% Main iteration 0089 0090 % iterate on each face 0091 for iFace = 1:nFaces 0092 0093 % check if face was already performed 0094 if ~flag(iFace) 0095 continue; 0096 end 0097 0098 % indices of faces with same normal 0099 ind = find(vectorNorm3d(crossProduct3d(normals(iFace, :), normals)) < acc); 0100 0101 % keep only coplanar faces (test coplanarity of points in both face) 0102 ind2 = false(size(ind)); 0103 for j = 1:length(ind) 0104 ind2(j) = isCoplanar(nodes([faces(iFace,:) faces(ind(j),:)], :), acc); 0105 end 0106 ind2 = ind(ind2); 0107 0108 % compute edges of all faces in the plane 0109 planeEdges = zeros(sum(Fn(ind2)), 2); 0110 pos = 1; 0111 for i = 1:length(ind2) 0112 face = faces(ind2(i), :); 0113 faceEdges = sort([face' face([2:end 1])'], 2); 0114 planeEdges(pos:sum(Fn(ind2(1:i))), :) = faceEdges; 0115 pos = sum(Fn(ind2(1:i)))+1; 0116 end 0117 planeEdges = unique(planeEdges, 'rows'); 0118 0119 % relabel plane edges 0120 [planeNodes, I, J] = unique(planeEdges(:)); %#ok<ASGLU> 0121 planeEdges2 = reshape(J, size(planeEdges)); 0122 0123 % The set of coplanar faces may not necessarily form a single connected 0124 % component. The following computes label of each connected component. 0125 component = grLabel(nodes(planeNodes, :), planeEdges2); 0126 0127 % compute degree (number of adjacent faces) of each edge. 0128 Npe = size(planeEdges, 1); 0129 edgeDegrees = zeros(Npe, 1); 0130 for i = 1:length(ind2) 0131 face = faces(ind2(i), :); 0132 faceEdges = sort([face' face([2:end 1])'], 2); 0133 for j = 1:size(faceEdges, 1) 0134 indEdge = find(sum(ismember(planeEdges, faceEdges(j,:)),2)==2); 0135 edgeDegrees(indEdge) = edgeDegrees(indEdge)+1; 0136 end 0137 end 0138 0139 % extract unique edges and nodes of the plane 0140 planeEdges = planeEdges(edgeDegrees==1, :); 0141 planeEdges2 = planeEdges2(edgeDegrees==1, :); 0142 0143 % find connected component of each edge 0144 planeEdgesComp = zeros(size(planeEdges, 1), 1); 0145 for iEdge = 1:size(planeEdges, 1) 0146 planeEdgesComp(iEdge) = component(planeEdges2(iEdge, 1)); 0147 end 0148 0149 % iterate on connected faces 0150 for c = 1:max(component) 0151 0152 % convert to chains of nodes 0153 loops = graph2Contours(nodes, planeEdges(planeEdgesComp==c, :)); 0154 0155 % add a simple Polygon for each loop 0156 facePolygon = loops{1}; 0157 for l = 2:length(loops) 0158 facePolygon = [facePolygon, NaN, loops{l}]; %#ok<AGROW> 0159 end 0160 faces2{length(faces2)+1, 1} = facePolygon; 0161 0162 % also add news edges 0163 edges2 = unique([edges2; planeEdges], 'rows'); 0164 end 0165 0166 % mark processed faces 0167 flag(ind2) = 0; 0168 end 0169 0170 0171 %% Additional processing on nodes 0172 0173 % select only nodes which appear in at least one edge 0174 indNodes = unique(edges2(:)); 0175 0176 % for each node, compute index of corresponding new node (or 0 if dropped) 0177 refNodes = zeros(nNodes, 1); 0178 for i = 1:length(indNodes) 0179 refNodes(indNodes(i)) = i; 0180 end 0181 0182 % changes indices of nodes in edges2 array 0183 for i = 1:length(edges2(:)) 0184 edges2(i) = refNodes(edges2(i)); 0185 end 0186 0187 % changes indices of nodes in faces2 array 0188 for iFace = 1:length(faces2) 0189 face = faces2{iFace}; 0190 for i = 1:length(face) 0191 if ~isnan(face(i)) 0192 face(i) = refNodes(face(i)); 0193 end 0194 end 0195 faces2{iFace} = face; 0196 end 0197 0198 % keep only boundary nodes 0199 nodes2 = nodes(indNodes, :); 0200 0201 0202 %% Process output arguments 0203 0204 if nargout == 1 0205 varargout{1} = faces2; 0206 elseif nargout == 2 0207 varargout{1} = nodes2; 0208 varargout{2} = faces2; 0209 elseif nargout == 3 0210 varargout{1} = nodes2; 0211 varargout{2} = edges2; 0212 varargout{3} = faces2; 0213 end 0214 0215 0216 function labels = grLabel(nodes, edges) 0217 %GRLABEL associate a label to each connected component of the graph 0218 % LABELS = grLabel(NODES, EDGES) 0219 % Returns an array with as many rows as the array NODES, containing index 0220 % number of each connected component of the graph. If the graph is 0221 % totally connected, returns an array of 1. 0222 % 0223 % Example 0224 % nodes = rand(6, 2); 0225 % edges = [1 2;1 3;4 6]; 0226 % labels = grLabel(nodes, edges); 0227 % labels = 0228 % 1 0229 % 1 0230 % 1 0231 % 2 0232 % 3 0233 % 2 0234 % 0235 % See also 0236 % getNeighbourNodes 0237 % 0238 % 0239 % ------ 0240 % Author: David Legland 0241 % e-mail: david.legland@grignon.inra.fr 0242 % Created: 2007-08-14, using Matlab 7.4.0.287 (R2007a) 0243 % Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas. 0244 0245 % init 0246 nNodes = size(nodes, 1); 0247 labels = (1:nNodes)'; 0248 0249 % iteration 0250 modif = true; 0251 while modif 0252 modif = false; 0253 0254 for i=1:nNodes 0255 neigh = getNeighbourNodes(i, edges); 0256 neighLabels = labels([i;neigh]); 0257 0258 % check for a modification 0259 if length(unique(neighLabels))>1 0260 modif = true; 0261 end 0262 0263 % put new labels 0264 labels(ismember(labels, neighLabels)) = min(neighLabels); 0265 end 0266 end 0267 0268 % change to have fewer labels 0269 labels2 = unique(labels); 0270 for i = 1:length(labels2) 0271 labels(labels==labels2(i)) = i; 0272 end 0273 0274 function nodes2 = getNeighbourNodes(node, edges) 0275 %GETNEIGHBOURNODES find nodes adjacent to a given node 0276 % 0277 % NEIGHS = getNeighbourNodes(NODE, EDGES) 0278 % NODE: index of the node 0279 % EDGES: the complete edges list 0280 % NEIGHS: the nodes adjacent to the given node. 0281 % 0282 % NODE can also be a vector of node indices, in this case the result is 0283 % the set of neighbors of any input node. 0284 % 0285 % 0286 % ----- 0287 % 0288 % author : David Legland 0289 % INRA - TPV URPOI - BIA IMASTE 0290 % created the 16/08/2004. 0291 % 0292 0293 % HISTORY 0294 % 10/02/2004 documentation 0295 % 13/07/2004 faster algorithm 0296 % 03/10/2007 can specify several input nodes 0297 0298 [i, j] = find(ismember(edges, node)); %#ok<ASGLU> 0299 nodes2 = edges(i,1:2); 0300 nodes2 = unique(nodes2(:)); 0301 nodes2 = sort(nodes2(~ismember(nodes2, node))); 0302 0303 function curves = graph2Contours(nodes, edges) %#ok<INUSL> 0304 %GRAPH2CONTOURS convert a graph to a set of contour curves 0305 % 0306 % CONTOURS = GRAPH2CONTOURS(NODES, EDGES) 0307 % NODES, EDGES is a graph representation (type "help graph" for details) 0308 % The algorithm assume every node has degree 2, and the set of edges 0309 % forms only closed loops. The result is a list of indices arrays, each 0310 % array containing consecutive point indices of a contour. 0311 % 0312 % To transform contours into drawable curves, please use : 0313 % CURVES{i} = NODES(CONTOURS{i}, :); 0314 % 0315 % 0316 % NOTE : contours are not oriented. To manage contour orientation, edges 0317 % also need to be oriented. So we must precise generation of edges. 0318 % 0319 % ----- 0320 % 0321 % author : David Legland 0322 % INRA - TPV URPOI - BIA IMASTE 0323 % created the 05/08/2004. 0324 % 0325 0326 0327 curves = {}; 0328 c = 0; 0329 0330 while size(edges,1)>0 0331 % find first point of the curve 0332 n0 = edges(1,1); 0333 curve = n0; 0334 0335 % second point of the curve 0336 n = edges(1,2); 0337 e = 1; 0338 0339 while true 0340 % add current point to the curve 0341 curve = [curve n]; %#ok<AGROW> 0342 0343 % remove current edge from the list 0344 edges = edges((1:size(edges,1))~=e,:); 0345 0346 % find index of edge containing reference to current node 0347 e = find(edges(:,1)==n | edges(:,2)==n); 0348 e = e(1); 0349 0350 % get index of next current node 0351 % (this is the other node of the current edge) 0352 if edges(e,1)==n 0353 n = edges(e,2); 0354 else 0355 n = edges(e,1); 0356 end 0357 0358 % if node is same as start node, loop is closed, and we stop 0359 % node iteration. 0360 if n==n0 0361 break; 0362 end 0363 end 0364 0365 % remove the last edge of the curve from edge list. 0366 edges = edges((1:size(edges,1))~=e,:); 0367 0368 % add the current curve to the list, and start a new curve 0369 c = c+1; 0370 curves{c} = curve; %#ok<AGROW> 0371 end