A well-tested C++
library for computational geometry problems. The code has bad style and conventions, but it’ll suffice for hackathons.
Types and constants
The numerical type involved in computation varies depending on the problem. So, for flexibility we’ll name this type gtype
that we’ll use throughout the code.
The macros equ
and inRange
should be self-explanatory, however, note that the comparison they make use some tolerance bound e
.
After that, follow some constants that are of type gtype
. len_before_point_digits
and ang_before_point_digits
are also problem dependent and they’re VERY critical. They define how many decimal digits before the decimal point we’re going to use. For example, if we know from our problem that distances will never exceed 100,000
, len_before_point_digits
should be 1E6
. The same applies for ang_before_point_digits
.
epsLen
, epsAng
and epsArea
are the tolerance bounds we’re likely to use for distances, radian angles and areas. Yep! It’s that simple.
typedef long double gtype;
#define equ(a, b, e) ( abs((a) - (b)) <= (e) )
#define inRange(a, b, p, e) ( (p) >= min(a, b) - (e) && (p) <= max(a, b) + (e))
// CONST
const gtype len_digits_before_point = 1E6;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_digits_before_point;
const gtype ang_digits_before_point = 1E2;
const gtype epsAng = numeric_limits<gtype>::epsilon() * ang_digits_before_point;
const gtype epsArea = epsLen * epsLen;
const gtype pi = M_PI;
Points
Let’s get down to business. Let’s define a point. What’s a point? It’s an ordered pair of some type. So, instead of writing it from scratch, we’ll use C++
‘s std::complex
which is also an ordered pair that can be given a numerical type. std::complex
spares us a lot of effort in overriding operators and has some very useful methods. We define the type std::complex <gtype>
as our point
.
To get a point
‘s x and y, we get std::complex
‘s real and imaginary parts respectively.
Another good feature of std::complex
is that it allows both Cartesian and polar constructors. But, the polar constructor must be given the correct type for the magnitude component, that’s why we defined the macro.
The last two macros do simple angle conversion.
typedef complex<gtype> point;
#define x real()
#define y imag()
#define polar(r, t) polar((gtype) (r), (gtype) (t))
bool cmp_point(point const& p1, point const& p2) {
return p1.x == p2.x ? (p1.y < p2.y) : (p1.x < p2.x);
}
#define deg(a) ( (a) * 180 /pi )
#define rad(a) ( (a) * pi /180 )
Vectors
Vectors are the building blocks of any geometrical structure. Rotations, products and transformations are very important in any problem. Fortunately, we don’t have to implement all of that. norm
, abs
, conj
, arg
are provided by std::complex
. We’ll use some of these functions and the rest we’ll code on our own. I’ve chosen the fastest and most precise functions I could lay my hands on:
// VECTOR
//#define norm(a) dot(a,a) // fast for [long double]
//#define abs(a) sqrt(norm(a)) // fast for [int]
#define arg(a) ( atan2((a).y, (a).x) )
#define unit(v) ( (v) / abs(v) )
#define dot(a, b) ( (conj(a) * (b)).x )
#define crs(a, b) ( (conj(a) * (b)).y )
#define rot(v, t) ( (v) * polar(1, t) )
#define rot2(v, t, o) ( rot((v) - (o), (t)) + (o) )
#define refelct(v, m) ( conj((v) / (m)) * (m) )
#define perp(v) ( point((v).y, -(v).x) )
#define vecVecProj(u, v) ( dot(u, v) * (v) / norm(v) )
Lines
Lines-point / line-line intersections are very frequently used. So, here are some useful macros:
// LINE
#define pntLinProj(a, b, p) ( (a) + vecVecProj((p) - (a), (b) - (a)) )
#define pntLinDist(a, b, p) ( abs(crs((b)-(a), (p)-(a)) / abs((b)-(a))) )
#define linIntrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define linIntrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )
#define linOnLin(a, b, p) ( abs(crs((b) - (a), (p) - (a))) <= epsLen )
#define linOnRay(a, b, p) ( linOnLin(a, b, p) && dot((b) - (a), (p) - (a)) >= -epsLen )
#define linOnSeg(a, b, p) ( linOnRay(b, a, p) && dot((b) - (a), (p) - (a)) >= -epsLen )
#define linInRect(a, b, p, e) ( inRange((a).x, (b).x, (p).x, e) && inRange((a).y, (b).y, (p).y, e) )
Triangles
Perimeter, surface, and cosine law:
// TRIANGLE
#define triPerim(a, b, c) ( (a) + (b) + (c) )
#define triSurface(a, b, c) ( (gtype) 0.5 * crs((b) - (a), (c) - (a)) )
#define triSurfaceH(s, a, b, c) ( sqrt((s) * ((s) - (a)) * ((s) - (b)) * ((s) - (c))) )
#define cosLaw(a, b, c) ( ((a) * (a) + (b) * (b) - (c) * (c)) / (2 * (a) * (b)) )
Circles
Arc length, chord length, sector surface area, and circle circle intersection:
// CIRCLE
#define cirArcLen(r, a) ((a) * abs(r))
#define cirSectorSurface(r, a) ((a) * norm(r) / 2)
#define cirChordLen(r, a) sqrt(2 * norm(r) * (1 - cos(a)))
int cirCircInter(point & c0, point & c1, gtype r0, gtype r1, vector<point> & out) {
point cT, v;
gtype rT, phi;
// init
if (r1 > r0)
rT = r0, r0 = r1, r1 = rT, cT = c0, c0 = point(c1), c1 = point(cT);
out.clear(), v = c1 - c0, phi = acos(cosLaw(r0, abs(v), r1));
// same center
if (equ(abs(v), 0, epsLen))
equ(r0, r1, epsLen) ? (equ(r0, 0, epsLen) ? out.push_back(c0) : out.resize(3)) : out.clear();
// none
else if (abs(abs(v) - r0) > r1 + epsLen)
out.clear();
// one point
else if (equ(abs(abs(v) - r0), r1, epsLen))
out.push_back(unit(v) * r0 + c0);
// two points
else
out.push_back(rot(unit(v) * r0, -phi) + c0), out.push_back(rot(unit(v) * r0, phi) + c0);
return out.size();
}
Polygons
Area, convexity check, and rotating calipers (min-width). Also includes convex hull creation:
// POLYGON (must be closed, start & end vertex overlap)
#define polAngle(n) (((n) - 2) * pi / (n))
long double polArea(vector<point> & p) {
long double area = 0;
for (size_t i = 0; i < p.size() - 1; ++i)
area += crs(p[i], p[i + 1]);
return abs(area * 0.5l);
}
bool polIsConvex(vector<point> & p) {
int a[2] = { 0, 0 }, n = p.size() - 1;
for (size_t i = 1; i < p.size(); ++i) {
long double cross = crs(p[i % n] - p[i - 1], p[(i + 1) % n] - p[i - 1]);
if (cross != 0)
++a[cross > 0];
}
return !(a[0] && a[1]);
}
// CONVEX HULL
// O(n.h) - jarvis's march
// TODO
// O(n.log(n)) - monotone chain
vector<point> mcH;
void monotoneChain(vector<point> &ps) {
vector<point> p(ps.begin(), ps.end() - 1);
int n = p.size(), k = 0;
mcH = vector<point>(2 * n);
sort(p.begin(), p.end(), cmp_point);
for (int i = 0; i < n; i++) {
while (k >= 2 && crs(mcH[k - 1] - mcH[k - 2], p[i] - mcH[k - 2]) <= 0)
k--;
mcH[k++] = p[i];
}
for (int i = n - 2, t = k + 1; i >= 0; i--) {
while (k >= t && crs(mcH[k - 1] - mcH[k - 2], p[i] - mcH[k - 2]) <= 0)
k--;
mcH[k++] = p[i];
}
mcH.resize(k);
}
// O(n) - rotating calipers (works on a ccw closed convex hull)
gtype rotatingCalipers(vector<point> &ps) {
int aI = 0, bI = 0;
for (size_t i = 1; i < ps.size(); ++i)
aI = (ps[i].y < ps[aI].y ? i : aI), bI = (ps[i].y > ps[bI].y ? i : bI);
gtype minWidth = ps[bI].y - ps[aI].y, aAng, bAng;
point aV = point(1, 0), bV = point(-1, 0);
for (gtype ang = 0; ang < pi; ang += min(aAng, bAng)) {
aAng = acos(dot(ps[aI + 1] - ps[aI], aV) / abs(aV) / abs(ps[aI + 1] - ps[aI]));
bAng = acos(dot(ps[bI + 1] - ps[bI], bV) / abs(bV) / abs(ps[bI + 1] - ps[bI]));
aV = rot(aV, min(aAng, bAng)), bV = rot(bV, min(aAng, bAng));
if (aAng < bAng)
minWidth = min(minWidth, pntLinDist(ps[aI], ps[aI] + aV, ps[bI])), aI = (aI + 1) % (ps.size() - 1);
else
minWidth = min(minWidth, pntLinDist(ps[bI], ps[bI] + bV, ps[aI])), bI = (bI + 1) % (ps.size() - 1);
}
return minWidth;
}