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.

   Example
       rad2deg(anglePoints3d([0 0 1],[1 1 0]))

   See also
   points3d, angles3d

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 %   Example
0014 %       rad2deg(anglePoints3d([0 0 1],[1 1 0]))
0015 %
0016 %   See also
0017 %   points3d, angles3d
0018 
0019 % ------
0020 % Author: David Legland
0021 % E-mail: david.legland@inrae.fr
0022 % Created: 2005-02-21
0023 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0024 
0025 p2 = [0 0 0];
0026 if isscalar(varargin)
0027     pts = varargin{1};
0028     if size(pts, 1)==2
0029         p1 = pts(1,:);
0030         p0 = [0 0 0];
0031         p2 = pts(2,:);
0032     else
0033         p1 = pts(1,:);
0034         p0 = pts(2,:);
0035         p2 = pts(3,:);
0036     end
0037     
0038 elseif length(varargin) == 2
0039     p1 = varargin{1};
0040     p0 = [0 0 0];
0041     p2 = varargin{2};
0042     
0043 elseif length(varargin) == 3
0044     p1 = varargin{1};
0045     p0 = varargin{2};
0046     p2 = varargin{3};
0047 end
0048 
0049 % ensure all data have same size
0050 n1 = size(p1, 1);
0051 n2 = size(p2, 1);
0052 n0 = size(p0, 1);
0053 
0054 if n1 ~= n0
0055     if n1 == 1
0056         p1 = repmat(p1, [n0 1]);
0057         n1 = n0;
0058     elseif n0==1
0059         p0 = repmat(p0, [n1 1]);
0060     else
0061         error('Arguments P1 and P0 must have the same size');
0062     end
0063 end
0064 
0065 if n1 ~= n2
0066     if n1 == 1
0067         p1 = repmat(p1, [n2 1]);
0068     elseif n2 == 1
0069         p2 = repmat(p2, [n1 1]);
0070     else
0071         error('Arguments P1 and P2 must have the same size');
0072     end
0073 end
0074 
0075 % normalized vectors
0076 p1 = normalizeVector3d(p1 - p0);
0077 p2 = normalizeVector3d(p2 - p0);
0078 
0079 % compute angle
0080 alpha = acos(dot(p1, p2, 2));