Zoj 1081 Points Within-improved arc-length method for determining the relationship between Points and simple polygon

The relationship between points and polygon through the corner method is known to all. The principle is relatively simple. The angle swept inside the polygon must be 360 degrees, but not necessarily on the boundary and outside.
It is also difficult to implement, and the floating point error is large, and some special cases should be considered.
I found an algorithm called the improved arc-length method on the Internet. Its principle is similar to that of the angle method, but it has made many important improvements. For example, the calculation corner is changed to the calculation cross product, which is determined by the cross product.
The rotation direction is determined based on the quadrant of the next point. Each deflection is a multiple of 90 degrees.
This algorithm can easily determine whether the point is inside the polygon, on the boundary, or outside the polygon.

This algorithm is described as follows:
First, extract an introduction to the arc-length method from the book: "The Arc-length method requires a polygon to be a directed polygon. It is generally stipulated that the forward side of the polygon is the inner area of the Polygon on the left side of the side.
The center of the measurement point is used as the unit circle, and all directed edges are projected to the unit circle in a radial direction. The algebraic sum of the arc lengths on the unit circle is calculated. If the algebra is 0, the point is outside the polygon;
If the algebra is 2 π, The point is inside the polygon. If the algebra is π, The point is on the polygon ."
According to the introduction in the book, the arc-length method is actually the corner method. However, the method of improvement is quite powerful: The coordinate origin is translated to the point P, and the new coordinate system divides the plane into four
Quadrant. For each polygon vertex P, only the quadrant of the polygon is taken into account. Then, the vertex P of the polygon is accessed in the adjacent order, and P [I + 1] are analyzed. there are three situations:
(1) P [I + 1] is in the next quadrant of P. At this time, the arc length and the addition of π/2;
(2) P [I + 1] is in the upper quadrant of P. In this case, the arc length and π/2 are reduced;
(3) P [I + 1] is in the relative quadrant of Pi. Calculate f = y [I + 1] * x-x [I + 1] * y (Cross Product). If f = 0, the point is on the polygon; if f <0, the arc length and π reduction; if f> 0, the arc length and π addition.
Finally, we can judge the calculated algebra just like the above.
There are two more points to note during implementation. The first one is that if a coordinate of P is 0, it will always be processed as a positive number; the second point is that special processing is required when vertices are tested and polygon vertices are duplicated.

The above is the content in the book, but there is still a problem. That is, when an edge of a polygon is on the coordinate axis and the two vertices are on both sides of the origin, an error occurs.
For example, edge (3, 0)-(-3, 0). According to the above processing, quadrant is the first and second, which makes algebra and addition π/2, it is possible that the final result is measured out of the polygon. In fact
The measured point is on a polygon (this side passes through this point ).
My solution is to calculate the cross product and dot product every time P and P [I + 1] are calculated to determine whether the vertex is on the edge, otherwise, continue the above process.
This sacrifices time, but ensures correctness.
In actual implementation, as you only need to know the Quadrant position of the current point and the previous point, the additional space only needs O (1 ). During implementation, you can change the above "π/2" to 1 and "π" to 2,
In this way, integers can be fully used for calculation. You don't have to consider the order of vertices. You can process them both counterclockwise and clockwise, but the last algebra and symbols are different. Write the entire Algorithm
It is very easy.

The Code is as follows:
# Include <stdio. h>
# Include <math. h>

# Define MAX (100 + 10)
Struct Point
{
Double x, y;
};

Point pts [MAX];
Const int OUT = 0;
Const int IN = 1;
Const int EDGE = 2;
Const double fPre = 1e-8;

Int DblCmp (double fD)
{
If (fabs (fD) <fPre)
{
Return 0;
}
Else
{
Return fD> 0? 1:-1;
}
}

Int GetQuadrant (Point p)
{
Return DblCmp (p. x)> = 0? (DblCmp (p. y)> = 0? 0: 3 ):
(DblCmp (p. y)> = 0? 1: 2 );
}

Double Det (double fX1, double fY1, double fX2, double fY2)
{
Return fX1 * fY2-fX2 * fY1;
}

Int PtInPolygon (Point * pts, int nN, Point p)
{
Int I, j, k;
For (j = 0; j <nN; ++ j)
{
Pts [j]. x-= p. x;
Pts [j]. y-= p. y;
}
Int nA1, nA2;
Int nSum = 0;
NA1 = GetQuadrant (pts [0]);
For (I = 0; I <nN; ++ I)
{
K = (I + 1) % nN;
If (DblCmp (pts [k]. x) = 0 & DblCmp (pts [k]. y) = 0)
{
Break; // coincides with the vertex
}
Int nC = DblCmp (Det (pts [I]. x, pts [I]. y,
Pts [k]. x, pts [k]. y ));
If (! NC & DblCmp (pts [I]. x * pts [k]. x) <= 0
& DblCmp (pts [I]. y * pts [k]. y) <= 0)
{
Break; // edge
}
NA2 = GetQuadrant (pts [k]);
If (nA1 + 1) % 4 = nA2)
{
NSum + = 1;
}
Else if (nA1 + 2) % 4 = nA2)
{
If (nC> 0)
{
NSum + = 2;
}
Else
{
NSum-= 2;
}
}
Else if (nA1 + 3) % 4 = nA2)
{
NSum-= 1;
}
NA1 = nA2;
}

For (j = 0; j <nN; ++ j)
{
Pts [j]. x + = p. x;
Pts [j]. y + = p. y;
}

If (I <nN)
{
Return EDGE;
}
Else if (nSum) // counterclockwise nSum = 4, clockwise nSum =-4
{
Return IN;
}
Else
{
Return OUT;
}
}

Int main ()
{
Int nN, nM;
Int nCase = 1;

While (scanf ("% d", & nN, & nM), nN)
{
If (nCase> 1)
{
Printf ("\ n ");
}

For (int I = 0; I <nN; ++ I)
{
Scanf ("% lf", & pts [I]. x, & pts [I]. y );
}
Printf ("Problem % d: \ n", nCase ++ );
For (int I = 0; I <nM; ++ I)
{
Point p;
Scanf ("% lf", & p. x, & p. y );
If (PtInPolygon (pts, nN, p ))
{
Printf ("Within \ n ");
}
Else
{
Printf ("Outside \ n ");
}
}
}

Return 0;
}