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
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))];