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.

   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 %   FACES = minConvexHull(PTS, PRECISION)
0012 %   Adjust the threshold for deciding if two faces are coplanar or
0013 %   parallel. Default value is 1e-14.
0014 %
0015 %   Example
0016 %     % extract square faces from a cube
0017 %     [n, e, f] = createCube;
0018 %     f2 = minConvexHull(n);
0019 %     drawMesh(n, f2);
0020 %
0021 %     % Subdivides and smooths a mesh rpresenting a cube
0022 %     [n, e, f] = createCube;
0023 %     [n2, f2] = subdivideMesh(n, triangulateFaces(f), 4);
0024 %     [n3, f3] = smoothMesh(n2, f2);
0025 %     figure; drawMesh(n3, f3);
0026 %     axis equal; view(3);
0027 %     % merge coplanar faces, making apparent the faces of the original cube
0028 %     f4 = minConvexHull(n3);
0029 %     figure; drawMesh(n3, f4);
0030 %     axis equal; view(3);
0031 %
0032 %
0033 %   See also
0034 %   meshes3d, mergeCoplanarFaces, drawMesh, convhull, convhulln
0035 %
0036 
0037 
0038 % ------
0039 % Author: David Legland
0040 % e-mail: david.legland@inra.fr
0041 % Created: 2006-07-05
0042 % Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).
0043 
0044 % HISTORY
0045 %   20/07/2006 add tolerance for coplanarity test
0046 %   21/08/2006 fix small bug due to difference of methods to test
0047 %       coplanarity, sometimes resulting in 3 points of a face being not
0048 %       coplanar! Also add control on precision
0049 %   18/09/2007 ensure faces are given as horizontal vectors
0050 
0051 % set up precision
0052 acc = 1e-14;
0053 if ~isempty(varargin)
0054     acc = varargin{1};
0055 end
0056 
0057 % triangulated convex hull. It is not uniquely defined.
0058 faces = convhulln(points);
0059 
0060 % compute centroid of the nodes
0061 pointsCentroid = centroid(points);
0062 
0063 % number of base triangular faces
0064 N = size(faces, 1);
0065 
0066 % compute normals of given faces
0067 normals = planeNormal(createPlane(...
0068     points(faces(:,1),:), points(faces(:,2),:), points(faces(:,3),:)));
0069 
0070 % initialize empty faces
0071 newFaces = {};
0072 
0073 
0074 % Processing flag for each triangle
0075 % 1 : triangle to process, 0 : already processed
0076 % in the beginning, every triangle face need to be processed
0077 flag = ones(N, 1);
0078 
0079 % iterate on each triangular face of the convex hull
0080 for iFace = 1:N
0081     
0082     % check if face was already performed
0083     if ~flag(iFace)
0084         continue;
0085     end
0086 
0087     % indices of faces with same normal
0088     ind = find(abs(vectorNorm3d(cross(repmat(normals(iFace, :), [N 1]), normals)))<acc);
0089     ind = ind(ind~=iFace);
0090     
0091     % keep only coplanar faces (test coplanarity of points in both face)
0092     ind2 = iFace;
0093     for j = 1:length(ind)
0094         if isCoplanar(points([faces(iFace,:) faces(ind(j),:)], :), acc)
0095             ind2 = [ind2 ind(j)]; %#ok<AGROW>
0096         end
0097     end
0098     
0099     
0100     % compute order of the vertices in current face
0101     faceVertices = unique(faces(ind2, :));
0102     [tmp, I]  = angleSort3d(points(faceVertices, :)); %#ok<ASGLU>
0103     
0104     % create the new face, ensuring it is a row vector
0105     face = faceVertices(I);
0106     face = face(:)';
0107     
0108     % ensure face has normal pointing outwards
0109     outerNormal = meshFaceCentroids(points, face) - pointsCentroid;
0110     if dot(meshFaceNormals(points, face), outerNormal, 2) < 0
0111         face = face([1 end:-1:2]);
0112     end
0113     
0114     % add a new face to the list
0115     newFaces = [newFaces {face}]; %#ok<AGROW>
0116     
0117     % mark processed faces
0118     flag(ind2) = 0;
0119 end
0120