1## Copyright (C) 2021 David Legland 2## All rights reserved. 3## 4## Redistribution and use in source and binary forms, with or without 5## modification, are permitted provided that the following conditions are met: 6## 7## 1 Redistributions of source code must retain the above copyright notice, 8## this list of conditions and the following disclaimer. 9## 2 Redistributions in binary form must reproduce the above copyright 10## notice, this list of conditions and the following disclaimer in the 11## documentation and/or other materials provided with the distribution. 12## 13## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' 14## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 15## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 16## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR 17## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 18## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 19## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 20## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 21## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 22## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 23## 24## The views and conclusions contained in the software and documentation are 25## those of the authors and should not be interpreted as representing official 26## policies, either expressed or implied, of the copyright holders. 27 28function newFaces = minConvexHull(points, varargin) 29%MINCONVEXHULL Return the unique minimal convex hull of a set of 3D points. 30% 31% FACES = minConvexHull(PTS) 32% NODES is a set of 3D points (as a Nx3 array). The function computes 33% the convex hull, and merge contiguous coplanar faces. The result is a 34% set of polygonal faces, such that there are no coplanar faces. 35% FACES is a cell array, each cell containing the vector of indices of 36% nodes given in NODES for the corresponding face. 37% 38% FACES = minConvexHull(PTS, PRECISION) 39% Adjust the threshold for deciding if two faces are coplanar or 40% parallel. Default value is 1e-14. 41% 42% Example 43% % extract square faces from a cube 44% [n, e, f] = createCube; 45% f2 = minConvexHull(n); 46% drawMesh(n, f2); 47% 48% % Subdivides and smooths a mesh rpresenting a cube 49% [n, e, f] = createCube; 50% [n2, f2] = subdivideMesh(n, triangulateFaces(f), 4); 51% [n3, f3] = smoothMesh(n2, f2); 52% figure; drawMesh(n3, f3); 53% axis equal; view(3); 54% % merge coplanar faces, making apparent the faces of the original cube 55% f4 = minConvexHull(n3); 56% figure; drawMesh(n3, f4); 57% axis equal; view(3); 58% 59% 60% See also 61% meshes3d, mergeCoplanarFaces, drawMesh, convhull, convhulln 62% 63 64 65% ------ 66% Author: David Legland 67% e-mail: david.legland@inra.fr 68% Created: 2006-07-05 69% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas). 70 71% HISTORY 72% 20/07/2006 add tolerance for coplanarity test 73% 21/08/2006 fix small bug due to difference of methods to test 74% coplanarity, sometimes resulting in 3 points of a face being not 75% coplanar! Also add control on precision 76% 18/09/2007 ensure faces are given as horizontal vectors 77 78% set up precision 79acc = 1e-14; 80if ~isempty(varargin) 81 acc = varargin{1}; 82end 83 84% triangulated convex hull. It is not uniquely defined. 85faces = convhulln(points); 86 87% compute centroid of the nodes 88pointsCentroid = centroid(points); 89 90% number of base triangular faces 91N = size(faces, 1); 92 93% compute normals of given faces 94normals = planeNormal(createPlane(... 95 points(faces(:,1),:), points(faces(:,2),:), points(faces(:,3),:))); 96 97% initialize empty faces 98newFaces = {}; 99 100 101% Processing flag for each triangle 102% 1 : triangle to process, 0 : already processed 103% in the beginning, every triangle face need to be processed 104flag = ones(N, 1); 105 106% iterate on each triangular face of the convex hull 107for iFace = 1:N 108 109 % check if face was already performed 110 if ~flag(iFace) 111 continue; 112 end 113 114 % indices of faces with same normal 115 ind = find(abs(vectorNorm3d(cross(repmat(normals(iFace, :), [N 1]), normals)))<acc); 116 ind = ind(ind~=iFace); 117 118 % keep only coplanar faces (test coplanarity of points in both face) 119 ind2 = iFace; 120 for j = 1:length(ind) 121 if isCoplanar(points([faces(iFace,:) faces(ind(j),:)], :), acc) 122 ind2 = [ind2 ind(j)]; %#ok<AGROW> 123 end 124 end 125 126 127 % compute order of the vertices in current face 128 faceVertices = unique(faces(ind2, :)); 129 [tmp, I] = angleSort3d(points(faceVertices, :)); %#ok<ASGLU> 130 131 % create the new face, ensuring it is a row vector 132 face = faceVertices(I); 133 face = face(:)'; 134 135 % ensure face has normal pointing outwards 136 outerNormal = meshFaceCentroids(points, face) - pointsCentroid; 137 if dot(meshFaceNormals(points, face), outerNormal, 2) < 0 138 face = face([1 end:-1:2]); 139 end 140 141 % add a new face to the list 142 newFaces = [newFaces {face}]; %#ok<AGROW> 143 144 % mark processed faces 145 flag(ind2) = 0; 146end 147 148