isCounterClockwise

 Compute the relative orientation of 3 points.

   CCW = isCounterClockwise(P1, P2, P3);
   Computes the orientation of the 3 points. The returns is:
   +1 if the path P1->P2->P3 turns Counter-Clockwise (i.e., the point P3
       is located "on the left" of the line P1-P2)
   -1 if the path turns Clockwise (i.e., the point P3 lies "on the right"
       of the line P1-P2) 
   0  if the point P3 is located on the line segment [P1 P2].

   This function can be used in more complicated algorithms: detection of
   line segment intersections, convex hulls, point in triangle...

   CCW = isCounterClockwise(P1, P2, P3, EPS);
   Specifies the threshold used for detecting colinearity of the 3 points.
   Default value is 1e-12 (absolute).

   Example
   isCounterClockwise([0 0], [10 0], [10 10])
   ans = 
       1
   isCounterClockwise([0 0], [0 10], [10 10])
   ans = 
       -1
   isCounterClockwise([0 0], [10 0], [5 0])
   ans = 
       0

   See also
   points2d, isPointOnLine, isPointInTriangle, polygonArea

   References
     Algorithm adapated from Sedgewick's book.

0001 function res = isCounterClockwise(p1, p2, p3, varargin)
0002 % Compute the relative orientation of 3 points.
0003 %
0004 %   CCW = isCounterClockwise(P1, P2, P3);
0005 %   Computes the orientation of the 3 points. The returns is:
0006 %   +1 if the path P1->P2->P3 turns Counter-Clockwise (i.e., the point P3
0007 %       is located "on the left" of the line P1-P2)
0008 %   -1 if the path turns Clockwise (i.e., the point P3 lies "on the right"
0009 %       of the line P1-P2)
0010 %   0  if the point P3 is located on the line segment [P1 P2].
0011 %
0012 %   This function can be used in more complicated algorithms: detection of
0013 %   line segment intersections, convex hulls, point in triangle...
0014 %
0015 %   CCW = isCounterClockwise(P1, P2, P3, EPS);
0016 %   Specifies the threshold used for detecting colinearity of the 3 points.
0017 %   Default value is 1e-12 (absolute).
0018 %
0019 %   Example
0020 %   isCounterClockwise([0 0], [10 0], [10 10])
0021 %   ans =
0022 %       1
0023 %   isCounterClockwise([0 0], [0 10], [10 10])
0024 %   ans =
0025 %       -1
0026 %   isCounterClockwise([0 0], [10 0], [5 0])
0027 %   ans =
0028 %       0
0029 %
0030 %   See also
0031 %   points2d, isPointOnLine, isPointInTriangle, polygonArea
0032 %
0033 %   References
0034 %     Algorithm adapated from Sedgewick's book.
0035 %
0036 
0037 % ------
0038 % Author: David Legland
0039 % e-mail: david.legland@inrae.fr
0040 % Created: 2010-04-09
0041 % Copyright 2011 INRAE - Cepia Software Platform.
0042 
0043 
0044 %   HISTORY
0045 %   2011-05-16 change variable names, add support for point arrays
0046 
0047 
0048 % get threshold value
0049 eps = 1e-12;
0050 if ~isempty(varargin)
0051     eps = varargin{1};
0052 end
0053 
0054 % ensure all data have same size
0055 np = max([size(p1, 1) size(p2, 1) size(p3,1)]);
0056 if np > 1
0057     if size(p1,1) == 1
0058         p1 = repmat(p1, np, 1);
0059     end
0060     if size(p2,1) == 1
0061         p2 = repmat(p2, np, 1);
0062     end
0063     if size(p3,1) == 1
0064         p3 = repmat(p3, np, 1);
0065     end    
0066 end
0067 
0068 % init with 0
0069 res = zeros(np, 1);
0070 
0071 % extract vector coordinates
0072 x0  = p1(:, 1);
0073 y0  = p1(:, 2);
0074 dx1 = p2(:, 1) - x0;
0075 dy1 = p2(:, 2) - y0;
0076 dx2 = p3(:, 1) - x0;
0077 dy2 = p3(:, 2) - y0;
0078 
0079 % check non colinear cases
0080 res(dx1 .* dy2 > dy1 .* dx2) =  1;
0081 res(dx1 .* dy2 < dy1 .* dx2) = -1;
0082 
0083 % case of colinear points
0084 ind = abs(dx1 .* dy2 - dy1 .* dx2) < eps;
0085 res(ind( (dx1(ind) .* dx2(ind) < 0) | (dy1(ind) .* dy2(ind) < 0) )) = -1;
0086 res(ind(  hypot(dx1(ind), dy1(ind)) <  hypot(dx2(ind), dy2(ind)) )) =  1;