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, F] = createCube;
P = createPlane([.5 .5 .5], [1 1 1]);
[N2, F2] = clipConvexPolyhedronHP(N, F, P);
figure('color','w'); view(3); axis equal
drawPolyhedron(N2, F2);
[v, f] = createSoccerBall;
p = createPlane([-.5 .5 -.5], [1 1 1]);
[v2, f2] = clipConvexPolyhedronHP(v, f, p);
figure('color','w'); view(3); axis equal
drawMesh(v, f, 'faceColor', 'none');
drawMesh(v2, f2);
See also
meshes3d, polyhedra, planes3d0001 function [nodes2, faces2] = clipConvexPolyhedronHP(nodes, faces, plane) 0002 %CLIPCONVEXPOLYHEDRONHP Clip a convex polyhedron by a plane. 0003 % 0004 % [NODES2, FACES2] = clipConvexPolyhedronHP(NODES, FACES, PLANE) 0005 % 0006 % return the new (convex) polyhedron whose vertices are 'below' the 0007 % specified plane, and with faces clipped accordingly. NODES2 contains 0008 % clipped vertices and new created vertices, FACES2 contains references 0009 % to NODES2 vertices. 0010 % 0011 % Example 0012 % [N, F] = createCube; 0013 % P = createPlane([.5 .5 .5], [1 1 1]); 0014 % [N2, F2] = clipConvexPolyhedronHP(N, F, P); 0015 % figure('color','w'); view(3); axis equal 0016 % drawPolyhedron(N2, F2); 0017 % 0018 % [v, f] = createSoccerBall; 0019 % p = createPlane([-.5 .5 -.5], [1 1 1]); 0020 % [v2, f2] = clipConvexPolyhedronHP(v, f, p); 0021 % figure('color','w'); view(3); axis equal 0022 % drawMesh(v, f, 'faceColor', 'none'); 0023 % drawMesh(v2, f2); 0024 % 0025 % See also 0026 % meshes3d, polyhedra, planes3d 0027 % 0028 0029 % ------ 0030 % Author: David Legland 0031 % e-mail: david.legland@nantes.inra.fr 0032 % Created: 2007-09-14, using Matlab 7.4.0.287 (R2007a) 0033 % Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas. 0034 0035 0036 %% Preprocessing 0037 0038 % used for identifying identical vertices 0039 tol = 1e-10; 0040 0041 % if faces is a numeric array, convert it to cell array 0042 if isnumeric(faces) 0043 faces2 = cell(size(faces, 1), 1); 0044 for f = 1:length(faces2) 0045 faces2{f} = faces(f,:); 0046 end 0047 faces = faces2; 0048 end 0049 0050 % find vertices below the plane 0051 b = isBelowPlane(nodes, plane); 0052 0053 % initialize results 0054 Nn = size(nodes, 1); 0055 nodes2 = zeros(0, 3); % list of new nodes 0056 faces2 = faces; % list of new faces. Start with initial list, and remove some of them 0057 0058 0059 %% Main iteration on faces 0060 0061 % iterate on each face, and test if either: 0062 % - all points below plane -> keep all face 0063 % - all points up plane -> remove face 0064 % - both -> clip the polygon 0065 keep = true(length(faces), 1); 0066 for f = 1:length(faces) 0067 % current face 0068 face = faces{f}; 0069 bf = b(face); 0070 0071 % face totally above plane 0072 if sum(bf) == 0 0073 keep(f) = false; 0074 continue; 0075 end 0076 0077 % face totally below plane 0078 if sum(bf == 1) == length(bf) 0079 continue; 0080 end 0081 0082 % clip polygon formed by face 0083 poly = nodes(face, :); 0084 clipped = clipConvexPolygon3dHP(poly, plane); 0085 0086 % identify indices of polygon vertices 0087 inds = zeros(1, size(clipped, 1)); 0088 faceb = face(bf==1); % indices of vertices still in clipped face 0089 [minDists, I] = minDistancePoints(nodes(faceb,:), clipped); %#ok<ASGLU> 0090 for i = 1:length(I) 0091 inds(I(i)) = faceb(i); 0092 end 0093 0094 % indices of new points in clipped polygon 0095 indNews = find(inds == 0); 0096 0097 if size(nodes2, 1) < 2 0098 nodes2 = [nodes2; clipped(indNews, :)]; %#ok<AGROW> 0099 inds(indNews(1)) = Nn + 1; 0100 inds(indNews(2)) = Nn + 2; 0101 faces2{f} = inds; 0102 continue; 0103 end 0104 0105 % distances from new vertices to already added vertices 0106 [minDists, I] = minDistancePoints(clipped(indNews, :), nodes2); 0107 0108 % compute index of first vertex 0109 if minDists(1) < tol 0110 inds(indNews(1)) = Nn + I(1); 0111 else 0112 nodes2 = [nodes2; clipped(indNews(1), :)]; %#ok<AGROW> 0113 inds(indNews(1)) = Nn + size(nodes2, 1); 0114 end 0115 0116 % compute index of second vertex 0117 if minDists(2) < tol 0118 inds(indNews(2)) = Nn + I(2); 0119 else 0120 nodes2 = [nodes2; clipped(indNews(2), :)]; %#ok<AGROW> 0121 inds(indNews(2)) = Nn + size(nodes2, 1); 0122 end 0123 0124 % stores the modified face 0125 faces2{f} = inds; 0126 end 0127 0128 0129 %% Postprocessing 0130 0131 % creates a new face formed by the added nodes 0132 [sortedNodes, I] = angleSort3d(nodes2); 0133 0134 % compute normal vector of new face, and reverse order if face points in 0135 % the opposite direction as plane normal 0136 newFaceNormal = meshFaceNormals(sortedNodes, 1:length(sortedNodes)); 0137 if dot(newFaceNormal, planeNormal(plane)) < 0 0138 I(2:end) = I(end:-1:2); 0139 end 0140 0141 % compute vertex indices of new face 0142 newFace = I' + Nn; 0143 0144 % remove faces outside plane and add the new face 0145 faces2 = {faces2{keep}, newFace}; 0146 0147 % remove clipped nodes, and add new nodes to list of nodes 0148 N2 = size(nodes2, 1); 0149 nodes2 = [nodes(b, :); nodes2]; 0150 0151 % new nodes are inside half-space by definition 0152 b = [b; ones(N2, 1)]; 0153 0154 % create look up table between old indices and new indices 0155 inds = zeros(size(nodes2, 1), 1); 0156 indb = find(b); 0157 for i = 1:length(indb) 0158 inds(indb(i)) = i; 0159 end 0160 0161 % update indices of faces 0162 for f = 1:length(faces2) 0163 face = faces2{f}; 0164 faces2{f} = inds(face)'; 0165 end