isCounterClockwise

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 %ISCOUNTERCLOCKWISE 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 2010-2024 INRAE - Cepia Software Platform
0042 
0043 % get threshold value
0044 eps = 1e-12;
0045 if ~isempty(varargin)
0046     eps = varargin{1};
0047 end
0048 
0049 % ensure all data have same size
0050 np = max([size(p1, 1) size(p2, 1) size(p3,1)]);
0051 if np > 1
0052     if size(p1,1) == 1
0053         p1 = repmat(p1, np, 1);
0054     end
0055     if size(p2,1) == 1
0056         p2 = repmat(p2, np, 1);
0057     end
0058     if size(p3,1) == 1
0059         p3 = repmat(p3, np, 1);
0060     end    
0061 end
0062 
0063 % init with 0
0064 res = zeros(np, 1);
0065 
0066 % extract vector coordinates
0067 x0  = p1(:, 1);
0068 y0  = p1(:, 2);
0069 dx1 = p2(:, 1) - x0;
0070 dy1 = p2(:, 2) - y0;
0071 dx2 = p3(:, 1) - x0;
0072 dy2 = p3(:, 2) - y0;
0073 
0074 % check non colinear cases
0075 res(dx1 .* dy2 > dy1 .* dx2) =  1;
0076 res(dx1 .* dy2 < dy1 .* dx2) = -1;
0077 
0078 % case of colinear points
0079 ind = abs(dx1 .* dy2 - dy1 .* dx2) < eps;
0080 res(ind( (dx1(ind) .* dx2(ind) < 0) | (dy1(ind) .* dy2(ind) < 0) )) = -1;
0081 res(ind(  hypot(dx1(ind), dy1(ind)) <  hypot(dx2(ind), dy2(ind)) )) =  1;