Code covered by the BSD License

# geom3d

### David Legland (view profile)

19 Jun 2009 (Updated )

Library to handle 3D geometric primitives: create, intersect, display, and make basic computations

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

clipPolygon3dHP(poly, plane)
```function poly2 = clipPolygon3dHP(poly, plane)
%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'
%
%   poygons3d, polyhedra, clipConvexPolygon3dHP
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 02/08/2005.
%

%   HISTORY
%   2011-08-17 rewrite algo, that works for convex polygons, but is slower
%       than function clipConvexPolgon3dHP

%% Pre-Processing

% ensure last point is the same as the first one (makes computation easier)
if sum(poly(end, :) == poly(1,:)) ~= 3
poly = [poly; poly(1,:)];
end

% compute index of position wrt plane for each vertex
below = isBelowPlane(poly, plane);

% in the case of a polygon totally over the plane, return empty array
if sum(below) == 0
poly2 = zeros(0, 3);
return;
end

% in the case of a polygon totally over the plane, return original polygon
if sum(~below) == 0
poly2 = poly;
return;
end

% number of intersections
nInter = sum(abs(diff(below)));

% number of vertices of new polygon
N   = size(poly, 1);
% N2  = sum(below(1:end-1)) + nInter;
N2  = sum(below) + nInter;
poly2 = zeros(N2, 3);

%% Iteration on polygon vertices

% vertex index in current polygon
% initialized with first vertex below the plane (vertices before are drop)
i = find(below, 1, 'first');

% vertex index in result polygon
j = 1;

while i <= N

if below(i)
% keep points located below the plane
poly2(j, :) = poly(i,:);
i = i + 1;
j = j + 1;

else
% current vertex is above the plane. We know that previous vertex
% was below. We compute intersection of supporting line, find the
% next vertex below, and find next intersection.

% compute intersection of current edge with plane
line = createLine3d(poly(i-1, :), poly(i, :));
inter1 = intersectLinePlane(line, plane);
poly2(j, :) = inter1;
j = j + 1;

% find index of next vertex below the plane, possibily re-starting
% from the beginning of the polygon
while ~below(mod(i - 1, N) + 1)
i = i + 1;
end

% compute intersection of current line with plane
i2 = mod(i - 1, N) + 1;
line = createLine3d(poly(i2-1, :), poly(i2, :));
inter2 = intersectLinePlane(line, plane);
poly2(j, :) = inter2;
j = j + 1;

% add also the current vertex
poly2(j, :) = poly(i2, :);
j = j + 1;
i = i + 1;
end
end

% remove last point if it is the same as the first one
if sum(poly2(end, :) == poly2(1,:)) == 3
poly2(end, :) = [];
end
```