intersectEdgePlane

INTERSECTEDGEPLANE Return intersection point between a plane and a edge.

   PT = intersectEdgePlane(edge, PLANE) return the intersection point of
   the given edge and the given plane.
   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
   edge :  [x1 y1 z1 x2 y2 z2]
   PT :    [xi yi zi]
   If EDGE and PLANE are parallel, return [NaN NaN NaN].
   If EDGE (or PLANE) is a matrix with 6 (or 9) columns and N rows, result
   is an array of points with N rows and 3 columns.
   
   Example:
   edge = [5 5 -1 5 5 1];
   plane = [0 0 0 1 0 0 0 1 0];
   intersectEdgePlane(edge, plane)     % should return [5 5 0].
   ans =
       5   5   0

   See Also:
   planes3d, intersectLinePlane, createLine3d, createPlane


   ---------

   author : David Legland
   INRA - TPV URPOI - BIA IMASTE
   created the 24/04/2007 from intersectLinePlane.

0001 function point = intersectEdgePlane(edge, plane, varargin)
0002 %INTERSECTEDGEPLANE Return intersection point between a plane and a edge.
0003 %
0004 %   PT = intersectEdgePlane(edge, PLANE) return the intersection point of
0005 %   the given edge and the given plane.
0006 %   PLANE : [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
0007 %   edge :  [x1 y1 z1 x2 y2 z2]
0008 %   PT :    [xi yi zi]
0009 %   If EDGE and PLANE are parallel, return [NaN NaN NaN].
0010 %   If EDGE (or PLANE) is a matrix with 6 (or 9) columns and N rows, result
0011 %   is an array of points with N rows and 3 columns.
0012 %
0013 %   Example:
0014 %   edge = [5 5 -1 5 5 1];
0015 %   plane = [0 0 0 1 0 0 0 1 0];
0016 %   intersectEdgePlane(edge, plane)     % should return [5 5 0].
0017 %   ans =
0018 %       5   5   0
0019 %
0020 %   See Also:
0021 %   planes3d, intersectLinePlane, createLine3d, createPlane
0022 %
0023 %
0024 %   ---------
0025 %
0026 %   author : David Legland
0027 %   INRA - TPV URPOI - BIA IMASTE
0028 %   created the 24/04/2007 from intersectLinePlane.
0029 %
0030 
0031 %   HISTORY
0032 %
0033 %   17/06/2011 E. J. Payton - fixed indexing error that caused incorrect
0034 %              points to be returned
0035 
0036 % extract tolerance for determination of parallel edges and planes
0037 tol = 1e-14;
0038 if ~isempty(varargin)
0039     tol = varargin{1};
0040 end
0041 
0042 % number of planes and edges
0043 np = size(plane, 1);
0044 ne = size(edge, 1);
0045 
0046 % unify sizes of data
0047 if np ~= ne
0048     if ne == 1
0049         % one edge and many planes
0050         edge = edge(ones(np, 1), :);
0051     elseif np == 1
0052         % one plane possible many edges
0053         plane = plane(ones(ne, 1), :);
0054     else
0055         % N planes and M edges, not allowed for now.
0056         error('Should have the same number of planes and edges');
0057     end
0058 end
0059 
0060 % initialize empty arrays
0061 point = zeros(size(plane, 1), 3);
0062 t = zeros(size(plane,1),3);
0063 
0064 % plane normal
0065 n = cross(plane(:,4:6), plane(:,7:9), 2);
0066 
0067 % create line supporting edge
0068 line = createLine3d(edge(:,1:3), edge(:,4:6));
0069 
0070 % get indices of edge and plane which are parallel
0071 par = abs(dot(n, line(:,4:6), 2)) < tol;
0072 point(par,:) = NaN;
0073 t(par) = NaN;
0074 
0075 % difference between origins of plane and edge
0076 dp = plane(:, 1:3) - line(:, 1:3);
0077 
0078 % relative position of intersection on line
0079 %t = dot(n(~par,:), dp(~par,:), 2)./dot(n(~par,:), line(~par,4:6), 2);
0080 t(~par) = dot(n(~par,:), dp(~par,:), 2) ./ dot(n(~par,:), line(~par,4:6), 2);
0081 
0082 % compute coord of intersection point
0083 %point(~par, :) = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);
0084 point(~par, :) = line(~par,1:3) + repmat(t(~par),1,3) .* line(~par,4:6);
0085 
0086 % set points outside of edge to [NaN NaN NaN]
0087 point(t<0, :) = NaN;
0088 point(t>1, :) = NaN;