dihedralAngle

DIHEDRALANGLE Compute dihedral angle between 2 planes.

   THETA = dihedralAngle(PLANE1, PLANE2)
   PLANE1 and PLANE2 are plane representations given in the following
   format:
   [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
   THETA is the angle between the two vectors given by plane normals,
   given between 0 and PI.

   References
   http://en.wikipedia.org/wiki/Dihedral_angle
   http://mathworld.wolfram.com/DihedralAngle.html

   See also:
   planes3d, lines3d, angles3d, planesBisector

   ---------
   author : David Legland
   INRA - TPV URPOI - BIA IMASTE
   created the 21/02/2005.

0001 function theta = dihedralAngle(plane1, plane2)
0002 %DIHEDRALANGLE Compute dihedral angle between 2 planes.
0003 %
0004 %   THETA = dihedralAngle(PLANE1, PLANE2)
0005 %   PLANE1 and PLANE2 are plane representations given in the following
0006 %   format:
0007 %   [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
0008 %   THETA is the angle between the two vectors given by plane normals,
0009 %   given between 0 and PI.
0010 %
0011 %   References
0012 %   http://en.wikipedia.org/wiki/Dihedral_angle
0013 %   http://mathworld.wolfram.com/DihedralAngle.html
0014 %
0015 %   See also:
0016 %   planes3d, lines3d, angles3d, planesBisector
0017 %
0018 %   ---------
0019 %   author : David Legland
0020 %   INRA - TPV URPOI - BIA IMASTE
0021 %   created the 21/02/2005.
0022 %
0023 
0024 % HISTORY
0025 % 2009-06-19 change convention for dihedral angle
0026 % 2011-03-20 improve computation precision
0027 
0028 % compute plane normals
0029 v1 = planeNormal(plane1);
0030 v2 = planeNormal(plane2);
0031 
0032 % number of vectors
0033 n1 = size(v1, 1);
0034 n2 = size(v2, 1);
0035 
0036 % ensures vectors have same dimension
0037 if n1 ~= n2
0038     if n1 == 1
0039         v1 = repmat(v1, [n2 1]);
0040     elseif n2 == 1
0041         v2 = repmat(v2, [n1 1]);
0042     else
0043         error('Arguments V1 and V2 must have the same size');
0044     end
0045 end
0046 
0047 % compute dihedral angle(s)
0048 theta = atan2(vectorNorm3d(cross(v1, v2, 2)), dot(v1, v2, 2));
0049 
0050 % % equivalent to following formula, but more precise for small angles:
0051 % n1 = normalizeVector3d(planeNormal(plane1));
0052 % n2 = normalizeVector3d(planeNormal(plane2));
0053 % theta = acos(dot(n1, n2, 2));
0054