polyhedronSlice

POLYHEDRONSLICE Intersect a convex polyhedron with a plane.

   SLICE = polyhedronSlice(NODES, FACES, PLANE)
   NODES: a Nx3 array
   FACES: either a cell array or a Nf*3 or Nf*4 array
   PLANE: a plane representation [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2].
   return the intersection polygon of the polyhedra with the plane, in the
   form of a set of ordered points.

   Works only for convex polyhedra.

   Example
   polyhedronSlice

   See also
   polyhedra, clipConvexPolyhedronHP, intersectPlaneMesh

0001 function points = polyhedronSlice(nodes, faces, plane)
0002 %POLYHEDRONSLICE Intersect a convex polyhedron with a plane.
0003 %
0004 %   SLICE = polyhedronSlice(NODES, FACES, PLANE)
0005 %   NODES: a Nx3 array
0006 %   FACES: either a cell array or a Nf*3 or Nf*4 array
0007 %   PLANE: a plane representation [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2].
0008 %   return the intersection polygon of the polyhedra with the plane, in the
0009 %   form of a set of ordered points.
0010 %
0011 %   Works only for convex polyhedra.
0012 %
0013 %   Example
0014 %   polyhedronSlice
0015 %
0016 %   See also
0017 %   polyhedra, clipConvexPolyhedronHP, intersectPlaneMesh
0018 %
0019 
0020 % ------
0021 % Author: David Legland
0022 % e-mail: david.legland@nantes.inra.fr
0023 % Created: 2007-09-18,    using Matlab 7.4.0.287 (R2007a)
0024 % Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.
0025 
0026 % if faces is a numeric array, convert it to cell array
0027 if isnumeric(faces)
0028     faces2 = cell(size(faces, 1), 1);
0029     for f = 1:length(faces2)
0030         faces2{f} = faces(f,:);
0031     end
0032     faces = faces2;
0033 else
0034     % ensure we have face with horizontal vectors...
0035     for f = 1:length(faces)
0036         face = faces{f};
0037         faces{f} = face(:)';
0038     end
0039 end
0040 
0041 % compute edges of the polyhedron
0042 inds = zeros(0, 2);
0043 for f = 1:length(faces)
0044     face = faces{f}';
0045     inds = [inds ; sort([face face([2:end 1])], 2)]; %#ok<AGROW>
0046 end
0047 inds = unique(inds, 'rows');
0048 edges = [nodes(inds(:,1), :) nodes(inds(:,2), :)];
0049 
0050 % intersection of edges with plane
0051 points = intersectEdgePlane(edges, plane);
0052 points = points(sum(isfinite(points), 2)==3, :);
0053 
0054 if ~isempty(points)
0055     points = angleSort3d(points);
0056 end