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 %   INRA - TPV URPOI - BIA IMASTE
0020 %   created the 17/02/2005.
0021 %
0022 
0023 %   HISTORY
0024 
0025 tol = 1e-14;
0026 if ~isempty(varargin)
0027     tol = varargin{1};
0028 end
0029 
0030 % plane normal
0031 n1 = normalizeVector3d(cross(plane1(:,4:6), plane1(:, 7:9), 2));
0032 n2 = normalizeVector3d(cross(plane2(:,4:6), plane2(:, 7:9), 2));
0033 
0034 % test if planes are parallel
0035 if abs(cross(n1, n2, 2)) < tol
0036     line = [NaN NaN NaN NaN NaN NaN];
0037     return;
0038 end
0039 
0040 % Uses Hessian form, ie : N.p = d
0041 % I this case, d can be found as : -N.p0, when N is normalized
0042 d1 = dot(n1, plane1(:,1:3), 2);
0043 d2 = dot(n2, plane2(:,1:3), 2);
0044 
0045 % compute dot products
0046 dot1 = dot(n1, n1, 2);
0047 dot2 = dot(n2, n2, 2);
0048 dot12 = dot(n1, n2, 2);
0049 
0050 % intermediate computations
0051 det = dot1*dot2 - dot12*dot12;
0052 c1  = (d1*dot2 - d2*dot12)./det;
0053 c2  = (d2*dot1 - d1*dot12)./det;
0054 
0055 % compute line origin and direction
0056 p0  = c1*n1 + c2*n2;
0057 dp  = cross(n1, n2, 2);
0058 
0059 line = [p0 dp];