Skip to content

Commit 7681b3e

Browse files
Chilleesimonlindholm
authored andcommitted
Updated inPolygon (#81)
1 parent 9b21ef1 commit 7681b3e

File tree

9 files changed

+175
-43
lines changed

9 files changed

+175
-43
lines changed

content/geometry/InsidePolygon.h

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
/**
2+
* Author: Victor Lecomte, chilli
3+
* Date: 2019-04-26
4+
* License: CC0
5+
* Source: https://vlecomte.github.io/cp-geo.pdf
6+
* Description: Returns true if p lies within the polygon. If strict is true,
7+
* it returns false for points on the boundary. The algorithm uses
8+
* products in intermediate steps so watch out for overflow.
9+
* Time: O(n)
10+
* Status: fuzz-tested and Kattis problem pointinpolygon
11+
* Usage:
12+
* vector<P> v = {P{4,4}, P{1,2}, P{2,1}};
13+
* bool in = inPolygon(v, P{3, 3}, false);
14+
*/
15+
#pragma once
16+
17+
#include "Point.h"
18+
#include "onSegment.h"
19+
#include "SegmentDistance.h"
20+
21+
template<class P>
22+
bool inPolygon(vector<P> &p, P a, bool strict = true) {
23+
int cnt = 0, n = sz(p);
24+
rep(i,0,n) {
25+
P q = p[(i + 1) % n];
26+
if (onSegment(p[i], q, a)) return !strict;
27+
//or: if (segDist(p[i], q, a) <= eps) return !strict;
28+
cnt ^= ((a.y<p[i].y) - (a.y<q.y)) * a.cross(p[i], q) > 0;
29+
}
30+
return cnt;
31+
}

content/geometry/chapter.tex

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ \section{Circles}
1919
\kactlimport{MinimumEnclosingCircle.h}
2020

2121
\section{Polygons}
22-
\kactlimport{insidePolygon.h}
22+
\kactlimport{InsidePolygon.h}
2323
\kactlimport{PolygonArea.h}
2424
\kactlimport{PolygonCenter.h}
2525
\kactlimport{PolygonCut.h}

content/geometry/insidePolygon.h

Lines changed: 0 additions & 34 deletions
This file was deleted.

content/geometry/onSegment.h

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,8 @@
33
* Date: 2019-04-26
44
* License: CC0
55
* Source: https://vlecomte.github.io/cp-geo.pdf
6-
* Description: Returns true iff p lies on the line segment from s to e. Intended for use with e.g. Point<ll> where overflow is an issue. Use \texttt{(segDist(s,e,p)<=epsilon)} instead when using Point<double>.
6+
* Description: Returns true iff p lies on the line segment from s to e.
7+
* Use \texttt{(segDist(s,e,p)<=epsilon)} instead when using Point<double>.
78
* Status:
89
*/
910
#pragma once

fuzz-tests/geometry/PointInsideHull.cpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ typedef vector<int> vi;
1212

1313
#include "../../content/geometry/ConvexHull.h"
1414
#include "../../content/geometry/PointInsideHull.h"
15-
#include "../../content/geometry/insidePolygon.h"
15+
#include "../../content/geometry/InsidePolygon.h"
1616

1717
int main() {
1818
rep(it,0,100000) {
@@ -25,8 +25,8 @@ int main() {
2525
int x = rand() % 22 - 11;
2626
int y = rand() % 22 - 11;
2727
P p{x,y};
28-
assert(insidePolygon(all(ps2), p, true) == (insideHull(ps2, p) == 2));
29-
assert(insidePolygon(all(ps2), p, false) == (insideHull(ps2, p) >= 1));
28+
assert(inPolygon(all(ps2), p, true) == (insideHull(ps2, p) == 2));
29+
assert(inPolygon(all(ps2), p, false) == (insideHull(ps2, p) >= 1));
3030
}
3131
}
3232
}

fuzz-tests/geometry/PolygonArea.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ typedef vector<int> vi;
1212

1313
#include "../../content/geometry/PolygonArea.h"
1414
#include "../../content/geometry/PolygonCenter.h"
15-
#include "../../content/geometry/insidePolygon.h"
15+
#include "../../content/geometry/InsidePolygon.h"
1616

1717
int main() {
1818
srand((int)time(0));
@@ -25,7 +25,7 @@ int main() {
2525
double y = rand() / (RAND_MAX + 1.0);
2626
x *= 10;
2727
y *= 10;
28-
if (!insidePolygon(all(ps), P{x,y}, true)) continue;
28+
if (!inPolygon(ps, P{x,y}, true)) continue;
2929
count++;
3030
su = su + P{x,y};
3131
}

fuzz-tests/geometry/PolygonCut.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ typedef vector<int> vi;
1313
#include "../../content/geometry/PolygonArea.h"
1414
#include "../../content/geometry/PolygonCut.h"
1515
#include "../../content/geometry/sideOf.h"
16-
#include "../../content/geometry/insidePolygon.h"
16+
#include "../../content/geometry/InsidePolygon.h"
1717
#include "../../content/geometry/SegmentIntersection.h"
1818

1919
typedef Point<double> P;
@@ -41,7 +41,7 @@ int main() {
4141
rep(it,0,ITS) {
4242
double x = rand() / (RAND_MAX + 1.0) * 10 - 5;
4343
double y = rand() / (RAND_MAX + 1.0) * 10 - 5;
44-
if (!insidePolygon(all(ps), P{x,y}, true)) continue;
44+
if (!inPolygon(ps, P{x,y}, true)) continue;
4545
if (sideOf(p, q, P{x,y}) > 0) continue;
4646
count++;
4747
}
Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
#include <bits/stdc++.h>
2+
using namespace std;
3+
4+
#define rep(i, a, b) for(int i = a; i < int(b); ++i)
5+
#define trav(a, v) for(auto& a : v)
6+
#define all(x) x.begin(), x.end()
7+
#define sz(x) (int)(x).size()
8+
9+
typedef long long ll;
10+
typedef pair<int, int> pii;
11+
typedef vector<int> vi;
12+
13+
14+
const double EPS =1e-8;
15+
#include "../utilities/genPolygon.h"
16+
#include "../../content/geometry/InsidePolygon.h"
17+
namespace old {
18+
19+
#include "../../content/geometry/onSegment.h"
20+
#include "../../content/geometry/SegmentDistance.h"
21+
22+
template<class It, class P>
23+
bool insidePolygon(It begin, It end, const P& p,
24+
bool strict = true) {
25+
int n = 0; //number of isects with line from p to (inf,p.y)
26+
for (It i = begin, j = end-1; i != end; j = i++) {
27+
//if p is on edge of polygon
28+
if (onSegment(*i, *j, p)) return !strict;
29+
//or: if (segDist(*i, *j, p) <= epsilon) return !strict;
30+
//increment n if segment intersects line from p
31+
n += (max(i->y,j->y) > p.y && min(i->y,j->y) <= p.y &&
32+
((*j-*i).cross(p-*i) > 0) == (i->y <= p.y));
33+
}
34+
return n&1; //inside if odd number of intersections
35+
}
36+
}
37+
typedef Point<double> P;
38+
ostream &operator<<(ostream &os, P p) { return cout << "(" << p.x << "," << p.y << ")"; }
39+
bool eq(P a, P b) {
40+
return (a-b).dist()<EPS;
41+
}
42+
const int NUMPOLY=100;
43+
const int PTPERPOLY=100;
44+
void test(int numPts, int range) {
45+
rep(i,0,NUMPOLY) {
46+
vector<P> poly;
47+
rep(j,0, numPts)
48+
poly.push_back(P(rand()%range, rand()%range));
49+
poly = genPolygon(poly);
50+
rep(i,0,PTPERPOLY){
51+
P p(rand()%range, rand()%range);
52+
assert(inPolygon(poly, p, true) == old::insidePolygon(all(poly), p, true));
53+
assert(inPolygon(poly, p, false) == old::insidePolygon(all(poly), p, false));
54+
}
55+
}
56+
57+
}
58+
int main() {
59+
test(20,5);
60+
test(1001,100);
61+
test(1000,1000);
62+
cout<<"Tests passed!"<<endl;
63+
}

fuzz-tests/utilities/genPolygon.h

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
#pragma once
2+
3+
// Translated from Python code posted here: https://codeforces.com/blog/entry/63058?#comment-472788
4+
// May generate polygons with colinear points. Won't generate polygons with duplicate points.
5+
#include "../../content/geometry/Point.h"
6+
7+
mt19937_64 rng(time(0));
8+
template<class P> pair<bool, vector<P>> conquer(vector<P> pts, int depth) {
9+
if (depth>100) {
10+
return {false, {}};
11+
}
12+
if (sz(pts) <= 2) return {true,pts};
13+
if (sz(pts) == 3) {
14+
swap(pts[1], pts[2]);
15+
return {true,pts};
16+
}
17+
18+
uniform_int_distribution<int> uni(2, sz(pts)-1);
19+
int divideId = uni(rng);
20+
P p1 = pts[divideId];
21+
uniform_real_distribution<> uniFloat(0.01, 0.99);
22+
double divideK = uniFloat(rng);
23+
P p2(divideK*(pts[1].x-pts[0].x) + pts[0].x, divideK*(pts[1].y - pts[0].y) + pts[0].y);
24+
vector<double> line = {p2.y - p1.y, p1.x - p2.x, -p1.x*p2.y + p1.y*p2.x};
25+
int idx0 = ((line[0]*pts[0].x + line[1]*pts[0].y + line[2]) >=0);
26+
int idx1 = ((line[0]*pts[1].x + line[1]*pts[1].y + line[2]) >=0);
27+
if (idx0==idx1)
28+
return conquer(pts, depth+1);
29+
array<vector<P>, 2> S;
30+
S[idx0].push_back(pts[0]);
31+
S[idx0].push_back(p1);
32+
S[!idx0].push_back(p1);
33+
S[!idx0].push_back(pts[1]);
34+
rep(i,2,sz(pts)) {
35+
if (i == divideId) continue;
36+
int idx = ((line[0]*pts[i].x + line[1]*pts[i].y + line[2]) >=0);
37+
S[idx].push_back(pts[i]);
38+
}
39+
auto pa = conquer(S[idx0], depth+1);
40+
auto pb = conquer(S[!idx0], depth+1);
41+
if (!pa.first || !pb.first) return {false, {}};
42+
pb.second.erase(pb.second.begin());
43+
pa.second.insert(pa.second.end(), all(pb.second));
44+
return pa;
45+
}
46+
template<class P> vector<P> genPolygon(vector<P> pts) {
47+
sort(all(pts));
48+
pts.resize(unique(all(pts)) - pts.begin());
49+
random_shuffle(all(pts));
50+
if (sz(pts) <=3) return pts;
51+
vector<double> line ={pts[1].y-pts[0].y, pts[0].x - pts[1].x, -pts[0].x*pts[1].y + pts[0].y*pts[1].x};
52+
array<vector<P>, 2> S;
53+
S[0].push_back(pts[0]);
54+
S[0].push_back(pts[1]);
55+
S[1].push_back(pts[1]);
56+
S[1].push_back(pts[0]);
57+
rep(i,2,sz(pts)) {
58+
int idx = (line[0]*pts[i].x + line[1]*pts[i].y + line[2]) >=0;
59+
S[idx].push_back(pts[i]);
60+
}
61+
auto ta = conquer(S[0],0);
62+
auto tb = conquer(S[1],0);
63+
auto pa=ta.second, pb=tb.second;
64+
if (!ta.first || !tb.first) {
65+
return genPolygon(pts);
66+
}
67+
pa.erase(pa.begin());
68+
pb.erase(pb.begin());
69+
pa.insert(pa.end(), all(pb));
70+
return pa;
71+
}

0 commit comments

Comments
 (0)