intersectPlanes

INTERSECTPLANES Return intersection line between 2 planes in space.

   LINE = intersectPlanes(PLANE1, PLANE2)
   Returns the straight line belonging to both planes.
   PLANE:  [x0 y0 z0  dx1 dy1 dz1  dx2 dy2 dz2]
   LINE:   [x0 y0 z0  dx dy dz]
   In case of parallel planes, returns line with NaN values.

   LINE = intersectPlanes(PLANE1, PLANE2, TOL)
   Also specifies the tolerance for detecting parallel planes.

   See also
   planes3d, lines3d, intersectLinePlane

0001 function line = intersectPlanes(plane1, plane2, varargin)
0002 %INTERSECTPLANES Return intersection line between 2 planes in space.
0003 %
0004 %   LINE = intersectPlanes(PLANE1, PLANE2)
0005 %   Returns the straight line belonging to both planes.
0006 %   PLANE:  [x0 y0 z0  dx1 dy1 dz1  dx2 dy2 dz2]
0007 %   LINE:   [x0 y0 z0  dx dy dz]
0008 %   In case of parallel planes, returns line with NaN values.
0009 %
0010 %   LINE = intersectPlanes(PLANE1, PLANE2, TOL)
0011 %   Also specifies the tolerance for detecting parallel planes.
0012 %
0013 %   See also
0014 %   planes3d, lines3d, intersectLinePlane
0015 %
0016 
0017 % ------
0018 % Author: David Legland
0019 % E-mail: david.legland@inrae.fr
0020 % Created: 2005-02-17
0021 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0022 
0023 tol = 1e-14;
0024 if ~isempty(varargin)
0025     tol = varargin{1};
0026 end
0027 
0028 % plane normal
0029 n1 = normalizeVector3d(cross(plane1(:,4:6), plane1(:, 7:9), 2));
0030 n2 = normalizeVector3d(cross(plane2(:,4:6), plane2(:, 7:9), 2));
0031 
0032 % test if planes are parallel
0033 if abs(cross(n1, n2, 2)) < tol
0034     line = [NaN NaN NaN NaN NaN NaN];
0035     return;
0036 end
0037 
0038 % Uses Hessian form, ie : N.p = d
0039 % I this case, d can be found as : -N.p0, when N is normalized
0040 d1 = dot(n1, plane1(:,1:3), 2);
0041 d2 = dot(n2, plane2(:,1:3), 2);
0042 
0043 % compute dot products
0044 dot1 = dot(n1, n1, 2);
0045 dot2 = dot(n2, n2, 2);
0046 dot12 = dot(n1, n2, 2);
0047 
0048 % intermediate computations
0049 det = dot1*dot2 - dot12*dot12;
0050 c1  = (d1*dot2 - d2*dot12)./det;
0051 c2  = (d2*dot1 - d1*dot12)./det;
0052 
0053 % compute line origin and direction
0054 p0  = c1*n1 + c2*n2;
0055 dp  = cross(n1, n2, 2);
0056 
0057 line = [p0 dp];