meshBoundary

MESHBOUNDARY Boundary of a mesh as a collection of 3D line strings.

   CURVES = meshBoundary(V, F)

   Example
     % Example 1
     % create centered icosahedron
     [v, f] = createIcosahedron;
     v(:,3) = v(:,3) - mean(v(:,3));
     % convert to simili-sphere
     [v2, f2] = subdivideMesh(v, f, 3);
     v3 = normalizeVector3d(v2);
     % clip with plane
     plane = createPlane([0 0 0], [-1 -2 3]);
     [vc, fc] = clipMeshByPlane(v3, f2, plane);
     figure; drawMesh(vc, fc); axis equal; view([80 -10]);
     % draw mesh boundary
     curves = meshBoundary(vc, fc);
     hold on; drawPolygon3d(curves, 'linewidth', 2, 'color', 'b');

     % Example 2
     mesh = readMesh('mushroom.off');
     plane = createPlane([0 0 0.7], [1 2 -3]);
     [vc, fc] = clipMeshByPlane(mesh, plane);
     bnd = meshBoundary(vc, fc);
     figure; drawMesh(vc, fc); axis equal; view([-40 50]);
     drawPolygon3d(bnd, 'linewidth', 2, 'color', 'b');

   See also
     meshes3d, meshBoundaryEdgeIndices, meshBoundaryVertexIndices,
     clipMeshByPlane, clipMeshVertices

0001 function polyList = meshBoundary(varargin)
0002 %MESHBOUNDARY Boundary of a mesh as a collection of 3D line strings.
0003 %
0004 %   CURVES = meshBoundary(V, F)
0005 %
0006 %   Example
0007 %     % Example 1
0008 %     % create centered icosahedron
0009 %     [v, f] = createIcosahedron;
0010 %     v(:,3) = v(:,3) - mean(v(:,3));
0011 %     % convert to simili-sphere
0012 %     [v2, f2] = subdivideMesh(v, f, 3);
0013 %     v3 = normalizeVector3d(v2);
0014 %     % clip with plane
0015 %     plane = createPlane([0 0 0], [-1 -2 3]);
0016 %     [vc, fc] = clipMeshByPlane(v3, f2, plane);
0017 %     figure; drawMesh(vc, fc); axis equal; view([80 -10]);
0018 %     % draw mesh boundary
0019 %     curves = meshBoundary(vc, fc);
0020 %     hold on; drawPolygon3d(curves, 'linewidth', 2, 'color', 'b');
0021 %
0022 %     % Example 2
0023 %     mesh = readMesh('mushroom.off');
0024 %     plane = createPlane([0 0 0.7], [1 2 -3]);
0025 %     [vc, fc] = clipMeshByPlane(mesh, plane);
0026 %     bnd = meshBoundary(vc, fc);
0027 %     figure; drawMesh(vc, fc); axis equal; view([-40 50]);
0028 %     drawPolygon3d(bnd, 'linewidth', 2, 'color', 'b');
0029 %
0030 %   See also
0031 %     meshes3d, meshBoundaryEdgeIndices, meshBoundaryVertexIndices,
0032 %     clipMeshByPlane, clipMeshVertices
0033 %
0034 
0035 % ------
0036 % Author: David Legland
0037 % E-mail: david.legland@inrae.fr
0038 % Created: 2019-05-01, using Matlab 8.6.0.267246 (R2015b)
0039 % Copyright 2019-2024 INRA - Cepia Software Platform
0040 
0041 [vertices, edges, faces] = parseMeshData(varargin{:});
0042 
0043 % Compute edge-vertex map if not specified
0044 if isempty(edges)
0045     edges = meshEdges(vertices, faces);
0046 end
0047 
0048 bndEdgeInds = meshBoundaryEdgeIndices(vertices, edges, faces);
0049 bndEdges = edges(bndEdgeInds, :);
0050 
0051 if isempty(bndEdgeInds)
0052     polyList = {};
0053     return;
0054 end
0055 
0056 % allocate empty array
0057 polyList = {};
0058 
0059 nPolys = 0;
0060 while ~isempty(bndEdges)
0061     nPolys = nPolys + 1;
0062     
0063     % current edge
0064     edge = bndEdges(1, :);
0065 
0066     % initialize new polyline at first vertex
0067     ind0 = edge(1);
0068     vertexInds = ind0;
0069     
0070     % current vertex
0071     index = edge(2);
0072     bndEdges(1, :) = [];
0073     
0074     % iterate over edges until current vertex becomes ind0
0075     while index ~= ind0
0076         % append current vertex to list of indices for current polygon
0077         vertexInds = [vertexInds ; index]; %#ok<AGROW>
0078         
0079         % index of the next edge containing current vertex
0080         edgeInd = find(sum(bndEdges == index, 2) > 0);
0081         
0082         % check validity
0083         if isempty(edgeInd)
0084             error('could not find next edge for vertex index %d', index);
0085         end
0086         if length(edgeInd) > 1
0087             error('two many edges contains vertex index %d', index);
0088         end
0089         
0090         % remove current edge from the list of edges to process
0091         edge = bndEdges(edgeInd, :);
0092         bndEdges(edgeInd, :) = [];
0093         
0094 %         % check if current edge closes current polygon
0095 %         if index == ind0
0096 %             break;
0097 %         end
0098         
0099         % identify the next index
0100         index = edge(edge ~= index);
0101     end
0102     
0103     % create the 3D polyline
0104     polyList{nPolys} = vertices(vertexInds, :); %#ok<AGROW>
0105 end