Home > matGeom > meshes3d > intersectPlaneMesh.m

intersectPlaneMesh

PURPOSE ^

INTERSECTPLANEMESH Compute the polylines resulting from plane-mesh intersection.

SYNOPSIS ^

function [polys, closedFlag] = intersectPlaneMesh(plane, v, f)

DESCRIPTION ^

INTERSECTPLANEMESH Compute the polylines resulting from plane-mesh intersection.

   POLYS = intersectPlaneMesh(P, V, F)
   [POLYS, CLOSED] = intersectPlaneMesh(P, V, F)
   Computes the intersection between a plane and a mesh. 
   The plane P is given as:
   P = [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2]
   The mesh is given as numeric array V of vertex coordinates and an array
   of (triangular) face vertex indices.
   The output POLYS is a cell array of polylines, where each cell contains
   a NV-by-3 numeric array of coordinates. The (optional) output CLOSED is
   a logical array the same size as the POLYS, indicating whether the
   corresponding polylines are closed (true), or open (false). 
   Use the functions 'drawPolygon3d' to display closed polylines, and
   'drawPolyline3d' to display open polylines.


   Example
     % Intersect a cube by a plane
     [v, f] = createCube; v = v * 10;
     plane = createPlane([5 5 5], [3 4 5]);
     % draw the primitives
     figure; hold on; set(gcf, 'renderer', 'opengl');
     axis([-10 20 -10 20 -10 20]); view(3);
     drawMesh(v, f); drawPlane3d(plane);
     % compute intersection polygon
     polys = intersectPlaneMesh(plane, v, f);
     drawPolygon3d(polys, 'LineWidth', 2);

     % Intersect a torus by a set of planes, and draw the results
     % first creates a torus slightly shifted and rotated
     torus = [.5 .6 .7   30 10   3 4];
     figure('color','w');
     % convert to mesh representation
     [v, f] = torusMesh(torus, 'nTheta', 64, 'nPhi', 64);
     f = triangulateFaces(f);
     drawMesh(v, f);
     hold on; view (3); axis equal; light;
     % compute intersections with collection of planes
     xList = -50:5:50;
     polySet = cell(length(xList), 1);
     for i = 1:length(xList)
         x0 = xList(i);
         plane = createPlane([x0 .5 .5], [1 .2 .3]);
         polySet{i} = intersectPlaneMesh(plane, v, f);
     end
     % draw the resulting 3D polygons
     drawPolygon3d(polySet, 'lineWidth', 2, 'color', 'y')

     % Demonstrate ability to draw open mesh intersections
     poly = circleArcToPolyline([10 0 5 90 180], 33);
     [x, y, z] = revolutionSurface(poly, linspace(-pi, pi, 65));
     [v, f] = surfToMesh(x, y, z);
     f = triangulateFaces(f);
     plane = createPlane([0 0 0], [5 2 -4]);
     figure; hold on; axis equal; view(3);
     drawMesh(v, f, 'linestyle', 'none', 'facecolor', [0.0 0.8 0.0], 'faceAlpha', 0.7);
     drawPlane3d(plane, 'facecolor', 'm', 'faceAlpha', 0.5);
     % compute and display intersection
     [curves, closed] = intersectPlaneMesh(plane, v, f);
     drawPolyline3d(curves(~closed), 'linewidth', 2, 'color', 'b')


   See also
     meshes3d, intersectPlanes, intersectEdgePlane

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [polys, closedFlag] = intersectPlaneMesh(plane, v, f)
0002 %INTERSECTPLANEMESH Compute the polylines resulting from plane-mesh intersection.
0003 %
0004 %   POLYS = intersectPlaneMesh(P, V, F)
0005 %   [POLYS, CLOSED] = intersectPlaneMesh(P, V, F)
0006 %   Computes the intersection between a plane and a mesh.
0007 %   The plane P is given as:
0008 %   P = [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2]
0009 %   The mesh is given as numeric array V of vertex coordinates and an array
0010 %   of (triangular) face vertex indices.
0011 %   The output POLYS is a cell array of polylines, where each cell contains
0012 %   a NV-by-3 numeric array of coordinates. The (optional) output CLOSED is
0013 %   a logical array the same size as the POLYS, indicating whether the
0014 %   corresponding polylines are closed (true), or open (false).
0015 %   Use the functions 'drawPolygon3d' to display closed polylines, and
0016 %   'drawPolyline3d' to display open polylines.
0017 %
0018 %
0019 %   Example
0020 %     % Intersect a cube by a plane
0021 %     [v, f] = createCube; v = v * 10;
0022 %     plane = createPlane([5 5 5], [3 4 5]);
0023 %     % draw the primitives
0024 %     figure; hold on; set(gcf, 'renderer', 'opengl');
0025 %     axis([-10 20 -10 20 -10 20]); view(3);
0026 %     drawMesh(v, f); drawPlane3d(plane);
0027 %     % compute intersection polygon
0028 %     polys = intersectPlaneMesh(plane, v, f);
0029 %     drawPolygon3d(polys, 'LineWidth', 2);
0030 %
0031 %     % Intersect a torus by a set of planes, and draw the results
0032 %     % first creates a torus slightly shifted and rotated
0033 %     torus = [.5 .6 .7   30 10   3 4];
0034 %     figure('color','w');
0035 %     % convert to mesh representation
0036 %     [v, f] = torusMesh(torus, 'nTheta', 64, 'nPhi', 64);
0037 %     f = triangulateFaces(f);
0038 %     drawMesh(v, f);
0039 %     hold on; view (3); axis equal; light;
0040 %     % compute intersections with collection of planes
0041 %     xList = -50:5:50;
0042 %     polySet = cell(length(xList), 1);
0043 %     for i = 1:length(xList)
0044 %         x0 = xList(i);
0045 %         plane = createPlane([x0 .5 .5], [1 .2 .3]);
0046 %         polySet{i} = intersectPlaneMesh(plane, v, f);
0047 %     end
0048 %     % draw the resulting 3D polygons
0049 %     drawPolygon3d(polySet, 'lineWidth', 2, 'color', 'y')
0050 %
0051 %     % Demonstrate ability to draw open mesh intersections
0052 %     poly = circleArcToPolyline([10 0 5 90 180], 33);
0053 %     [x, y, z] = revolutionSurface(poly, linspace(-pi, pi, 65));
0054 %     [v, f] = surfToMesh(x, y, z);
0055 %     f = triangulateFaces(f);
0056 %     plane = createPlane([0 0 0], [5 2 -4]);
0057 %     figure; hold on; axis equal; view(3);
0058 %     drawMesh(v, f, 'linestyle', 'none', 'facecolor', [0.0 0.8 0.0], 'faceAlpha', 0.7);
0059 %     drawPlane3d(plane, 'facecolor', 'm', 'faceAlpha', 0.5);
0060 %     % compute and display intersection
0061 %     [curves, closed] = intersectPlaneMesh(plane, v, f);
0062 %     drawPolyline3d(curves(~closed), 'linewidth', 2, 'color', 'b')
0063 %
0064 %
0065 %   See also
0066 %     meshes3d, intersectPlanes, intersectEdgePlane
0067 %
0068 
0069 % ------
0070 % Author: David Legland
0071 % E-mail: david.legland@inrae.fr
0072 % Created: 2012-07-31, using Matlab 7.9.0.529 (R2009b)
0073 % Copyright 2012-2024 INRA - Cepia Software Platform
0074 
0075 e = [];
0076 if isstruct(v)
0077     f = v.faces;
0078     if isfield(v, 'edges')
0079         e = v.edges;
0080     end
0081     v = v.vertices;
0082 end
0083 
0084 
0085 %% Computation of crossing edges
0086 
0087 % compute the edge list
0088 if isempty(e)
0089     e = meshEdges(f);
0090 end
0091 edges = [ v(e(:,1), :) v(e(:,2), :) ];
0092 
0093 % identify which edges cross the mesh
0094 inds = isBelowPlane(v, plane);
0095 edgeCrossInds = find(sum(inds(e), 2) == 1);
0096 
0097 % compute one intersection point for each edge
0098 intersectionPoints = intersectEdgePlane(edges(edgeCrossInds, :), plane);
0099 
0100 
0101 
0102 %% mapping edges <-> faces
0103 % identify for each face the indices of edges that intersect the plane, as
0104 % well as for each edge, the indices of the two faces around it.
0105 % We expect each face to contain either 0 or 2 intersecting edges.
0106 %
0107 
0108 nFaces = length(f);
0109 faceEdges = cell(1, nFaces);
0110 nCrossEdges = length(edgeCrossInds);
0111 crossEdgeFaces = cell(nCrossEdges, 1);
0112 
0113 for iEdge = 1:length(edgeCrossInds)
0114     % identify index of faces adjacent to edge
0115     edge = e(edgeCrossInds(iEdge), :);
0116     indFaces = find(sum(ismember(f, edge), 2) == 2);
0117     
0118     % assert mesh is manifold (no edge connected to more than three faces)
0119     if length(indFaces) > 2
0120         error('crossing edge %d (%d,%d) is associated to %d faces', ...
0121             iEdge, edge(1), edge(2), length(indFaces));
0122     end
0123     
0124     crossEdgeFaces{iEdge} = indFaces;
0125     
0126     % add current edge to list of edges associated to each face
0127     for iFace = 1:length(indFaces)
0128         indEdges = faceEdges{indFaces(iFace)};
0129         indEdges = [indEdges iEdge]; %#ok<AGROW>
0130         faceEdges{indFaces(iFace)} = indEdges;
0131     end
0132 end
0133 
0134 % initialize an array indicating which edges need to be processed
0135 nCrossEdges = length(edgeCrossInds);
0136 remainingCrossEdges = true(nCrossEdges, 1);
0137 
0138 
0139 %% Iterate on edges and faces to form open poylines
0140 
0141 % create empty cell array of open polylines
0142 openPolys = {};
0143 
0144 % identify crossing edges at extremity of open polylines
0145 extremityEdgeInds = find(cellfun(@length, crossEdgeFaces) == 1);
0146 remainingExtremities = true(length(extremityEdgeInds), 1);
0147 
0148 % iterate while there are remaining extremity crossing edges
0149 while any(remainingExtremities)
0150     % start from arbitrary remaining extremity
0151     extremityIndex = find(remainingExtremities, 1, 'first');
0152     remainingExtremities(extremityIndex) = false;
0153 
0154     % use extremity as current edge
0155     startEdgeIndex = extremityEdgeInds(extremityIndex);
0156     currentEdgeIndex = startEdgeIndex;
0157     
0158     % mark current edge as processed
0159     remainingCrossEdges(currentEdgeIndex) = false;
0160     
0161     % initialize new set of edge indices
0162     polyEdgeInds = currentEdgeIndex;
0163 
0164     % find the unique face adjacent to current edge
0165     edgeFaces = crossEdgeFaces{currentEdgeIndex};
0166     currentFace = edgeFaces(1);
0167 
0168     % iterate along current face-edge couples until back to first edge
0169     while true
0170         % find the index of next crossing edge
0171         inds = faceEdges{currentFace};
0172         currentEdgeIndex = inds(inds ~= currentEdgeIndex);
0173         
0174         % add index of current edge
0175         polyEdgeInds = [polyEdgeInds currentEdgeIndex]; %#ok<AGROW>
0176 
0177         % mark current edge as processed
0178         remainingCrossEdges(currentEdgeIndex) = false;
0179     
0180         % find the index of the other face containing current edge
0181         inds = crossEdgeFaces{currentEdgeIndex};
0182 
0183         % check if we found an extremity edge
0184         if isscalar(inds)
0185             ind = extremityEdgeInds == currentEdgeIndex;
0186             remainingExtremities(ind) = false;
0187             break;
0188         end
0189 
0190         % switch to next face
0191         currentFace = inds(inds ~= currentFace);
0192     end
0193     
0194     % create polygon, and add it to list of polygons
0195     poly = intersectionPoints(polyEdgeInds, :);
0196     openPolys = [openPolys, {poly}]; %#ok<AGROW>
0197 end
0198 
0199 
0200 %% Iterate on edges and faces to form closed polylines
0201 
0202 % create empty cell array of polygons
0203 rings = {};
0204 
0205 % iterate while there are some crossing edges to process
0206 while any(remainingCrossEdges)
0207     
0208     % start at any edge, mark it as current
0209     startEdgeIndex = find(remainingCrossEdges, 1, 'first');
0210     currentEdgeIndex = startEdgeIndex;
0211     
0212     % mark current edge as processed
0213     remainingCrossEdges(currentEdgeIndex) = false;
0214     
0215     % initialize new set of edge indices
0216     polyEdgeInds = currentEdgeIndex;
0217 
0218     % choose one of the two faces around the edge
0219     edgeFaces = crossEdgeFaces{currentEdgeIndex};
0220     currentFace = edgeFaces(1);
0221 
0222     % iterate along current face-edge couples until back to first edge
0223     while true
0224         % find the index of next crossing edge
0225         inds = faceEdges{currentFace};
0226         currentEdgeIndex = inds(inds ~= currentEdgeIndex);
0227      
0228         % mark current edge as processed
0229         remainingCrossEdges(currentEdgeIndex) = false;
0230     
0231         % check end of current loop
0232         if currentEdgeIndex == startEdgeIndex
0233             break;
0234         end
0235         
0236         % add index of current edge
0237         polyEdgeInds = [polyEdgeInds currentEdgeIndex]; %#ok<AGROW>
0238 
0239         % find the index of the other face containing current edge
0240         inds = crossEdgeFaces{currentEdgeIndex};
0241         currentFace = inds(inds ~= currentFace);
0242     end
0243     
0244     % create polygon, and add it to list of polygons
0245     poly = intersectionPoints(polyEdgeInds, :);
0246     rings = [rings, {poly}]; %#ok<AGROW>
0247 end
0248 
0249 
0250 %% Format output array
0251 polys = [rings, openPolys];
0252 closedFlag = [true(1, length(rings)), false(1, length(openPolys))];

Generated on Thu 21-Nov-2024 11:30:22 by m2html © 2003-2022