Compute the polygons resulting from plane-mesh intersection.
POLYS = intersectPlaneMesh(P, V, F)
Computes the intersection between a plane and a mesh. The result is a
cell array of polygons.
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.
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; drawTorus(torus, 'nTheta', 180, 'nPhi', 180);
hold on; view (3); axis equal; light;
% convert to mesh representation
[v, f] = torusMesh(torus, 'nTheta', 64, 'nPhi', 64);
% 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} = intersectPlaneMesh2(plane, v, f);
end
% draw the resulting 3D polygons
drawPolygon3d(polySet, 'lineWidth', 2, 'color', 'k')
See also
meshes3d, intersectPlanes, intersectEdgePlane0001 function polys = intersectPlaneMesh(plane, v, f) 0002 % Compute the polygons resulting from plane-mesh intersection. 0003 % 0004 % POLYS = intersectPlaneMesh(P, V, F) 0005 % Computes the intersection between a plane and a mesh. The result is a 0006 % cell array of polygons. 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 % 0012 % 0013 % Example 0014 % % Intersect a cube by a plane 0015 % [v f] = createCube; v = v * 10; 0016 % plane = createPlane([5 5 5], [3 4 5]); 0017 % % draw the primitives 0018 % figure; hold on; set(gcf, 'renderer', 'opengl'); 0019 % axis([-10 20 -10 20 -10 20]); view(3); 0020 % drawMesh(v, f); drawPlane3d(plane); 0021 % % compute intersection polygon 0022 % polys = intersectPlaneMesh(plane, v, f); 0023 % drawPolygon3d(polys, 'LineWidth', 2); 0024 % 0025 % % Intersect a torus by a set of planes, and draw the results 0026 % % first creates a torus slightly shifted and rotated 0027 % torus = [.5 .6 .7 30 10 3 4]; 0028 % figure; drawTorus(torus, 'nTheta', 180, 'nPhi', 180); 0029 % hold on; view (3); axis equal; light; 0030 % % convert to mesh representation 0031 % [v, f] = torusMesh(torus, 'nTheta', 64, 'nPhi', 64); 0032 % % compute intersections with collection of planes 0033 % xList = -50:5:50; 0034 % polySet = cell(length(xList), 1); 0035 % for i = 1:length(xList) 0036 % x0 = xList(i); 0037 % plane = createPlane([x0 .5 .5], [1 .2 .3]); 0038 % polySet{i} = intersectPlaneMesh2(plane, v, f); 0039 % end 0040 % % draw the resulting 3D polygons 0041 % drawPolygon3d(polySet, 'lineWidth', 2, 'color', 'k') 0042 % 0043 % 0044 % See also 0045 % meshes3d, intersectPlanes, intersectEdgePlane 0046 % 0047 0048 % ------ 0049 % Author: David Legland 0050 % e-mail: david.legland@inra.fr 0051 % Created: 2012-07-31, using Matlab 7.9.0.529 (R2009b) 0052 % Copyright 2012 INRA - Cepia Software Platform. 0053 0054 e = []; 0055 if isstruct(v) 0056 f = v.faces; 0057 if isfield(v, 'edges') 0058 e = v.edges; 0059 end 0060 v = v.vertices; 0061 end 0062 0063 0064 %% Computation of crossing edges 0065 0066 % compute the edge list 0067 if isempty(e) 0068 e = meshEdges(f); 0069 end 0070 edges = [ v(e(:,1), :) v(e(:,2), :) ]; 0071 0072 % identify which edges cross the mesh 0073 inds = isBelowPlane(v, plane); 0074 edgeCrossInds = find(sum(inds(e), 2) == 1); 0075 0076 % compute one intersection point for each edge 0077 intersectionPoints = intersectEdgePlane(edges(edgeCrossInds, :), plane); 0078 0079 0080 0081 %% mapping edges <-> faces 0082 % identify for each face the indices of edges that intersect the plane, as 0083 % well as for each edge, the indices of the two faces around it. 0084 % We expect each face to contain either 0 or 2 intersecting edges. 0085 % 0086 0087 nFaces = length(f); 0088 faceEdges = cell(1, nFaces); 0089 nCrossEdges = length(edgeCrossInds); 0090 crossEdgeFaces = zeros(nCrossEdges, 2); 0091 0092 for iEdge = 1:length(edgeCrossInds) 0093 edge = e(edgeCrossInds(iEdge), :); 0094 indFaces = find(sum(ismember(f, edge), 2) == 2); 0095 0096 if length(indFaces) ~= 2 0097 error('crossing edge %d (%d,%d) is associated to %d faces', ... 0098 iEdge, edge(1), edge(2), length(indFaces)); 0099 end 0100 0101 crossEdgeFaces(iEdge, :) = indFaces; 0102 0103 for iFace = 1:length(indFaces) 0104 indEdges = faceEdges{indFaces(iFace)}; 0105 indEdges = [indEdges iEdge]; %#ok<AGROW> 0106 faceEdges{indFaces(iFace)} = indEdges; 0107 end 0108 end 0109 0110 0111 %% Iterate on edges and faces to form polygons 0112 0113 % initialize an array indicating which indices need to be processed 0114 nCrossEdges = length(edgeCrossInds); 0115 remainingCrossEdges = true(nCrossEdges, 1); 0116 0117 % create empty cell array of polygons 0118 polys = {}; 0119 0120 % iterate while there are some crossing edges to process 0121 while any(remainingCrossEdges) 0122 0123 % start at any edge, mark it as current 0124 startEdgeIndex = find(remainingCrossEdges, 1, 'first'); 0125 currentEdgeIndex = startEdgeIndex; 0126 0127 % mark current edge as processed 0128 remainingCrossEdges(currentEdgeIndex) = false; 0129 0130 % initialize new set of edge indices 0131 polyEdgeInds = currentEdgeIndex; 0132 0133 % choose one of the two faces around the edge 0134 currentFace = crossEdgeFaces(currentEdgeIndex, 1); 0135 0136 % iterate along current face-edge couples until back to first edge 0137 while true 0138 % find the index of next crossing edge 0139 inds = faceEdges{currentFace}; 0140 currentEdgeIndex = inds(inds ~= currentEdgeIndex); 0141 0142 % mark current edge as processed 0143 remainingCrossEdges(currentEdgeIndex) = false; 0144 0145 % find the index of the other face containing current edge 0146 inds = crossEdgeFaces(currentEdgeIndex, :); 0147 currentFace = inds(inds ~= currentFace); 0148 0149 % check end of current loop 0150 if currentEdgeIndex == startEdgeIndex 0151 break; 0152 end 0153 0154 % add index of current edge 0155 polyEdgeInds = [polyEdgeInds currentEdgeIndex]; %#ok<AGROW> 0156 end 0157 0158 % create polygon, and add it to list of polygons 0159 poly = intersectionPoints(polyEdgeInds, :); 0160 polys = [polys, {poly}]; %#ok<AGROW> 0161 end