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

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 % Author: David Legland
0026 % E-mail: david.legland@inrae.fr
0027 % Created: 2007-04-24, from intersectLinePlane
0028 % Copyright 2007-2024 INRA - TPV URPOI - BIA IMASTE
0029 
0030 % extract tolerance for determination of parallel edges and planes
0031 tol = 1e-14;
0032 if ~isempty(varargin)
0033     tol = varargin{1};
0034 end
0035 
0036 % number of planes and edges
0037 np = size(plane, 1);
0038 ne = size(edge, 1);
0039 
0040 % unify sizes of data
0041 if np ~= ne
0042     if ne == 1
0043         % one edge and many planes
0044         edge = edge(ones(np, 1), :);
0045     elseif np == 1
0046         % one plane possible many edges
0047         plane = plane(ones(ne, 1), :);
0048     else
0049         % N planes and M edges, not allowed for now.
0050         error('Should have the same number of planes and edges');
0051     end
0052 end
0053 
0054 % initialize empty arrays
0055 point = zeros(size(plane, 1), 3);
0056 t = zeros(size(plane,1),3);
0057 
0058 % plane normal
0059 n = cross(plane(:,4:6), plane(:,7:9), 2);
0060 
0061 % create line supporting edge
0062 line = createLine3d(edge(:,1:3), edge(:,4:6));
0063 
0064 % get indices of edge and plane which are parallel
0065 par = abs(dot(n, line(:,4:6), 2)) < tol;
0066 point(par,:) = NaN;
0067 t(par) = NaN;
0068 
0069 % difference between origins of plane and edge
0070 dp = plane(:, 1:3) - line(:, 1:3);
0071 
0072 % relative position of intersection on line
0073 %t = dot(n(~par,:), dp(~par,:), 2)./dot(n(~par,:), line(~par,4:6), 2);
0074 t(~par) = dot(n(~par,:), dp(~par,:), 2) ./ dot(n(~par,:), line(~par,4:6), 2);
0075 
0076 % compute coord of intersection point
0077 %point(~par, :) = line(~par,1:3) + repmat(t,1,3).*line(~par,4:6);
0078 point(~par, :) = line(~par,1:3) + repmat(t(~par),1,3) .* line(~par,4:6);
0079 
0080 % set points outside of edge to [NaN NaN NaN]
0081 point(t<0, :) = NaN;
0082 point(t>1, :) = NaN;