Point in Polygon algorithm - Why does it work? I'm using the following algorithm (in C) to find if a point lays within a given polygon
typedef struct {
   int h,v;
} Point;

int InsidePolygon(Point *polygon,int n,Point p)
{
   int i;
   double angle=0;
   Point p1,p2;

   for (i=0;i<n;i++) {
      p1.h = polygon[i].h - p.h;
      p1.v = polygon[i].v - p.v;
      p2.h = polygon[(i+1)%n].h - p.h;
      p2.v = polygon[(i+1)%n].v - p.v;
      angle += Angle2D(p1.h,p1.v,p2.h,p2.v);
   }

   if (ABS(angle) < PI)
      return(FALSE);
   else
      return(TRUE);
}

/*
   Return the angle between two vectors on a plane
   The angle is from vector 1 to vector 2, positive anticlockwise
   The result is between -pi -> pi
*/
double Angle2D(double x1, double y1, double x2, double y2)
{
   double dtheta,theta1,theta2;

   theta1 = atan2(y1,x1);
   theta2 = atan2(y2,x2);
   dtheta = theta2 - theta1;
   while (dtheta > PI)
      dtheta -= TWOPI;
   while (dtheta < -PI)
      dtheta += TWOPI;

   return(dtheta);
}

I found the algorithm while searching online for such an algorithm, Point in Polygon algorithm, Solution 2 (2D).


*

*The explanation in the link says it checks if the sum of angles is $2\pi$ but the algorithm (which works) checks it with $\pi$.

*How are the angles even calculated?

 A: Here is a counter-proposal for efficient computation, using only additions and multiplications.
The idea is to count the intersections of the edges of the polygon with the horizontal line through the test point that lie on its right (like in "Solution 1"). An even number means outside.

For the test point $P$ and the edge $QR$, such an intersection is detected by the condition
$$(Q_y\le P_y\land R_y\ge P_y\land\Delta_{PQR}\ge0)\lor(Q_y\ge P_y\land R_y\le P_y\land\Delta_{PQR}\le0),$$ where $\Delta_{PQR}$ is the signed area of the triangle $PQR$, which tells on what side of $QR$ the point $P$ lies.
$$\Delta_{PQR}=(Q_x-P_x)(R_y-P_y)-(Q_y-P_y)(R_x-P_x).$$ Shortcut evaluation is recommended.
Inside= False
q= N-1; r= 0
while r < N:
    if V[q].Y <= P.Y:
        if V[r].Y > Y and Delta(P, V[q], V[r]) >= 0:
            Inside= not Inside
    else:
        if V[r].Y <= Y and Delta(P, V[q], V[r]) <= 0:
            Inside= not Inside
    q= r; r++

The total cost is $2N$ ordinate comparisons and $M$ tests on the sign of $\Delta$, where $M$ is the number of edges straddling the ordinates of $P$, plus loop overhead.
A: The angle is computed with the atan2 function (see
https://en.wikipedia.org/wiki/Arctan2). But the code is not computationally
robust. The mathematical idea may be OK, but since $\pi$ and $2\pi$ are not exactly
representable as floating point numbers, there will be uncertainties produced by using floating point arithmetic near the interesting sum $2\pi$. Testing equality of two floating point variables is a questionable concept in many situations (since the angle is in the interval  $[-\pi, +\pi],\,$ the test $|\, . |<\pi\,$ is actually an equality test in disguise).
And there is a bug, if the sum happens to be $\pm\pi,$ the code should handle this separately.
