clipPolygon3dHP

CLIPPOLYGON3DHP clip a 3D polygon with Half-space.

   usage
   POLY2 = clipPolygon3dHP(POLY, PLANE)
   POLY is a [Nx3] array of points, and PLANE is given as :
   [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2].
   The result POLY2 is also an array of 3d points, sometimes smaller than
   poly, and that can be [0x3] (empty polygon).

   POLY2 = clipPolygon3dHP(POLY, PT0, NORMAL)
   uses plane with normal NORMAL and containing point PT0.

   TODO: not yet implemented

   There is a problem for non-convex polygons, as they can be clipped in
   several polygons. Possible solutions:
   * create another function 'clipConvexPolygon3dPlane' or
       'clipConvexPolygon3d', using a simplified algorithm
   * returns a list of polygons instead of a single polygon,
   * in the case of one polygon as return decide what to return
   * and rename this function to 'clipPolygon3d'

   See also:
   poygons3d, polyhedra, clipConvexPolygon3dHP

   ---------

   author : David Legland
   INRA - TPV URPOI - BIA IMASTE
   created the 02/08/2005.

0001 function poly2 = clipPolygon3dHP(poly, plane)
0002 %CLIPPOLYGON3DHP clip a 3D polygon with Half-space.
0003 %
0004 %   usage
0005 %   POLY2 = clipPolygon3dHP(POLY, PLANE)
0006 %   POLY is a [Nx3] array of points, and PLANE is given as :
0007 %   [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2].
0008 %   The result POLY2 is also an array of 3d points, sometimes smaller than
0009 %   poly, and that can be [0x3] (empty polygon).
0010 %
0011 %   POLY2 = clipPolygon3dHP(POLY, PT0, NORMAL)
0012 %   uses plane with normal NORMAL and containing point PT0.
0013 %
0014 %   TODO: not yet implemented
0015 %
0016 %   There is a problem for non-convex polygons, as they can be clipped in
0017 %   several polygons. Possible solutions:
0018 %   * create another function 'clipConvexPolygon3dPlane' or
0019 %       'clipConvexPolygon3d', using a simplified algorithm
0020 %   * returns a list of polygons instead of a single polygon,
0021 %   * in the case of one polygon as return decide what to return
0022 %   * and rename this function to 'clipPolygon3d'
0023 %
0024 %   See also:
0025 %   poygons3d, polyhedra, clipConvexPolygon3dHP
0026 %
0027 %   ---------
0028 %
0029 %   author : David Legland
0030 %   INRA - TPV URPOI - BIA IMASTE
0031 %   created the 02/08/2005.
0032 %
0033 
0034 %   HISTORY
0035 %   2007-01-04 add todo flag
0036 %   2011-08-17 rewrite algo, that works for convex polygons, but is slower
0037 %       than function clipConvexPolgon3dHP
0038 
0039 
0040 %% Pre-Processing
0041 
0042 % ensure last point is the same as the first one (makes computation easier)
0043 if sum(poly(end, :) == poly(1,:)) ~= 3
0044     poly = [poly; poly(1,:)];
0045 end
0046 
0047 % compute index of position wrt plane for each vertex
0048 below = isBelowPlane(poly, plane);
0049 
0050 % in the case of a polygon totally over the plane, return empty array
0051 if sum(below) == 0
0052     poly2 = zeros(0, 3);
0053     return;
0054 end
0055 
0056 % in the case of a polygon totally over the plane, return original polygon
0057 if sum(~below) == 0
0058     poly2 = poly;
0059     return;
0060 end
0061 
0062 % number of intersections
0063 nInter = sum(abs(diff(below)));
0064 
0065 % number of vertices of new polygon
0066 N   = size(poly, 1);
0067 % N2  = sum(below(1:end-1)) + nInter;
0068 N2  = sum(below) + nInter;
0069 poly2 = zeros(N2, 3);
0070 
0071 
0072 %% Iteration on polygon vertices
0073 
0074 % vertex index in current polygon
0075 % initialized with first vertex below the plane (vertices before are drop)
0076 i = find(below, 1, 'first');
0077 
0078 % vertex index in result polygon
0079 j = 1;
0080 
0081 while i <= N
0082     
0083     if below(i)
0084         % keep points located below the plane
0085         poly2(j, :) = poly(i,:);
0086         i = i + 1;
0087         j = j + 1;
0088 
0089     else
0090         % current vertex is above the plane. We know that previous vertex
0091         % was below. We compute intersection of supporting line, find the
0092         % next vertex below, and find next intersection.
0093         
0094         % compute intersection of current edge with plane
0095         line = createLine3d(poly(i-1, :), poly(i, :));
0096         inter1 = intersectLinePlane(line, plane);
0097         poly2(j, :) = inter1;
0098         j = j + 1;
0099         
0100         % find index of next vertex below the plane, possibily re-starting
0101         % from the beginning of the polygon
0102         while ~below(mod(i - 1, N) + 1)
0103             i = i + 1;
0104         end
0105         
0106         % compute intersection of current line with plane
0107         i2 = mod(i - 1, N) + 1;
0108         line = createLine3d(poly(i2-1, :), poly(i2, :));
0109         inter2 = intersectLinePlane(line, plane);
0110         poly2(j, :) = inter2;
0111         j = j + 1;
0112         
0113         % add also the current vertex
0114         poly2(j, :) = poly(i2, :);
0115         j = j + 1;
0116         i = i + 1;
0117     end
0118 end
0119 
0120 % remove last point if it is the same as the first one
0121 if sum(poly2(end, :) == poly2(1,:)) == 3
0122     poly2(end, :) = [];
0123 end