minConvexHull

MINCONVEXHULL Return the unique minimal convex hull of a set of 3D points.

   FACES = minConvexHull(PTS)
   NODES is a set of 3D points  (as a Nx3 array). The function computes
   the convex hull, and merge contiguous coplanar faces. The result is a
   set of polygonal faces, such that there are no coplanar faces.
   FACES is a cell array, each cell containing the vector of indices of
   nodes given in NODES for the corresponding face.

   NOTE: minConvexHull does not remove unreferenced vertices, etc. The
   function trimMesh can be subsequently applied to the mesh to adress
   this.

   FACES = minConvexHull(PTS, PRECISION)
   Adjust the threshold for deciding if two faces are coplanar or
   parallel. Default value is 1e-14.

   Example
     % extract square faces from a cube
     [n, e, f] = createCube;
     f2 = minConvexHull(n);
     drawMesh(n, f2);

     % Subdivides and smooths a mesh rpresenting a cube
     [n, e, f] = createCube;
     [n2, f2] = subdivideMesh(n, triangulateFaces(f), 4);
     [n3, f3] = smoothMesh(n2, f2);
     figure; drawMesh(n3, f3);
     axis equal; view(3);
     % merge coplanar faces, making apparent the faces of the original cube
     f4 = minConvexHull(n3);
     figure; drawMesh(n3, f4);
     axis equal; view(3);


   See also
   meshes3d, mergeCoplanarFaces, drawMesh, convhull, convhulln

0001 function newFaces = minConvexHull(points, varargin)
0002 %MINCONVEXHULL Return the unique minimal convex hull of a set of 3D points.
0003 %
0004 %   FACES = minConvexHull(PTS)
0005 %   NODES is a set of 3D points  (as a Nx3 array). The function computes
0006 %   the convex hull, and merge contiguous coplanar faces. The result is a
0007 %   set of polygonal faces, such that there are no coplanar faces.
0008 %   FACES is a cell array, each cell containing the vector of indices of
0009 %   nodes given in NODES for the corresponding face.
0010 %
0011 %   NOTE: minConvexHull does not remove unreferenced vertices, etc. The
0012 %   function trimMesh can be subsequently applied to the mesh to adress
0013 %   this.
0014 %
0015 %   FACES = minConvexHull(PTS, PRECISION)
0016 %   Adjust the threshold for deciding if two faces are coplanar or
0017 %   parallel. Default value is 1e-14.
0018 %
0019 %   Example
0020 %     % extract square faces from a cube
0021 %     [n, e, f] = createCube;
0022 %     f2 = minConvexHull(n);
0023 %     drawMesh(n, f2);
0024 %
0025 %     % Subdivides and smooths a mesh rpresenting a cube
0026 %     [n, e, f] = createCube;
0027 %     [n2, f2] = subdivideMesh(n, triangulateFaces(f), 4);
0028 %     [n3, f3] = smoothMesh(n2, f2);
0029 %     figure; drawMesh(n3, f3);
0030 %     axis equal; view(3);
0031 %     % merge coplanar faces, making apparent the faces of the original cube
0032 %     f4 = minConvexHull(n3);
0033 %     figure; drawMesh(n3, f4);
0034 %     axis equal; view(3);
0035 %
0036 %
0037 %   See also
0038 %   meshes3d, mergeCoplanarFaces, drawMesh, convhull, convhulln
0039 %
0040 
0041 % ------
0042 % Author: David Legland
0043 % E-mail: david.legland@inrae.fr
0044 % Created: 2006-07-05
0045 % Copyright 2006-2024 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)
0046 
0047 % set up precision
0048 acc = 1e-14;
0049 if ~isempty(varargin)
0050     acc = varargin{1};
0051 end
0052 
0053 % triangulated convex hull. It is not uniquely defined.
0054 faces = convhulln(points);
0055 
0056 % compute centroid of the nodes
0057 pointsCentroid = centroid(points);
0058 
0059 % number of base triangular faces
0060 N = size(faces, 1);
0061 
0062 % compute normals of given faces
0063 normals = planeNormal(createPlane(...
0064     points(faces(:,1),:), points(faces(:,2),:), points(faces(:,3),:)));
0065 
0066 % initialize empty faces
0067 newFaces = {};
0068 
0069 
0070 % Processing flag for each triangle
0071 % 1 : triangle to process, 0 : already processed
0072 % in the beginning, every triangle face need to be processed
0073 flag = ones(N, 1);
0074 
0075 % iterate on each triangular face of the convex hull
0076 for iFace = 1:N
0077     
0078     % check if face was already performed
0079     if ~flag(iFace)
0080         continue;
0081     end
0082 
0083     % indices of faces with same normal
0084     ind = find(abs(vectorNorm3d(cross(repmat(normals(iFace, :), [N 1]), normals)))<acc);
0085     ind = ind(ind~=iFace);
0086     
0087     % keep only coplanar faces (test coplanarity of points in both face)
0088     ind2 = iFace;
0089     for j = 1:length(ind)
0090         if isCoplanar(points([faces(iFace,:) faces(ind(j),:)], :), acc)
0091             ind2 = [ind2 ind(j)]; %#ok<AGROW>
0092         end
0093     end
0094     
0095     
0096     % compute order of the vertices in current face
0097     faceVertices = unique(faces(ind2, :));
0098     [tmp, I]  = angleSort3d(points(faceVertices, :)); %#ok<ASGLU>
0099     
0100     % create the new face, ensuring it is a row vector
0101     face = faceVertices(I);
0102     face = face(:)';
0103     
0104     % ensure face has normal pointing outwards
0105     outerNormal = meshFaceCentroids(points, face) - pointsCentroid;
0106     if dot(meshFaceNormals(points, face), outerNormal, 2) < 0
0107         face = face([1 end:-1:2]);
0108     end
0109     
0110     % add a new face to the list
0111     newFaces = [newFaces {face}]; %#ok<AGROW>
0112     
0113     % mark processed faces
0114     flag(ind2) = 0;
0115 end