anglePoints3d

ANGLEPOINTS3D Compute angle between three 3D points.

   ALPHA = anglePoints3d(P1, P2)
   Computes angle (P1, O, P2), in radians, between 0 and PI.

   ALPHA = anglePoints3d(P1, P2, P3)
   Computes angle (P1, P2, P3), in radians, between 0 and PI.

   ALPHA = anglePoints3d(PTS)
   PTS is a 3x3 or 2x3 array containing coordinate of points.

   See also
   points3d, angles3d

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

0001 function alpha = anglePoints3d(varargin)
0002 %ANGLEPOINTS3D Compute angle between three 3D points.
0003 %
0004 %   ALPHA = anglePoints3d(P1, P2)
0005 %   Computes angle (P1, O, P2), in radians, between 0 and PI.
0006 %
0007 %   ALPHA = anglePoints3d(P1, P2, P3)
0008 %   Computes angle (P1, P2, P3), in radians, between 0 and PI.
0009 %
0010 %   ALPHA = anglePoints3d(PTS)
0011 %   PTS is a 3x3 or 2x3 array containing coordinate of points.
0012 %
0013 %   See also
0014 %   points3d, angles3d
0015 %
0016 %   ---------
0017 %   author : David Legland
0018 %   INRA - TPV URPOI - BIA IMASTE
0019 %   created the 21/02/2005.
0020 %
0021 
0022 %   HISTORY
0023 %   20/09/2005: add case of single argument for all points
0024 %   04/01/2007: check typo
0025 %   27/05/2014: adjust known vector sizes n1, n0, n2 once corrected for
0026 
0027 
0028 p2 = [0 0 0];
0029 if length(varargin) == 1
0030     pts = varargin{1};
0031     if size(pts, 1)==2
0032         p1 = pts(1,:);
0033         p0 = [0 0 0];
0034         p2 = pts(2,:);
0035     else
0036         p1 = pts(1,:);
0037         p0 = pts(2,:);
0038         p2 = pts(3,:);
0039     end
0040     
0041 elseif length(varargin) == 2
0042     p1 = varargin{1};
0043     p0 = [0 0 0];
0044     p2 = varargin{2};
0045     
0046 elseif length(varargin) == 3
0047     p1 = varargin{1};
0048     p0 = varargin{2};
0049     p2 = varargin{3};
0050 end
0051 
0052 % ensure all data have same size
0053 n1 = size(p1, 1);
0054 n2 = size(p2, 1);
0055 n0 = size(p0, 1);
0056 
0057 if n1 ~= n0
0058     if n1 == 1
0059         p1 = repmat(p1, [n0 1]);
0060         n1 = n0;
0061     elseif n0==1
0062         p0 = repmat(p0, [n1 1]);
0063     else
0064         error('Arguments P1 and P0 must have the same size');
0065     end
0066 end
0067 
0068 if n1 ~= n2
0069     if n1 == 1
0070         p1 = repmat(p1, [n2 1]);
0071     elseif n2 == 1
0072         p2 = repmat(p2, [n1 1]);
0073     else
0074         error('Arguments P1 and P2 must have the same size');
0075     end
0076 end
0077 
0078 % normalized vectors
0079 p1 = normalizeVector3d(p1 - p0);
0080 p2 = normalizeVector3d(p2 - p0);
0081 
0082 % compute angle
0083 alpha = acos(dot(p1, p2, 2));