planesBisector

 PLANESBISECTOR  Bisector plane between two other planes.
 
   BIS = planesBisector(PL1, PL2);
   Returns the planes that contains the intersection line between PL1 and
   PL2 and that bisect the dihedral angle of PL1 and PL2. 
   Note that computing the bisector of PL2 and PL1 (in that order) returns
   the same plane but with opposite orientation.

   Example
     % Draw two planes together with their bisector
     pl1 = createPlane([3 4 5], [1 2 3]);
     pl2 = createPlane([3 4 5], [2 -3 4]);
     % compute bisector
     bis = planesBisector(pl1, pl2);
     % setup display
     figure; hold on; axis([0 10 0 10 0 10]);
     set(gcf, 'renderer', 'opengl')
     view(3);
     % draw the planes
     drawPlane3d(pl1, 'g');
     drawPlane3d(pl2, 'g');
     drawPlane3d(bis, 'b');

   See also
   planes3d, dihedralAngle, intersectPlanes

   Author: Ben X. Kang
   Dept. Orthopaedics & Traumatology
   Li Ka Shing Faculty of Medicine
   The University of Hong Kong
   Pok Fu Lam, Hong Kong

0001 function out = planesBisector(plane1, plane2)
0002 % PLANESBISECTOR  Bisector plane between two other planes.
0003 %
0004 %   BIS = planesBisector(PL1, PL2);
0005 %   Returns the planes that contains the intersection line between PL1 and
0006 %   PL2 and that bisect the dihedral angle of PL1 and PL2.
0007 %   Note that computing the bisector of PL2 and PL1 (in that order) returns
0008 %   the same plane but with opposite orientation.
0009 %
0010 %   Example
0011 %     % Draw two planes together with their bisector
0012 %     pl1 = createPlane([3 4 5], [1 2 3]);
0013 %     pl2 = createPlane([3 4 5], [2 -3 4]);
0014 %     % compute bisector
0015 %     bis = planesBisector(pl1, pl2);
0016 %     % setup display
0017 %     figure; hold on; axis([0 10 0 10 0 10]);
0018 %     set(gcf, 'renderer', 'opengl')
0019 %     view(3);
0020 %     % draw the planes
0021 %     drawPlane3d(pl1, 'g');
0022 %     drawPlane3d(pl2, 'g');
0023 %     drawPlane3d(bis, 'b');
0024 %
0025 %   See also
0026 %   planes3d, dihedralAngle, intersectPlanes
0027 %
0028 %   Author: Ben X. Kang
0029 %   Dept. Orthopaedics & Traumatology
0030 %   Li Ka Shing Faculty of Medicine
0031 %   The University of Hong Kong
0032 %   Pok Fu Lam, Hong Kong
0033 %
0034 
0035 % Let the two planes be defined by equations
0036 %
0037 %  a1*x + b1*y + c1*z + d1 = 0
0038 %
0039 % and
0040 %
0041 %  a2*x + b2*y + c2*z + d2 = 0
0042 %
0043 % in which vectors [a1,b1,c1] and [a2,b2,c2] are normalized to be of unit
0044 % length (a^2+b^2+c^2 = 1). Then
0045 %
0046 %  (a1+a2)*x + (b1+b2)*y + (c1+c2)*z + (d1+d2) = 0
0047 %
0048 % is the equation of the desired plane which bisects the dihedral angle
0049 % between the two planes.  These coefficients cannot be all zero because
0050 % the two given planes are not parallel.
0051 %
0052 % Notice that there is a second solution to this problem
0053 %
0054 %  (a1-a2)*x + (b1-b2)*y + (c1-c2)*z + (d1-d2) = 0
0055 %
0056 % which also is a valid plane and orthogonal to the first solution. One of
0057 % these planes bisects the acute dihedral angle, and the other the
0058 % supplementary obtuse dihedral angle, between the two given planes.
0059 
0060 
0061 P1 = plane1(1:3);            % a point on the plane
0062 n1 = planeNormal(plane1);    % the normal of the plane
0063 % d1 = -dot(n1, P1);        % for line equation
0064 
0065 P2 = plane2(1:3);
0066 n2 = planeNormal(plane2);
0067 % d2 = -dot(n2, P2);
0068 
0069 if ~isequal(P1(1:3), P2(1:3))
0070     L = intersectPlanes(plane1, plane2);    % intersection of the given two planes
0071     Pt = L(1:3);                            % a point on the line intersection
0072 %     v2 = cross(n1-n2, L(4:6));                % another vector lie on the bisect plane
0073 %     out = [v1, v2]';
0074 else
0075     Pt = P1(1:3);
0076 end
0077 
0078 % use column-wise vector
0079 out = createPlane(Pt, n1 - n2);
0080 
0081 
0082 %%  EOF  %%