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

clipConvexPolyhedronHP(nodes, faces, plane)
```function [nodes2, faces2] = clipConvexPolyhedronHP(nodes, faces, plane)
%CLIPCONVEXPOLYHEDRONHP Clip a convex polyhedron by a plane
%
%   [NODES2, FACES2] = clipConvexPolyhedronHP(NODES, FACES, PLANE)
%
%   return the new (convex) polyhedron whose vertices are 'below' the
%   specified plane, and with faces clipped accordingly. NODES2 contains
%   clipped vertices and new created vertices, FACES2 contains references
%   to NODES2 vertices.
%
%   Example
%   [N E F] = createCube;
%   P = createPlane([.5 .5 .5], [1 1 1]);
%   [N2 F2] = clipConvexPolyhedronHP(N, F, P);
%   drawPolyhedra(N2, F2);
%
%   meshes3d, polyhedra, planes3d
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2007-09-14,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Preprocessing

% used for identifying identical vertices
tol = 1e-10;

% if faces is a numeric array, convert it to cell array
if isnumeric(faces)
faces2 = cell(size(faces, 1), 1);
for f = 1:length(faces2)
faces2{f} = faces(f,:);
end
faces = faces2;
end

% find vertices below the plane
b = isBelowPlane(nodes, plane);

% initialize results
Nn  = size(nodes, 1);
nodes2 = zeros(0, 3);   % list of new nodes
faces2 = faces;         % list of new faces. Start with initial list, and remove some of them

%% Main iteration on faces

% iterate on each face, and test if either:
%   - all points below plane -> keep all face
%   - all points up plane -> remove face
%   - both -> clip the polygon
keep = true(length(faces), 1);
for f = 1:length(faces)
% current face
face = faces{f};
bf = b(face);

% face totally above plane
if sum(bf) == 0
keep(f) = false;
continue;
end

% face totally below plane
if sum(bf == 1) == length(bf)
continue;
end

% clip polygon formed by face
poly = nodes(face, :);
clipped = clipConvexPolygon3dHP(poly, plane);

% identify indices of polygon vertices
inds = zeros(1, size(clipped, 1));
faceb = face(bf==1); % indices of vertices still in clipped face
[minDists, I] = minDistancePoints(nodes(faceb,:), clipped); %#ok<ASGLU>
for i = 1:length(I)
inds(I(i)) = faceb(i);
end

% indices of new points in clipped polygon
indNews = find(inds == 0);

if size(nodes2, 1) < 2
nodes2 = [nodes2; clipped(indNews, :)]; %#ok<AGROW>
inds(indNews(1)) = Nn + 1;
inds(indNews(2)) = Nn + 2;
faces2{f} = inds;
continue;
end

[minDists, I] = minDistancePoints(clipped(indNews, :), nodes2);

% compute index of first vertex
if minDists(1) < tol
inds(indNews(1)) = Nn + I(1);
else
nodes2 = [nodes2; clipped(indNews(1), :)]; %#ok<AGROW>
inds(indNews(1)) = Nn + size(nodes2, 1);
end

% compute index of second vertex
if minDists(2) < tol
inds(indNews(2)) = Nn + I(2);
else
nodes2 = [nodes2; clipped(indNews(2), :)]; %#ok<AGROW>
inds(indNews(2)) = Nn + size(nodes2, 1);
end

% stores the modified face
faces2{f} = inds;
end

%% Postprocessing

% creates a new face formed by the added nodes
[tmp, I] = angleSort3d(nodes2); %#ok<ASGLU>
newFace = I' + Nn;

% remove faces outside plane and add the new face
faces2 = {faces2{keep}, newFace};

% remove clipped nodes, and add new nodes to list of nodes
N2 = size(nodes2, 1);
nodes2 = [nodes(b, :); nodes2];

% new nodes are inside half-space by definition
b = [b; ones(N2, 1)];

% create look up table between old indices and new indices
inds = zeros(size(nodes2, 1), 1);
indb = find(b);
for i = 1:length(indb)
inds(indb(i)) = i;
end

% update indices of faces
for f = 1:length(faces2)
face = faces2{f};
faces2{f} = inds(face)';
end
```