Polygons: Area, Winding and Point-in-Polygon
Exact polygon area with the shoelace formula, a convexity test that rejects stars, and point-in-polygon that stays right when the ray hits a vertex.
Three questions every polygon gets asked
A land registry stores each parcel as a list of corners and needs its area. A map application has to decide whether the point a user tapped lies inside a delivery zone. A physics engine accepts only convex shapes and has to reject the others before they reach the solver. Each of these is a few lines of code, and each has a version that is almost right.
The almost-right versions fail in recognizable ways. The area comes out negative because the corners were listed clockwise. The tap test says "outside" for a point that is plainly inside, because the horizontal line through the point passes through a corner of the polygon. A five-pointed star passes the convexity test because every corner turns left. A small grid of integer coordinates is enough to hit all three.
This lesson builds the three answers on the exact orientation test from lesson 1, each with the case that breaks the naive version, and checks point-in-polygon against an independent reference on a quarter of a million queries.
What counts as a polygon in this course
A polygon here is a std::vector<Point> of its corners (its vertices), in order. Consecutive vertices are joined by edges, and the last vertex is joined back to the first. Three conventions come with that sentence, and most polygon bugs come from code that silently assumes a different one:
- No closing vertex. The first vertex is not repeated at the end. Some file formats repeat it (GeoJSON requires it), so strip it when you read such data.
- Counter-clockwise order. Walking along the edges, the interior is on your left. The next section shows how to detect and fix the order with one function call.
- Integer coordinates.
struct Point { long long x, y; };as in the rest of the course, so that every predicate is exact.
A polygon is simple if its edges meet only where consecutive edges share a vertex: no edge crosses another, and no vertex touches an edge it does not belong to. A bow tie, whose two halves meet in the middle, is not simple. Neither is a star drawn with five strokes. Most algorithms in this and the next lesson assume a simple polygon. Point-in-polygon is the exception, and it is where the difference between the two kinds becomes visible.
Area from a fan of triangles: the shoelace formula
Lesson 1 showed that cross(a, b, c) is twice the signed area of triangle abc, positive when the corners run counter-clockwise. A polygon's area can be assembled from such triangles.
Pick any point o, for convenience the origin. For each edge from p to q, form the triangle o, p, q. Walk around a counter-clockwise polygon and add up the signed areas of these triangles. Where an edge is on the far side of the polygon as seen from o, its triangle runs counter-clockwise and counts positive. Where an edge is on the near side, its triangle runs clockwise and counts negative. The triangles overlap, but every point outside the polygon is covered by as many positive as negative triangles, so they cancel, and what remains is the area of the polygon. This works for every simple polygon, convex or not, and for any choice of o.
With o at the origin, cross(o, p, q) simplifies to p.x * q.y - q.x * p.y, and the whole formula is
2A = sum over edges (p, q) of (p.x * q.y - q.x * p.y).
The name "shoelace" comes from the crossing pattern you get when you write the coordinates in two columns and multiply diagonally.
Do one by hand. The L-shaped polygon (0,0), (4,0), (4,1), (1,1), (1,3), (0,3) has edges whose terms are 0, 4, 3, 2, 3, 0. The sum is 12, so the area is 6: a 4 by 1 strip plus a 1 by 2 strip. The program below computes the same, then shows what happens with clockwise order, a translated copy and a bow tie.
#include <algorithm>
#include <limits>
#include <print>
#include <vector>
struct Point { long long x, y; };
// The bounds this lesson promises: shoelace terms are at most 2 * kMaxCoord^2
// in absolute value, and the sum of kMaxVertices of them must still fit.
constexpr long long kMaxCoord = 1'000'000;
constexpr long long kMaxVertices = 1'000'000;
static_assert(2 * kMaxCoord * kMaxCoord * kMaxVertices
<= std::numeric_limits<long long>::max(),
"shoelace sum must not overflow");
// Twice the signed area: positive for counter-clockwise, negative for clockwise.
long long twiceArea(const std::vector<Point>& poly) {
long long sum = 0;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point p = poly[i], q = poly[(i + 1) % n];
sum += p.x * q.y - q.x * p.y; // cross(origin, p, q)
}
return sum;
}
void makeCounterClockwise(std::vector<Point>& poly) {
if (twiceArea(poly) < 0) std::ranges::reverse(poly);
}
int main() {
std::vector<Point> rect{{0, 0}, {4, 0}, {4, 3}, {0, 3}};
std::vector<Point> rectCw{{0, 0}, {0, 3}, {4, 3}, {4, 0}};
std::vector<Point> lShape{{0, 0}, {4, 0}, {4, 1}, {1, 1}, {1, 3}, {0, 3}};
std::vector<Point> moved = lShape;
for (Point& p : moved) p = {p.x + 1000, p.y - 500}; // same shape, elsewhere
std::vector<Point> bowTie{{0, 0}, {2, 2}, {2, 0}, {0, 2}};
std::println("rectangle, CCW: {}", twiceArea(rect));
std::println("rectangle, CW: {}", twiceArea(rectCw));
std::println("L-shape: {}", twiceArea(lShape));
std::println("L-shape, moved: {}", twiceArea(moved));
std::println("bow tie: {}", twiceArea(bowTie));
makeCounterClockwise(rectCw);
std::println("rectangle after makeCounterClockwise: {}", twiceArea(rectCw));
}Press Run to compile and run this example.
rectangle, CCW: 24
rectangle, CW: -24
L-shape: 12
L-shape, moved: 12
bow tie: 0
rectangle after makeCounterClockwise: 24The sign is the orientation: a clockwise rectangle has area -24, so if twiceArea is negative, reverse the vertices. Moving the polygon by (1000, -500) changes every term but not the sum, as the fan argument promises; the terms grow with the distance from the origin, which is why the overflow bound below matters.
The bow tie has area 0. Its left half runs counter-clockwise and its right half clockwise, and the formula adds them with opposite signs. For a polygon that crosses itself, the shoelace sum is not "the area" in any everyday sense. It counts each region as many times as the boundary winds around it, with a sign. That quantity, the winding number, comes back below.
The overflow bound. Each term p.x * q.y - q.x * p.y is at most 2 * M^2 in absolute value when every coordinate satisfies |x|, |y| <= M. With M = 10^9, the bound from lesson 1, one term can reach 2 * 10^18, and the running sum of a few terms can leave the range of long long (about 9.22 * 10^18) even if the final result would fit. Rather than reason about partial sums, this lesson promises a smaller box: |x|, |y| <= 10^6 and at most 10^6 vertices. Then the whole sum is at most 2 * 10^18, and the static_assert checks that arithmetic at compile time. If you need the larger coordinates, accumulate in a wider type such as __int128 (a GCC and Clang extension, not standard C++); lesson 2 discusses the trade-off.
Convex or not: left turns are not enough
A polygon is convex if the segment between any two of its points stays inside it. For a counter-clockwise polygon, the obvious test is that every corner turns left: cross(p[i], p[i+1], p[i+2]) > 0 for every i.
That test is necessary, and it is not sufficient. A five-pointed star drawn with five strokes, visiting every second corner of a pentagon, turns left at every vertex. Its edges go round the center twice instead of once, so it is not even simple. The left-turn test cannot see this, because it looks at one corner at a time.
The fix is to also count how many times the boundary goes round. As you walk along a convex polygon, the direction of travel rotates through exactly one full turn. During that turn, the x-component of the direction changes sign exactly twice: once at the rightmost point, where you stop heading right and start heading left, and once at the leftmost point. A star that goes round twice changes sign four times. Counting those switches is exact, because it only compares integer coordinates.
Two degenerate cases need a decision:
- Collinear vertices, where three consecutive vertices lie on a line and the turn is
0. A square with an extra vertex in the middle of one side is still a convex region. Some callers want it accepted, others (a physics engine computing face normals, for example) need every vertex to be a real corner. Make it a parameter. - Reversals, where the turn is
0because the boundary doubles back on itself, as in the zero-area polygon(0,0), (4,0), (2,0). A dot product tells the two apart: when the boundary keeps going forward,(b - a) . (c - b)is positive.
#include <print>
#include <vector>
struct Point { long long x, y; };
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
long long dot(Point a, Point b, Point c) { // (b - a) . (c - b): still going forward?
return (b.x - a.x) * (c.x - b.x) + (b.y - a.y) * (c.y - b.y);
}
enum class Collinear { Allow, Reject };
// Only the turns: the test most people write first.
bool allLeftTurns(const std::vector<Point>& poly) {
for (std::size_t i = 0, n = poly.size(); i < n; ++i)
if (cross(poly[i], poly[(i + 1) % n], poly[(i + 2) % n]) < 0) return false;
return true;
}
// Convex and counter-clockwise: every turn is a left turn (or straight on, if
// allowed), and the edge directions go round exactly once.
bool isConvex(const std::vector<Point>& poly, Collinear policy) {
const std::size_t n = poly.size();
if (n < 3) return false;
std::vector<int> dxSigns; // sign of dx for every edge that is not vertical
for (std::size_t i = 0; i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n], c = poly[(i + 2) % n];
long long turn = cross(a, b, c);
if (turn < 0) return false; // a right turn
if (turn == 0 && dot(a, b, c) <= 0) return false; // spike or repeated vertex
if (turn == 0 && policy == Collinear::Reject) return false;
if (b.x != a.x) dxSigns.push_back(b.x > a.x ? 1 : -1);
}
// Going round once, the edges switch between heading right and heading left
// exactly twice. A star that goes round k times switches 2k times.
int switches = 0;
for (std::size_t i = 0; i < dxSigns.size(); ++i)
switches += dxSigns[i] != dxSigns[(i + 1) % dxSigns.size()];
return switches == 2;
}
int main() {
struct Case { const char* name; std::vector<Point> poly; };
const Case cases[] = {
{"square", {{0, 0}, {4, 0}, {4, 4}, {0, 4}}},
{"square, clockwise", {{0, 0}, {0, 4}, {4, 4}, {4, 0}}},
{"square + midpoint", {{0, 0}, {2, 0}, {4, 0}, {4, 4}, {0, 4}}},
{"arrowhead", {{0, 0}, {4, 2}, {0, 4}, {1, 2}}},
{"flat (zero area)", {{0, 0}, {4, 0}, {2, 0}}},
{"pentagram", {{5, 10}, {2, 0}, {10, 6}, {0, 6}, {8, 0}}},
};
std::println("{:<18} {:<11} {:<16} {}", "", "left turns", "allow collinear",
"strict");
for (const auto& [name, poly] : cases) {
std::println("{:<18} {:<11} {:<16} {}", name, allLeftTurns(poly),
isConvex(poly, Collinear::Allow), isConvex(poly, Collinear::Reject));
}
}Press Run to compile and run this example.
left turns allow collinear strict
square true true true
square, clockwise false false false
square + midpoint true true false
arrowhead false false false
flat (zero area) true false false
pentagram true false falseThe left-turn test accepts two shapes that are not convex polygons: the flat one and the pentagram. The full test rejects both. Note also that the clockwise square is rejected: isConvex checks "convex and counter-clockwise", which is what code downstream relies on. Normalize the order first if the input may be clockwise.
The whole test is one pass, O(n). It was checked in the background against a brute-force definition (every vertex on or to the left of every edge's line, distinct vertices, positive area) on two million random polygons with 3 to 7 vertices on a 4 by 4 grid, under both collinear policies, with no disagreement.
Is the point inside? Counting crossings with a half-open rule
The idea of the classic test fits in one sentence: shoot a ray from the query point to the right and count how many edges it crosses. An odd count means inside, an even count means outside, because every crossing switches between the two. This is the crossing number or ray casting test.
The sentence hides two decisions.
The first is what to do when the ray passes exactly through a vertex. The vertex belongs to two edges, so a naive count sees two crossings where the boundary crossed once, or crosses where the boundary only touched the ray and turned back. The standard fix is the half-open rule: treat every vertex that lies on the ray as if it were slightly below it. In code, an edge counts only if exactly one of its endpoints is strictly above the query point, (a.y > p.y) != (b.y > p.y). A vertex where the boundary passes through the ray then counts once, a vertex where it touches and turns back counts zero or two times, and horizontal edges never count. No division is needed: to decide whether the crossing lies to the right of p, ask on which side of the edge p is.
The second decision is what to say for points exactly on the boundary. The parity count gives an arbitrary answer there (W. Randolph Franklin's widely copied PNPOLY function documents that boundary points come out as "either inside or outside"). Some applications are happy with that: if a map is split into polygons that share edges, a half-open rule assigns each point to exactly one of them. Others, like the land registry, need to know. So the function returns three values, Inside, Outside and Boundary, and checks the boundary first with the exact on-segment test from lesson 1.
#include <algorithm>
#include <print>
#include <vector>
struct Point { long long x, y; };
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool onSegment(Point a, Point b, Point p) {
return cross(a, b, p) == 0 && std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x)
&& std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}
enum class Where { Outside, Inside, Boundary };
// Crossing number: shoot a ray from p towards +x and count the edges it crosses.
Where crossingTest(const std::vector<Point>& poly, Point p) {
bool inside = false;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
if (onSegment(a, b, p)) return Where::Boundary;
// Half-open rule: the edge counts only if exactly one endpoint is
// strictly above p. A vertex on the ray belongs to the edge going up from it.
if ((a.y > p.y) != (b.y > p.y)) {
// The crossing is to the right of p if p is left of the upward edge.
int side = orientation(a, b, p);
if (b.y > a.y ? side > 0 : side < 0) inside = !inside;
}
}
return inside ? Where::Inside : Where::Outside;
}
const char* name(Where w) {
switch (w) {
case Where::Outside: return "outside";
case Where::Inside: return "inside";
case Where::Boundary: return "boundary";
}
return "?";
}
int main() {
// Counter-clockwise; the line y = 2 runs through three vertices and an edge.
std::vector<Point> poly{{0, 0}, {8, 0}, {8, 2}, {6, 4}, {4, 2}, {2, 2}, {0, 4}};
for (Point p : {Point{1, 2}, Point{5, 2}, Point{7, 2}, Point{9, 2},
Point{3, 2}, Point{4, 1}, Point{6, 4}, Point{3, 3}}) {
std::println("({}, {}) is {}", p.x, p.y, name(crossingTest(poly, p)));
}
}Press Run to compile and run this example.
(1, 2) is inside
(5, 2) is inside
(7, 2) is inside
(9, 2) is outside
(3, 2) is boundary
(4, 1) is inside
(6, 4) is boundary
(3, 3) is outsideThe polygon was chosen to be hostile. The horizontal line y = 2 contains the vertices (8, 2), (4, 2) and (2, 2) and the whole edge between the last two. Follow the ray from (1, 2) by hand. The edge from (8, 0) to (8, 2) has no endpoint strictly above y = 2, so it does not count. The edge from (8, 2) up to (6, 4) counts, and so does the edge from (6, 4) down to (4, 2): the boundary arrives at (4, 2) from above and leaves along the horizontal edge, which never counts. The edge from (2, 2) up to (0, 4) counts, because its crossing at x = 2 is to the right of x = 1. Three crossings: inside.
The test is O(n) per query with one orientation per straddling edge.
Here is the same polygon with every grid point from (-1, -1) to (9, 5) classified: blue inside, green on the boundary, small dark dots outside. The program draws in the usual y-up coordinates and flips only when writing the SVG.
#include <algorithm>
#include <format>
#include <print>
#include <string>
#include <vector>
struct Point { long long x, y; };
class Svg {
public:
Svg(double width, double height) : w_(width), h_(height) {}
void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
a.x, fy(a.y), b.x, fy(b.y), color, width);
}
void dot(Point p, const char* color = "#1f2937", double r = 3) {
add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
}
void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
const char* fill = "none") {
std::string list;
for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
}
void text(Point p, const std::string& s, double size = 8) {
add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
}
void print() const {
std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
std::print("{}", body_);
std::println("</svg>");
}
private:
double fy(double y) const { return h_ - y; } // SVG's y axis points down
template <class... Args>
void add(std::format_string<Args...> fmt, Args&&... args) {
body_ += std::format(fmt, std::forward<Args>(args)...);
body_ += '\n';
}
double w_, h_;
std::string body_;
};
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool onSegment(Point a, Point b, Point p) {
return cross(a, b, p) == 0 && std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x)
&& std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}
enum class Where { Outside, Inside, Boundary };
Where crossingTest(const std::vector<Point>& poly, Point p) {
bool inside = false;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
if (onSegment(a, b, p)) return Where::Boundary;
if ((a.y > p.y) != (b.y > p.y)) {
int side = orientation(a, b, p);
if (b.y > a.y ? side > 0 : side < 0) inside = !inside;
}
}
return inside ? Where::Inside : Where::Outside;
}
Point toSvg(Point p) { return {20 * p.x + 30, 20 * p.y + 30}; } // grid -> picture
int main() {
std::vector<Point> poly{{0, 0}, {8, 0}, {8, 2}, {6, 4}, {4, 2}, {2, 2}, {0, 4}};
Svg svg(220, 140);
std::vector<Point> drawn;
for (Point p : poly) drawn.push_back(toSvg(p));
svg.polygon(drawn, "#1f2937", "rgba(37,99,235,0.15)");
for (long long x = -1; x <= 9; ++x) {
for (long long y = -1; y <= 5; ++y) {
Where w = crossingTest(poly, {x, y});
const char* color = w == Where::Inside ? "#2563eb"
: w == Where::Boundary ? "#16a34a"
: "#1f2937";
svg.dot(toSvg({x, y}), color, w == Where::Outside ? 1.5 : 3);
}
}
svg.print();
}Press Run to compile and run this example.
Every green dot sits on an edge, including the ones on the horizontal edge and at the vertices of the valley, and no dot on the line y = 2 is misclassified.
Winding numbers and self-intersecting polygons
The crossing test answers "is the count odd?". The winding number answers a different question: how many times does the boundary go round the point counter-clockwise? It uses the same rays and the same half-open rule, but an edge crossing the ray upward adds one, and an edge crossing it downward subtracts one.
For a simple polygon the two tests agree: the winding number is 1 inside a counter-clockwise polygon, -1 inside a clockwise one and 0 outside, and the crossing count is odd exactly when the winding number is not zero. For a polygon that crosses itself they differ, and the difference has names. The even-odd rule calls a point inside when the crossing count is odd. The non-zero rule calls it inside when the winding number is not zero. SVG's fill-rule property offers both, with nonzero as the default.
#include <algorithm>
#include <print>
#include <vector>
struct Point { long long x, y; };
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
// Winding number of the closed polygon around p (p must not be on the boundary).
// Same half-open rule as the crossing test, but upward crossings add one and
// downward crossings subtract one instead of both toggling a flag.
int windingNumber(const std::vector<Point>& poly, Point p) {
int winding = 0;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
if (a.y <= p.y && b.y > p.y && orientation(a, b, p) > 0) ++winding;
if (a.y > p.y && b.y <= p.y && orientation(a, b, p) < 0) --winding;
}
return winding;
}
int crossings(const std::vector<Point>& poly, Point p) {
int count = 0;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
if ((a.y > p.y) != (b.y > p.y)) {
int side = orientation(a, b, p);
count += b.y > a.y ? side > 0 : side < 0;
}
}
return count;
}
int main() {
std::vector<Point> star{{5, 10}, {2, 0}, {10, 6}, {0, 6}, {8, 0}};
std::vector<Point> square{{0, 0}, {4, 0}, {4, 4}, {0, 4}};
std::vector<Point> squareCw{{0, 0}, {0, 4}, {4, 4}, {4, 0}};
auto report = [](const char* what, const std::vector<Point>& poly, Point p) {
int c = crossings(poly, p), w = windingNumber(poly, p);
std::println("{:<22} crossings {} even-odd: {:<7} winding {:>2} non-zero: {}",
what, c, c % 2 ? "inside" : "outside", w, w ? "inside" : "outside");
};
report("star, center (5, 5)", star, {5, 5});
report("star, top tip (5, 8)", star, {5, 8});
report("star, outside (1, 1)", star, {1, 1});
report("square, CCW (2, 2)", square, {2, 2});
report("square, CW (2, 2)", squareCw, {2, 2});
}Press Run to compile and run this example.
star, center (5, 5) crossings 2 even-odd: outside winding 2 non-zero: inside
star, top tip (5, 8) crossings 1 even-odd: inside winding 1 non-zero: inside
star, outside (1, 1) crossings 4 even-odd: outside winding 0 non-zero: outside
square, CCW (2, 2) crossings 1 even-odd: inside winding 1 non-zero: inside
square, CW (2, 2) crossings 1 even-odd: inside winding -1 non-zero: insideThe pentagram's center is wound round twice. The even-odd rule cuts it out, which is why a star drawn with evenodd has a hole in the middle, and the non-zero rule fills it. The point (1, 1) is outside by both rules, even though its ray crosses four edges. And the winding number carries the orientation: the clockwise square gives -1. That sign is how the shoelace formula treated the bow tie: the signed area it computes (half its sum) is the integral of the winding number over the plane, so a region wound +1 and a region wound -1 cancel.
Which test should you use? For simple polygons either one, and the crossing test is a little shorter. For polygons that may cross themselves, decide which rule your data means. Fonts and vector graphics usually specify one explicitly, and your test must match it.
Testing point-in-polygon against an independent reference
The examples above check eight points by hand. To trust the code, test it on thousands of polygons and every point of a grid around them, against a second implementation that shares no tie-breaking logic.
The reference uses a trick that removes the ties instead of breaking them. For grid points with coordinates in [0, R], aim the ray not horizontally but in the direction (2R + 1, 1). To pass through another grid vertex, the ray would have to rise by a whole unit, and by then it has moved more than R to the right, beyond every vertex. An edge cannot lie along the ray for the same reason. So every crossing is a proper crossing in the middle of an edge, and the four-orientation test of two segments decides it without any special cases.
The random polygons come from sorting random grid points by angle around the lowest one, with ties on the last ray reversed so the boundary does not overlap itself. Such polygons are star-shaped from that point, so they are not the most twisted simple polygons there are. On a 7 by 7 grid they have plenty of what matters here: horizontal edges, collinear vertices and vertices exactly on the query ray.
#include <algorithm>
#include <print>
#include <vector>
struct Point { long long x, y; };
struct Lcg {
unsigned long long s;
unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); } // inclusive
};
long long cross(Point a, Point b, Point c) {
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
long long v = cross(a, b, c);
return (v > 0) - (v < 0);
}
bool onSegment(Point a, Point b, Point p) {
return cross(a, b, p) == 0 && std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x)
&& std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}
enum class Where { Outside, Inside, Boundary };
Where crossingTest(const std::vector<Point>& poly, Point p) {
bool inside = false;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
if (onSegment(a, b, p)) return Where::Boundary;
if ((a.y > p.y) != (b.y > p.y)) {
int side = orientation(a, b, p);
if (b.y > a.y ? side > 0 : side < 0) inside = !inside;
}
}
return inside ? Where::Inside : Where::Outside;
}
Where windingTest(const std::vector<Point>& poly, Point p) {
int winding = 0;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
if (onSegment(a, b, p)) return Where::Boundary;
if (a.y <= p.y && b.y > p.y && orientation(a, b, p) > 0) ++winding;
if (a.y > p.y && b.y <= p.y && orientation(a, b, p) < 0) --winding;
}
return winding != 0 ? Where::Inside : Where::Outside;
}
// Independent reference for grid points in [0, R]^2. A ray with direction (2R+1, 1)
// cannot pass through a grid vertex or run along an edge, so every crossing is a
// proper crossing and no tie-breaking rule is needed.
Where reference(const std::vector<Point>& poly, Point p, long long R) {
Point far{p.x + (R + 1) * (2 * R + 1), p.y + (R + 1)}; // beyond the grid
int count = 0;
for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
Point a = poly[i], b = poly[(i + 1) % n];
long long dot = (a.x - p.x) * (b.x - p.x) + (a.y - p.y) * (b.y - p.y);
if (cross(a, b, p) == 0 && dot <= 0) return Where::Boundary;
count += orientation(a, b, p) * orientation(a, b, far) < 0
&& orientation(p, far, a) * orientation(p, far, b) < 0;
}
return count % 2 ? Where::Inside : Where::Outside;
}
// A random simple polygon: sort points by angle around the lowest one. Ties on a
// ray go outwards, except on the last ray, which must come back inwards.
std::vector<Point> randomPolygon(Lcg& rng, int n, int R) {
std::vector<Point> pts;
for (int i = 0; i < n; ++i) pts.push_back({rng.range(0, R), rng.range(0, R)});
auto byYx = [](Point a, Point b) { return a.y != b.y ? a.y < b.y : a.x < b.x; };
std::ranges::sort(pts, byYx);
auto same = [](Point a, Point b) { return a.x == b.x && a.y == b.y; };
pts.erase(std::unique(pts.begin(), pts.end(), same), pts.end());
Point o = pts[0];
auto dist = [o](Point p) {
return (p.x - o.x) * (p.x - o.x) + (p.y - o.y) * (p.y - o.y);
};
std::sort(pts.begin() + 1, pts.end(), [&](Point a, Point b) {
int turn = orientation(o, a, b);
return turn != 0 ? turn > 0 : dist(a) < dist(b);
});
auto last = pts.end() - 1;
while (last - 1 > pts.begin() && orientation(o, *(last - 1), pts.back()) == 0) --last;
std::reverse(last, pts.end());
return pts;
}
int main() {
Lcg rng{2024};
const int R = 6;
long polygons = 0, queries = 0, crossingWrong = 0, windingWrong = 0;
long tally[3] = {};
while (polygons < 3000) {
auto poly = randomPolygon(rng, rng.range(3, 12), R);
long long area = 0;
for (std::size_t i = 0; i < poly.size(); ++i)
area += cross({0, 0}, poly[i], poly[(i + 1) % poly.size()]);
if (area <= 0) continue; // all points collinear: not a polygon
++polygons;
for (long long x = -1; x <= R + 1; ++x)
for (long long y = -1; y <= R + 1; ++y) {
Where want = reference(poly, {x, y}, R + 2);
crossingWrong += crossingTest(poly, {x, y}) != want;
windingWrong += windingTest(poly, {x, y}) != want;
++tally[int(want)];
++queries;
}
}
std::println("polygons {}, queries {}", polygons, queries);
std::println("outside {}, inside {}, boundary {}", tally[0], tally[1], tally[2]);
std::println("crossing test disagrees with the reference: {}", crossingWrong);
std::println("winding test disagrees with the reference: {}", windingWrong);
}Press Run to compile and run this example.
polygons 3000, queries 243000
outside 188354, inside 22714, boundary 31932
crossing test disagrees with the reference: 0
winding test disagrees with the reference: 0All three answers occur many times, so a bug in any branch would have shown up. As a check of the check, replacing the half-open condition with the closed one (min(a.y, b.y) <= p.y <= max(a.y, b.y)) on a similar run produced about 143,000 disagreements out of 1.6 million queries.