CP Geo
CP Geo
CP Geo
programmers
Victor Lecomte
http://creativecommons.org/licenses/by-sa/4.0/
1
Contents
2 Basics 27
2.1 Points and vectors . . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.1 Complex numbers . . . . . . . . . . . . . . . . . . . . 27
2.1.2 Point representation . . . . . . . . . . . . . . . . . . . 30
2.2 Transformations . . . . . . . . . . . . . . . . . . . . . . . . . 33
2
2.2.1 Translation . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.2 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.3 Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.4 General linear transformation . . . . . . . . . . . . . . 35
2.3 Products and angles . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.1 Dot product . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.2 Cross product . . . . . . . . . . . . . . . . . . . . . . . 38
2.4 Lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.4.1 Line representation . . . . . . . . . . . . . . . . . . . . 45
2.4.2 Side and distance . . . . . . . . . . . . . . . . . . . . . 47
2.4.3 Perpendicular through a point . . . . . . . . . . . . . 48
2.4.4 Sorting along a line . . . . . . . . . . . . . . . . . . . 49
2.4.5 Translating a line . . . . . . . . . . . . . . . . . . . . . 49
2.4.6 Line intersection . . . . . . . . . . . . . . . . . . . . . 50
2.4.7 Orthogonal projection and reflection . . . . . . . . . . 51
2.4.8 Angle bisectors . . . . . . . . . . . . . . . . . . . . . . 52
2.5 Segments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.5.1 Point on segment . . . . . . . . . . . . . . . . . . . . . 53
2.5.2 Segment-segment intersection . . . . . . . . . . . . . . 54
2.5.3 Segment-point distance . . . . . . . . . . . . . . . . . 56
2.5.4 Segment-segment distance . . . . . . . . . . . . . . . . 57
2.6 Polygons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.6.1 Polygon area . . . . . . . . . . . . . . . . . . . . . . . 58
2.6.2 Cutting-ray test . . . . . . . . . . . . . . . . . . . . . 59
2.6.3 Winding number . . . . . . . . . . . . . . . . . . . . . 62
2.7 Circles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.7.1 Defining a circle . . . . . . . . . . . . . . . . . . . . . 66
2.7.2 Circumcircle . . . . . . . . . . . . . . . . . . . . . . . 67
2.7.3 Circle-line intersection . . . . . . . . . . . . . . . . . . 68
2.7.4 Circle-circle intersection . . . . . . . . . . . . . . . . . 69
2.7.5 Tangent lines . . . . . . . . . . . . . . . . . . . . . . . 71
3
3 3D geometry 74
3.1 Points, products and orientation . . . . . . . . . . . . . . . . 74
3.1.1 Point representation . . . . . . . . . . . . . . . . . . . 74
3.1.2 Dot product . . . . . . . . . . . . . . . . . . . . . . . . 75
3.1.3 Cross product . . . . . . . . . . . . . . . . . . . . . . . 76
3.1.4 Mixed product and orientation . . . . . . . . . . . . . 78
3.2 Planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.2.1 Defining planes . . . . . . . . . . . . . . . . . . . . . . 80
3.2.2 Side and distance . . . . . . . . . . . . . . . . . . . . . 81
3.2.3 Translating a plane . . . . . . . . . . . . . . . . . . . . 82
3.2.4 Orthogonal projection and reflection . . . . . . . . . . 83
3.2.5 Coordinate system based on a plane . . . . . . . . . . 83
3.3 Lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.3.1 Line representation . . . . . . . . . . . . . . . . . . . . 85
3.3.2 Distance from a line . . . . . . . . . . . . . . . . . . . 87
3.3.3 Sorting along a line . . . . . . . . . . . . . . . . . . . 87
3.3.4 Orthogonal projection and reflection . . . . . . . . . . 87
3.3.5 Plane-line intersection . . . . . . . . . . . . . . . . . . 88
3.3.6 Plane-plane intersection . . . . . . . . . . . . . . . . . 89
3.3.7 Line-line distance and nearest points . . . . . . . . . . 90
3.4 Angles between planes and lines . . . . . . . . . . . . . . . . . 92
3.4.1 Between two planes . . . . . . . . . . . . . . . . . . . 92
3.4.2 Between two lines . . . . . . . . . . . . . . . . . . . . 93
3.4.3 Between a plane and a line . . . . . . . . . . . . . . . 94
3.4.4 Perpendicular through a point . . . . . . . . . . . . . 94
3.5 Polyhedrons . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.5.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.5.2 Surface area . . . . . . . . . . . . . . . . . . . . . . . . 96
3.5.3 Face orientation . . . . . . . . . . . . . . . . . . . . . 98
3.5.4 Volume . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.6 Spherical geometry . . . . . . . . . . . . . . . . . . . . . . . . 102
3.6.1 Spherical coordinate system . . . . . . . . . . . . . . . 102
3.6.2 Sphere-line intersection . . . . . . . . . . . . . . . . . 103
4
3.6.3 Great-circle distance . . . . . . . . . . . . . . . . . . . 104
3.6.4 Spherical segment intersection . . . . . . . . . . . . . 105
3.6.5 Angles on a sphere . . . . . . . . . . . . . . . . . . . . 110
3.6.6 Spherical polygons and area . . . . . . . . . . . . . . . 111
3.6.7 Solid angle . . . . . . . . . . . . . . . . . . . . . . . . 112
3.6.8 3D winding number . . . . . . . . . . . . . . . . . . . 114
5
Chapter 1
In this section we explore two ways in which very small starting imprecisions
can become very large imprecisions in the final output of a program.
6
1.1.1 When doing numerically unstable computations
There are some types of computations which can transform small impreci-
sions into catastrophically large ones, and line intersection is one of them.
Imagine you have four points A, B, C, D which were obtained through an
previous imprecise process (for example, their position can vary by a dis-
tance of at most r = 10−6 ), and you have to compute the intersection of
lines AB and CD.
For illustration, we represent the imprecisions on the points with small
disk of radius r: the exact position is the black dot, while the small gray
disk contains the positions it could take because of imprecisions. The dashed
circle gives an idea of where point I might lie.
In the best case, when A, B and C, D aren’t too close together, and not
too far from the intersection point I, then the imprecision on I isn’t too big.
A
C
But if those conditions are not respected, the intersection I might vary in
a very wide range or even fail to exist, if given the imprecision lines AB and
CD end up being parallel, or if A and B (or C and D) end up coinciding.
7
B
A
C
D
This shows that finding the intersection of two lines defined by imprecise
points is a task that is inherently problematic for floating-point arithmetic,
as it can produce wildly incorrect results even if the starting imprecision is
quite small.
8
This shows that even when only very small precision mistakes are made
(≈ 1.11 × 10−16 ), the overal loss of precision can get very big, and carefully
checking the maximal imprecision of your program is very important.
In this section, we explore ways in which small imprecisions can modify the
behavior of an algorithm in ways other than just causing further impreci-
sions.
The first scenario we will explore is when we have to make clear-cut decisions,
such as deciding if two objects touch.
Let’s say we have a line l and a point P computed imprecisely, and we
want to figure out if the point lies on the line. Obviously, we cannot simply
check if the point we have computed lies on the line, as it might be just
slightly off due to imprecision. So the usual approach is to compute the
distance from P to l and then figure out if that distance is less than some
small value like cutoff = 10−9 .
While this approach tends to works pretty well in practice, to be sure
that this solution works in every case and choose cutoff properly,1 we need to
know two things. First, we need to know error , the biggest imprecision that
we might make while computing the distance. Secondly, and more critically,
we need to know chance , the smallest distance that point P might be from l
while not being on it, in other words, the closest distance that it might be
from l “by coincidence”.
Only once we have found those two values, and made sure that error <
chance , can we then choose the value of cutoff within [error , chance ).2 Indeed,
if cutoff < error , there is a risk that P is on l but we say it is not, while if
cutoff ≥ chance there is a risk that P is not on l but we say it is.
Even though error can be easily found with some basic knowledge of
floating-point arithmetic and a few multiplications (see next section), find-
ing chance is often very difficult. It depends directly on which geometric
operations were done to find P (intersections, tangents, etc.), and in most
cases where chance can be estimated, it is in fact possible to make the com-
parison entirely with integers, which is of course the preferred solution.
1
And motivated problem setters do tend to find the worst cases.
2
In practice you should try to choose cutoff so that there is a factor of safety on both
sides, in case you made mistakes while computing error or chance .
9
1.2.2 By violating basic assumptions
In this section, we try to build a basic model of precision errors that we can
use to obtain a rough but reliable estimate of a program’s precision.
When the output of a problem is some real value like a distance or an area,
the problem statement often specifies a constraint such as: “The answer
should be accurate to an absolute or relative error of at most 10−5 .” While
considering the relative accuracy of an answer can be a useful and conve-
nient way to specify the required precision of an answer in some cases (for
example in tasks where only addition and multiplication of positive values
are performed), we think that for most geometry problems it is unsuitable.
The reason for this is the need to subtract3 large values of similar mag-
nitude. For example, suppose that we are able to compute two values with
relative precision 10−6 , such as A = 1000 ± 10−3 and B = 999 ± 10−3 . If we
compute their difference, we obtain A−B = 1±2×10−3 . The absolute error
remains of a comparable size, being only multiplied by 2, but on the other
hand relative error increases drastically from 10−6 to 2×10−3 because of the
decrease in magnitude. This phenomenon is called catastrophic cancellation.
In fact, whenever a certain relative error can affect big numbers, catas-
trophic cancellation can cause the corresponding absolute error to appear
on very small values. The consequence is that if a problem statement has a
3
Strictly speaking, we mean both subtraction of values of the same sign, and addition
of values of opposite signs.
10
certain tolerance on the relative error of the answer, and a correct solution
has an error close to it for the biggest possible values, then the problem
statement also needs to specify a tolerance on the corresponding absolute
error in case catastrophic cancellation happens. And since that tolerance on
absolute error is at least as tolerant as the tolerance on relative error for all
possible values, it makes it redundant. This is why we think that tolerance
on “absolute or relative error” is misleading at best.4
Catastrophic cancellation shows that relative precision is not a reliable
way to think about precision whenever subtractions are involved — and that
includes the wide majority of geometry problems. In fact, the most common
geometric operations (distances, intersections, even dot/cross products) all
involve subtractions of values which could be very similar in magnitude.
Examples of this appear in two of the case studies of section 1.4: in prob-
lem “Keeping the Dogs Apart” and when finding the solution of a quadratic
equation.
Another example occurs when computing areas of polygons made of
imprecise points. Even when the area ends up being small, the imprecision
on it can be large if there were computations on large values in intermediate
steps, which is the case when the coordinates have large magnitudes.
area area
absolute error ≈ absolute error
relative error relative error
11
1.3.2 Precision guarantees from IEEE 754
We explained earlier why we need to work with absolute error, but since
IEEE 754 gives us guarantees in terms of relative errors, we need to consider
the biggest magnitude that will be reached during the computations. In
other words, if all computations are precise up to a relative error of , and
the magnitude of the values never goes over M , then the absolute error of
an operation is at most M .
This allows us to give good guarantees for numbers obtained after a
certain number of + and − operations: a value that is computed in n opera-
tions6 will have an absolute error of at most nM compared to the theoretical
5
This assumes the magnitudes do not go outside the allowable range (≈ 10±308 for
±4932
double and ≈ 10 for long double) which almost never happens for geometry problems.
6
Note that when we say a value is “computed in n operations” we mean that it is
computed by a single formula that contains n operations, and not that n operations are
necessary to actually compute it. For example (a+b)+(a+b) is considered to be “computed
in 3 operations” even though we can implement this with only 2 additions.
12
result.
We can prove the guarantee by induction: let’s imagine we have two
intermediate results a and b who were computed in na and nb operations
respectively. By the inductive hypothesis their imprecise computed values
a0 and b0 respect the following conditions.
|a0 − a| ≤ na M |b0 − b| ≤ nb M
| round(a0 + b0 ) − (a + b)|
= round(a0 + b0 ) − (a0 + b0 ) + (a0 + b0 ) − (a + b)
≤ | round(a0 + b0 ) − (a0 + b0 )| + |(a0 + b0 ) − (a + b)|
≤ M + |(a0 − a) + (b0 − b)|
≤ M + |a0 − a| + |b0 − b|
≤ M + na M + nb M
= (na + nb + 1)M
where the first two steps follow from the triangle inequality. Since the sum
is “computed in na + nb + 1 operations”, the bound of (na + nb + 1)M that
is obtained is small enough. The proof for subtraction is very similar.
The model above gives good guarantees but is very limited: it only works for
computations that use only addition and subtraction. Multiplication does
not give guarantees of the form nM . However, we can still say interesting
things if we take a closer look the different types of values we use in geometry:
• Adimensional “0D” values: e.g. angles, constant factors;
• 1D values: e.g. coordinates, lengths, distances, radii;
• 2D values: e.g. areas, dot products, cross products;
• 3D values: e.g. volumes, mixed products.
Usually, the problem statement gives guarantees on the magnitude of
coordinates, so we can find some constant M so that all 1D values that will
be computed in the code have a magnitude less than M . And since 2D and
3D values are usually created by products of 1D values, we can usually say
that 2D values are bounded in magnitude by M 2 and 3D values by M 3 (we
may need to multiply M by a constant factor).
13
It turns out that computations made of +, −, × and in which all d-
dimensional values are bounded in magnitude by M d have good precision
guarantees. In fact, we can prove that the absolute error of a d-dimensional
number computed in n operations is at most M d ((1 + )n − 1), which as-
suming n 1 is about nM d .
The proof is similar in spirit to what we did with only + and − earlier.
Since it is a bit long, we will not detail it here, but it can be found in
section B.1, along with a more precise definition of the precision guarantees
and its underlying assumptions.
Note that this does not cover multiplication by an adimensional factor
bigger than 1: this makes sense, since for example successive multiplication
by 2 of a small value could make the absolute error grow out of control even
if the magnitude remains under M d for a while.
In other cases, this formula nM d gives us a quick and reliable way to
estimate precision errors.
14
√ √
x x
δ x x
δ
close to 0 ⇒ big effect far from 0 ⇒ small effect
√
Therefore for
√ a given imprecision δ, the biggest imprecision on x it
might cause is δ. This is usually pretty bad: if the argument of the square
root had an imprecision of nM 2 then in the worst case the result will have
√ √
an imprecision of nM , instead of the nM bound that we have for
+, −, × operations.
For example let us consider a circle C of radius tangent to a line l. If C
gets closer to l by 10−6 , then the intersection points will move by about
p p √
12 − (1 − 10−6 )2 ≈ 2 × 10−6 = 2 × 10−3
away from the tangency point, as pictured below.
l C l C
√
Note that here we have only shown that 1/x and x perform poorly on
imprecise inputs. Please bear in mind that on exact inputs, the IEEE 754
guarantees that the result is the closest represented floating-point number.
So when the lines and circles are defined by integers, line intersections and
circle-line intersections have a relative precision error proportional to and
thus an absolute error proportional to M .
We will first talk about problem “Keeping the Dogs Apart”, which we
mentioned before, because it is a good example of accumulation of error
15
and how to deal with it. It was written by Markus Fanebust Dregi for
NCPC 2016. You can read the full statement and submit it at https:
//open.kattis.com/problems/dogs.
16
s = moveBy(b[j], b[j+1], end-sumB);
// We get to the end of the segment for one dog or the other,
// so move to the next and update the sum of lengths
if (endA < endB) {
i++;
sumA += dA;
} else {
j++;
sumB += dB;
}
}
// output ans
The first way we can limit the accumulation of error in sumA and sumB is to
realize that in fact, we only care about the difference between them: if we
add a certain constant to both variables, this doesn’t change the value of
start-sumA, end-sumA, start-sumB or end-sumB, so the value of p, q, r, s is
unchanged.
So we can adapt the code by adding these lines at the end of the while
loop:
double minSum = min(sumA, sumB);
sumA -= minSum;
sumB -= minSum;
17
After this, one of sumA and sumB becomes zero, while the other carries
the error on both. In total, at most n + m additions and n + m subtractions
are performed on them, for a total of k ≤ 4 × 105 . But since the difference
√ sumA 4and sumB never exceeds the length of one segment, that is,
between
M = 2 × 10 , the error is much lower than before:
√
kM = 4 × 105 × 2 × 104 × 2−53 ≈ 6.3 × 10−7
so it gives an AC verdict.
So here we managed to reduce the precision mistakes on our results by
reducing the magnitude of the numbers that we manipulate. Of course, this
is only possible if the problem allows it.
Now we present different way to reduce the precision mistake, based on the
fact that all the terms in the sum we’re considering are positive. This is a
good thing, because it avoids catastrophic cancellation (see section 1.3.1).
In fact, addition of nonnegative numbers conserves relative precision: if
you sum two nonnegative numbers a and b with relative errors of ka and kb
respectively, the worst-case relative error on a+b is about7 (max(ka , kb )+1).
Let’s say we need to compute the sum of n nonnegative numbers a1 , . . . , an .
We suppose they are exact. If we perform the addition in the conventional
order, like this:
!
· · · (a1 + a2 ) + a3 + a4 + · · · + an
then
• a1 + a2 will have a relative error of (max(0, 0) + 1) = ;
• (a1 + a2 ) + a3 will have a relative error of (max(1, 0) + 1) = 2;
• (a1 + a2 ) + a3 + a4 will have a relative error of (max(2, 0) + 1) = 3;
• . . . and so on.
So the complete sum will have an error of (n − 1), not better than what we
had before.
But what if we computed the additions in another order? For example,
with n = 8, we could do this:
(a1 + a2 ) + (a3 + a4 ) + (a5 + a6 ) + (a7 + a8 )
7
It could in fact go up to (max(ka , kb )(1 + ) + 1) but the difference is negligible for
our purposes.
18
then all additions of two numbers have error , all additions of 4 num-
bers have error (max(1, 1) + 1) = 2, and the complete addition has error
(max(2, 2) + 1) = 3, which is much better than (n − 1) = 7. In general,
for n = 2k , we can reach a relative precision of k.
We can use this grouping technique to create an accumulator such that
the relative error after adding n numbers is at most 2 log2 (n).8 Here is an
O(n) implementation:
struct stableSum {
int cnt = 0;
vector<double> v, pref{0};
void operator+=(double a) {
assert(a >= 0);
int s = ++cnt;
while (s % 2 == 0) {
a += v.back();
v.pop_back(), pref.pop_back();
s /= 2;
}
v.push_back(a);
pref.push_back(pref.back() + a);
}
double val() {return pref.back();}
};
Let’s break this code down. This structure provides two methods: add(
a) to add a number a to the sum, and val() to get the current value of the
sum. Array v contains the segment sums that currently form the complete
sum, similar to Fenwick trees: for example, if we have added 11 elements, v
would contain three elements:
while pref contains the prefix sums of v: pref[i] contains the sum of the i
first elements of v.
Function add() performs the grouping: when adding a new element a, it
will merge it with the last element of v while they contain the same number
of terms, then a is added to the end of v. For example, if we add the 12th
8
We could even get (log2 n+1) but we don’t know a way to do it faster than O(n log n).
19
element a12 , the following steps will happen:
The number of additions we have to make for the ith number is the
number of times it is divisible by 2. Since we only add one element to v
when adding an element to the sum, this is amortized constant time.
By simply changing the types of sumA and sumB to stableSum and adding
.val() whenever the value is read in the initial code, we can get down to an
error of about
√
2 log2 105 M = 2 log2 105 × 2 × 109 × 2−53 ≈ 5.2 × 10−6
As another example, we will study the precision problems that can occur
when computing the roots of an equation of the type ax2 + bx + c = 0 with
a 6= 0. We will see how some precision problems are unavoidable, while
others can be circumvented by
When working with full-precision reals, we can solve quadratic equations
in the following way. First we compute the discriminant ∆ = b2 − 4ac. If
∆ < 0, there is no solution, while if ∆ ≥ 0 there is are 1 or 2 solutions, given
by √
−b ± ∆
x=
2a
The first difficulty when working with floating-point numbers is the com-
putation of ∆: if ∆ ≈ 0, that is b2 ≈ 4ac, then the imprecisions can change
the sign of ∆, therefore changing the number of solutions.
Even if that does not happen, since we have to perform a square root,
the problems that we illustrated with line-circle intersection in section 1.3.5
can also happen here.10 Take the example of equation x2 − 2x + 1 = 0,
9
Theoretically we can’t really be sure though, since both sumA and sumB could have that
error, and we still have to take into account the other operations performed.
10
Which is not surprising, since the bottom of a parabola looks a lot like a circle.
20
which is a single root x = 0. If there is a small error on c, it can translate
into a large error on the roots. For example, if c = 1 − 10−6 , then the roots
become
√ p √
−b ± b2 − 4ac 2 ± 4 − 4(1 − 10−6 ) 2 ± 2 10−6
x= = = = 1 ± 10−3 .
2a 2 2
where the error 10−3 is much bigger than the initial error on c.
x2 − 2x + (1 + ) = 0 x2 − 2x + 1 = 0 x2 − 2x + (1 − ) = 0
√
no solution x=1 x=1±
√
Even if the √ computation of ∆ is very precise, a second problem can
occur. If b and ∆ have a similar magnitude, in other words when b2 ac,
then catastrophic cancellation will occur for one of the roots. For example
if a = 1, b = 104 , c = 1, then the roots will be:
√ √
−104 − 108 − 4 −104 + 108 − 4
x1 = ≈ −10 4
x2 = ≈ 10−4
2 2
√
The computation of x1 goes fine because √ −b and − ∆ have the same
sign. But because the√magnitude of −b + ∆ is 108 times smaller than the
magnitude of b and ∆, the relative √ error on x2 will be 108 times bigger
than the relative error on b and ∆.
Fortunately, in this case we can avoid catastrophic cancellation entirely
by rearranging the expression:
√ √ √
−b ± b2 − 4ac −b ± b2 − 4ac −b ∓ b2 − 4ac
= × √
2a 2a −b ∓ b2 − 4ac
√ 2
(−b2 ) − b2 − 4ac
= √
2a −b ∓ b2 − 4ac
4ac
= √
2a −b ∓ b2 − 4ac
2c
= √
−b ∓ b2 − 4ac
In this new expression, since the sign of the operation is opposite from the
sign in the original expression, catastrophic cancellation happens in only one
of the two.
21
√
−b− ∆ 2c√
So if b ≥ 0, we can use 2a for the first solution and −b− ∆
for
2c√
the second solution, while if b ≤ 0, we can use −b+ ∆
for the first solution
√
and −b+
2a
∆
for the second solution. We only need to be careful that the
denominator is never zero.
This gives a safer way to find the roots of a quadratic equation. This
function returns the number of solutions, and places them in out in no
particular order.
int quadRoots(double a, double b, double c, pair<double,double> &out
) {
assert(a != 0);
double disc = b*b - 4*a*c;
if (disc < 0) return 0;
double sum = (b >= 0) ? -b-sqrt(disc) : -b+sqrt(disc);
out = {sum/(2*a), sum == 0 ? 0 : (2*c)/sum};
return 1 + (disc > 0);
}
In many cases, there are several ways to write an expression, and they
can have very different behaviors when used with floating-point numbers. So
if you realize that the expression you are using can cause precision problems
in some cases, it can be a good idea to rearrange the expression to handle
them, as we did here.
This last case study will study one possible implementation for the inter-
section of two circles. It will show us why we shouldn’t rely too much on
mathematical truths when building our programs.
We want to know whether two circles of centers C1 , C2 and radii r1 , r2
touch, and if they do what are the intersection points. Here, we will solve
this problem with triangle inequalities and the cosine rule.11 Let d = |C1 C2 |.
The question of whether the circles touch is equivalent to the question
of whether there exists a (possibly degenerate) triangle with edge lengths
d, r1 , r2 .
11
The way we actually implement it in this book is completely different.
22
r1 r2
α
C1 C2
d
We know that such a triangle exists iff the triangle inequalities are re-
spected, that is:
|r2 − r1 | ≤ d ≤ r1 + r2
If this is true, then we can find the angle at C1 , which we’ll call α, thanks
to the cosine rule:
d2 + r12 − r22
cos α =
2dr1
Once we have α, we can find the intersection points in the following way:
# »
if we take vector C1 C2 , resize it to have length r1 , then rotate by α in either
direction, this gives the vectors from C1 to either intersection points.
C1 C2 C1 C2 C1 C2
This gives the following code. It uses a function abs() to compute the
length of a vector (see section 2.1.2) and a function rot() to rotate a vector
by a given angle (see section 2.2.3).
bool circleCircle(pt c1, double r1, pt c2, double r2, pair<pt,pt> &
out) {
double d = abs(c2-c1);
if (d < abs(r2-r1) || d > r1+r2) // triangle inequalities
return false;
double alpha = acos((d*d + r1*r1 - r2*r2)/(2*d*r1));
pt rad = (c2-c1)/d*r1; // vector C1C2 resized to have length d
out = {c1 + rot(rad, -alpha), c1 + rot(rad, alpha)};
return true;
}
23
This implementation is quite nice, but unfortunately it will sometimes
output nan values. In particular, if
then the triangle inequalities are respected, so the function returns true,
but the program computes
d2 + r12 − r22
>1
2dr1
24
double d = abs(c2-c1), co = (d*d + r1*r1 - r2*r2)/(2*d*r1);
if (abs(co) > 1) return false;
double alpha = acos(co);
pt rad = (c2-c1)/d*r1; // vector C1C2 resized to have length d
out = {c1 + rot(rad, -alpha), c1 + rot(rad, alpha)};
return true;
}
25
r2 − d2 and (r + d)(r − d) are not always exactly the same value. Be
extra careful when those values are then used in an operation that is
√
not defined everywhere (like x, arccos(x), tan(x), xy etc.).
In general, try to build programs that are resistant to the oddities of
floating-point numbers. Imagine that some evil demon is slightly modifying
every result you compute in the way that is most likely to make your program
fail. And try to write clean code that is clearly correct at first glance. If you
need long explanations to justify why your program will not fail, then it is
more likely that your program will in fact fail.
Finally, here is some general advice about precision issues when creating a
geometry problem and its datasets.
• Never use floating-point numbers as inputs, as this will already cause
imprecisions when first reading the input numbers, and completely
exclude the use of integers make it impossible to determine some things
with certainty, like whether two segmentns touch, whether some points
are collinear, etc.
• Make the magnitude of the input coordinates as small as possible to
avoid causing overflows or big imprecisions in the contestant’s codes.
• Favor problems where the important comparisons can be made entirely
with integers.
• Avoid situations in which imprecise points are used for numerically
unstable operations such as finding the intersection of two lines.
• In most cases, you should specify the tolerance in terms of absolute
error only (see subsection 1.3.1).
• Make sure to prove that all correct algorithm are able to reach the
precision that you require, and be careful about operations like circle-
line intersection which can greatly amplify imprecisions. Since error
analysis is more complicated than it seems at first sight and requires
a bit of expertise, you may want to ask a friend for a second opinion.
26
Chapter 2
Basics
In this section, we will first introduce complex numbers, because they are a
useful way to think about and represent 2D points, especially when rotations
are involved. We will then present two different ways to represent points
in code: one by creating our own structure, the other by using the C++
built-in complex type. Either can be used to run the code samples in this
book, though complex requires less typing.
Complex numbers are an extension of the real numbers with a new unit,
the imaginary unit, noted i. A complex number is usually written as a + bi
(for a, b ∈ R) and we can interpret it geometrically as point (a, b) in the
two-dimensional plane, or as a vector with components #» v = (a, b). We
will sometimes use all these notations interchangeably. The set of complex
numbers is written as C.
3i 5 + 3i
3)
(5,
#»v =
27
Basic operations
#»w #» v#»−
w w#»
#»v + #»
#»
w w
#» #»
v
v
v#»−
w#»
addition subtraction
#»
v
#»
v
#»
2v
multiplication by scalar
Polar form
28
(r sin ϕ)i b2 a + bi
√ a2 +
r= ϕ
r cos ϕ
For a given complex number a + bi, we can compute its polar form as
p
r = |a + bi| = a2 + b2
ϕ = arg(a + bi) = atan2(b, a)
and conversely, a complex number with polar coordinates r, ϕ can be written
r cos ϕ + (r sin ϕ)i = r(cos ϕ + i sin ϕ) =: r cis ϕ
where cis ϕ = cos ϕ + i sin ϕ is the unit vector that forms an angle of ampli-
tude ϕ with the positive real axis.
Note that this is not a one-to-one mapping. Firstly, adding or subtracting
2π from ϕ doesn’t change the point being represented; to solve this problem,
ϕ is generally taken in (−π, π]. Secondly, when r = 0, all values of ϕ
represent the same point.
Multiplication
v#»
∗ w#» #»
w #»
v
29
Remarkably, multiplication is also very simple to compute from the co-
ordinates: it works a bit like polynomial multiplication, except that we
transform i2 into −1.
Exercise 1
Prove that (r1 cis ϕ1 ) ∗ (r2 cis ϕ2 ) = (r1 r2 ) cis(ϕ1 + ϕ2 ) using this new
definition of product.
[Go to solution]
#»
v
v ∗ 12 i
#»
ϕ+ π
2
ϕ
In this section we explain how to implement the point structure that we will
use throughout the rest of the book. The code is only available in C++ at
the moment, but should be easy to translate in most languages.
Let us first declare the basic operations: addition, subtraction, and multi-
plication/division by a scalar.
typedef double T;
struct pt {
30
T x,y;
pt operator+(pt p) {return {x+p.x, y+p.y};}
pt operator-(pt p) {return {x-p.x, y-p.y};}
pt operator*(T d) {return {x*d, y*d};}
pt operator/(T d) {return {x/d, y/d};} // only for floating-
point
};
For generality, we declare type T: the type of the the coordinates. Generally,
either double or long long (for exact computations with integers) is appro-
priate. long double can also be very useful if extra precision is required.
The cases where integers cannot be used are often quite clear (e.g. division
by scalar, rotation by arbitrary angle).
We define some comparators for convenience:
bool operator==(pt a, pt b) {return a.x == b.x && a.y == b.y;}
bool operator!=(pt a, pt b) {return !(a == b);}
The squared absolute value sq() can be used to compute and compare dis-
tances quickly and exactly if the coordinates are integers. We use double
for abs() because it will return floating-point values even for integer co-
ordinates (if you are using long double you should probably change it to
long double).
We also define a signum function, which will be useful for several appli-
cations. It returns -1 for negative numbers, 0 for zero, and 1 for positive
numbers.
31
template <typename T> int sgn(T x) {
return (T(0) < x) - (x < T(0));
}
Using the complex type in C++ can be a very practical choice in contests
such as ACM-ICPC where everything must be typed from scratch, as many
of the operations we need are already implemented and ready to use.
The code below defines a pt with similar functionality.
typedef double T;
typedef complex<T> pt;
#define x real()
#define y imag()
Warning
As with the custom structure, you should choose the appropriate
coordinate type for T. However, be warned that if you define it as
an integral type like long long, some functions which should always
return floating-point numbers (like abs() and arg()) will be truncated
to integers.
The macros x and y are shortcuts for accessing the real and imaginary
parts of a number, which are used as x- and y-coordinates:
pt p{3,-4};
cout << p.x << " " << p.y << "\n"; // 3 -4
// Can be printed out of the box
cout << p << "\n"; // (3,-4)
We can perform all the operations that we have with the custom struc-
ture and then some more. Of course, we can also use complex multiplication
and division. Note however that we can only multiply/divide by scalars of
type T (so if T is double, then int will not work).
pt a{3,1}, b{1,-2};
a += 2.0*b; // a = (5,-3)
32
cout << a*b << " " << a/-b << "\n"; // (-1,-13) (-2.2,-1.4)
There are also useful methods for dealing with polar coordinates:
pt p{4,3};
// Get the absolute value and argument of point (in [-pi,pi])
cout << abs(p) << " " << arg(p) << "\n"; // 5 0.643501
// Make a point from polar coordinates
cout << polar(2.0, -M_PI/2) << "\n"; // (1.41421,-1.41421)
Warning
The complex library provides function norm, which is mostly equivalent
to the sq that we defined earlier. However, it is not guaranteed to be
exact for double: for example, the following expression evaluates to
false.
norm(complex<double>(2.0,1.0)) == 5.0
2.2 Transformations
2.2.1 Translation
#»v
33
The implementation is self-explanatory:
pt translate(pt v, pt p) {return p+v;}
2.2.2 Scaling
2.2.3 Rotation
ϕ
c ϕ
34
pt rot(pt p, double a) {
return {p.x*cos(a) - p.y*sin(a), p.x*sin(a) + p.y*cos(a)};
}
And among those, we will use the rotation by 90◦ quite often:
It is easy to check that all those transformations are of the form f (p) =
a∗p+b as claimed in the beginning of this section. In fact, all transformations
of this type can be obtained as combinations of translations, scalings and
rotations.1
Just like for real numbers, to determine a linear transformation such as
this one, we only need to know the image of two points to know the complete
function. Indeed, if we know f (p) = a ∗ p + b and f (q) = a ∗ q + b, then we
can find a as f (q)−f
q−p
(p)
, and then b as f (p) − a ∗ p.
And thus if we want to know a new point f (r) of that transformation,
we can then compute it as:
f (q) − f (p)
f (r) = f (p) + (r − p) ∗
q−p
f (q)
f (r)
r
f (p)
p q
35
pt linearTransfo(pt p, pt q, pt r, pt fp, pt fq) {
return fp + (r-p) * (fq-fp) / (q-p);
}
Otherwise, you can use the cryptic but surprisingly short solution from
[2] (see the next sections for dot() and cross()):
pt linearTransfo(pt p, pt q, pt r, pt fp, pt fq) {
pt pq = q-p, num{cross(pq, fq-fp), dot(pq, fq-fp)};
return fp + pt{cross(r-p, num), dot(r-p, num)} / sq(pq);
}
Besides complex multiplication, which is nice to have but is not useful so of-
ten, there are two products involving vectors that are of critical importance:
dot product and cross product. In this section, we’ll look at their definition,
properties and some basic use cases.
#» #»
w
w θ θ
θ
#»
v #»
v #»
w #»
v
θ < π/2 θ = π/2 θ > π/2
#»
v ·w#» = 5 #»
v ·w#» = 0 #»
v ·w#» = −5
36
If we fix k #» #» as above, the dot product is maximal when the
v k and k wk
vectors point in the same direction, because cos θ = cos(0) = 1, and minimal
when they point in opposite directions, because cos θ = cos(π) = −1.
Math insight
Because of the definition of cosine in right triangles, dot prod-
uct can be interpreted in an interesting way (assuming #» v 6= 0):
#» #» #» #» #»
v · w = k v kπ v ( w) where π v ( w)
#» #» := #»
k wk cos θ is the signed length
of the projection of w #» onto the line that contains #»
v (see examples be-
low). In particular, this means that the dot product does not change
if one of the vectors moves perpendicular to the other.
#»
w
#» )
w #»
v #»
v
π #»v ( θ
θ
#» »)
#w
w π #»v (
#» #»
π #»
v ( w) = 1.58 v ( w) = −1.26
π #»
Dot product is often used for testing if two vectors are perpendicular,
since we just need to test whether #» #» = 0:
v ·w
bool isPerp(pt v, pt w) {return dot(v,w) == 0;}
It can also be used for finding the angle between two vectors, in [0, π].
Because of precision errors, we need to be careful not to call acos with a
value that is out of the allowable range [−1, 1].
double angle(pt v, pt w) {
double cosTheta = dot(v,w) / abs(v) / abs(w);
return acos(max(-1.0, min(1.0, cosTheta)));
}
37
2.3.2 Cross product
#»
v
θ
#» #» θ
w v #»
w
θ
#»
v #»
w
0<θ<π θ=π −π < θ < 0
#» #» = 5
v ×w #»
v ×w#» = 0 #»
v ×w#» = −7
38
Math insight
Because of the definition of sine in right triangles, cross product can
also be interpreted in an interesting way (assuming v 6= 0): #» v ×w #» =
k #»
v kd #» #» #» #»
v ( w) = k wk sin θ is the signed distance from the
v ( w), where d #»
line that contains v , with positive values on the left side of #»
#» v . In
particular, this means that the cross product doesn’t change if one of
the vectors moves parallel to the other.
d v#»( #»
w)
#»
v #»
v
#»
w θ
θ
#»
w
d v#»( #»
w)
#» #»
d #»
v ( w) = 2.69 v ( w) = −2.37
d #»
Like dot product, cross product has a very simple expression in cartesian
coordinates: if #» #» = (w , w ), then #»
v = (vx , vy ) and w x y
#» = v w − v w :
v ×w x y y x
Implementation trick
When using complex, we can implement both dot() and cross() with
this trick, which is admittedly quite cryptic, but requires less typing
and is less prone to typos:
T dot(pt v, pt w) {return (conj(v)*w).x;}
T cross(pt v, pt w) {return (conj(v)*w).y;}
and see that the real and imaginary parts are indeed the dot product
and the cross product.
39
Orientation
One of the main uses of cross product is in determining the relative position
of points and other objects. For this, we define the function orient(A, B, C) =
# » # » # »
AB × AC. It is positive if C is on the left side of AB, negative on the right
# »
side, and zero if C is on the line containing AB. It is straightforward to
implement:
T orient(pt a, pt b, pt c) {return cross(b-a,c-a);}
C B
C
A B A
B A
C
left turn collinear right turn
orient(A, B, C) > 0 orient(A, B, C) = 0 orient(A, B, C) < 0
C B
B C
A A
orient(A, B, C) > 0 orient(A, B, C) < 0
no swap swap
40
P P
C C C
B B B
A A P A
orient(A, B, P ) ≥ 0 orient(A, B, P ) < 0 orient(A, B, P ) ≥ 0
orient(A, C, P ) > 0 orient(A, C, P ) ≤ 0 orient(A, C, P ) ≤ 0
KO KO OK
bool inAngle(pt a, pt b, pt c, pt p) {
assert(orient(a,b,c) != 0);
if (orient(a,b,c) < 0) swap(b,c);
return orient(a,b,p) >= 0 && orient(a,c,p) <= 0;
}
B
C
B A
A
C
double orientedAngle(pt a, pt b, pt c) {
if (orient(a,b,c) >= 0)
return angle(b-a, c-a);
else
return 2*M_PI - angle(b-a, c-a);
}
41
Yet another use case is checking if a polygon P1 · · · Pn is convex: we com-
pute the n orientations of three consecutive vertices orient(Pi , Pi+1 , Pi+2 ),
wrapping around from n to 1 when necessary. The polygon is convex if they
are all ≥ 0 or all ≤ 0, depending on the order in which the vertices are given.
bool isConvex(vector<pt> p) {
bool hasPos=false, hasNeg=false;
for (int i=0, n=p.size(); i<n; i++) {
int o = orient(p[i], p[(i+1)%n], p[(i+2)%n]);
if (o > 0) hasPos = true;
if (o < 0) hasNeg = true;
}
return !(hasPos && hasNeg);
}
Polar sort
42
v#»4
v#»5
v#»6 v#»3
v#»0
v#»1 v#»2
In general #»
v should go before w #» when #»
v ×w #» > 0, because that means
w#» is to the left of #»
v when looking from the origin. This test works well for
directions that are sufficiently close: for example, v#»2 × v#»3 > 0. But when
they are more than 180◦ apart in the order, it stops working: for example
v#»1 × v#»5 < 0. So we first need to split the points in two halves according to
their argument:
v#»4
v#»5
v#»6
v#»3
v#»0
v#»2
v#»1
43
• Right now, points that are in the exact same direction are considered
equal, and thus will be sorted arbitrarily. If we want, we can use their
magnitude as a tie breaker:
void polarSort(vector<pt> &v) {
sort(v.begin(), v.end(), [](pt v, pt w) {
return make_tuple(half(v), 0, sq(v)) <
make_tuple(half(w), cross(v,w), sq(w));
});
}
With this tweak, if two points are in the same direction, the point that
is further from the origin will appear later.
• We can perform a polar sort around some point O other than the
origin: we just have to subtract that point O from the vectors #» v and
#» when comparing them. This as if we translated the whole plane so
w
that O is moved to (0, 0):
void polarSortAround(pt o, vector<pt> &v) {
sort(v.begin(), v.end(), [](pt v, pt w) {
return make_tuple(half(v-o), 0)) <
make_tuple(half(w-o), cross(v-o, w-o));
});
}
#»
v
44
2.4 Lines
In this section we will discuss how to represent lines and a wide variety of
applications.
Lines are sets of points (x, y) in the plane which obey an equation of the
form ax + by = c, with at least one of a and b nonzero. a and b determine
the direction of the line, while c determines its position.
c +2
y=
+b +1
ax c
y=
+b c
ax y=
b
ax+ −1
c
y=
+b −2
ax c
y=
+b
ax
(a, b)
c
b y=
ax+
45
a)
(b, −
#»v =
Q
P y =c
+b
ax
In this way, finding the equation of a line going through two points P
# »
and Q is easy: define the direction vector #»
v = (b, −a) = P Q, then find c as
#»
v × P.
struct line {
pt v; T c;
// From direction vector v and offset c
line(pt v, T c) : v(v), c(c) {}
// From equation ax+by=c
line(T a, T b, T c) : v({b,-a}), c(c) {}
// From points P and Q
line(pt p, pt q) : v(q-p), c(cross(v,p)) {}
A line will always have some implicit orientation, with two sides: the
“positive side” of the line (ax + by − c > 0) is on the left of #»
v , while the
“negative side” (ax+by −c < 0) is on the right of #»
v . In our implementation,
this orientation is determined by the points that were used to create the line.
We will represent it by an arrow at the end of the line. The figure below
shows the differences that occur when creating a line as line(p,q) or line(
q,p).
46
“positive side” “negative side”
ax + by − c > 0 ax + by − c < 0
Q Q
P P
“positive side”
sidel (P ) > 0
=0
P)
e l(
sid
“negative side”
sidel (P ) < 0
47
A
distl (A)
distl (B)
Two lines are perpendicular if and only if their direction vectors are per-
pendicular. Let’s say we have a line l of direction vector #»v . To find a line
perpendicular to line l and which goes through a certain point P , we could
define its direction vector as perp( #» v rotated by 90◦ counter-
v ) (that is, #»
clockwise, see section 2.2.3) and then try to work out c. However, it’s simpler
to just compute it as the line from P to P + perp( #» v ).
#»
v
l P + perp( #»
v)
perp( #»
v)
P
48
2.4.4 Sorting along a line
One subtask that often needs to be done in geometry problems is, given
points on a line l, to sort them in the order they appear on the line, following
the direction of #»
v.
#»
v
4
3
2
1
0
We can use the dot product to figure out the order of two points: a point
A comes before a point B if #»
v · A < #»v · B. So we can a comparator out of
it.
bool cmpProj(pt p, pt q) {
return dot(v,p) < dot(v,q);
}
0
c l0
y= l
+b
#» ax
P+ t #»
t
y =c
+b
ax
P
To find its new value c0 , we can see that for some point P on l, then
#»
P + t must be on the translated line. So we have
sidel (P ) = #»
v ×P −c=0
#» #»
sidel0 P + t = v × P + t − c0 = 0
#»
3
If you don’t know what an orthogonal projection is, read section 2.4.7.
49
which allows us to find c0 :
#» #» #»
c0 = #»
v × (P + t ) = #»
v × P + #»
v × t = c + #»
v × t
0
= c l0
+ by
ax l
c
δ b y=
ax+
There is a unique intersection point between two lines l1 and l2 if and only
if v# l»1 × v# l»2 6= 0. If it exists, we will show that it is equal to
cl1 v# l»2 − cl2 v# l»1
P =
v# l»1 × v# l»2
l1
P
l2
50
We only show that it lies on l1 (it should be easy to see that the expression
is actually symmetric in l1 and l2 ). It suffices to see that sidel1 (P ) = 0:
#»
# » cl1 vl2 − cl2 v# l»1
sidel1 (P ) = vl1 × − cl1
v# l»1 × v# l»2
cl (v# l» × v# l»2 ) − cl2 (v# l»1 × v# l»1 ) − cl1 (v# l»1 × v# l»2 )
= 1 1
v# l»1 × v# l»2
−cl2 (v# l»1 × v# l»1 )
=
v# » × v# »
l1 l2
=0
P
P
l
l
P0
P0
projection reflection
51
so we find k = − sidel (P )/k #»
v k2 .
pt proj(pt p) {return p - perp(v)*side(p)/sq(v);}
To find the reflection, we need to move P in the same direction but twice
the distance:
pt refl(pt p) {return p - perp(v)*2*side(p)/sq(v);}
lext (l1 , l2 )
l2
lint (l1 , l2 )
l1
sidel1 (P ) sidel2 (P )
# » =−
kvl1 k kv# l»2 k
52
This is exactly an expression of the form sidel (P ) = #»
v × P − c = 0 which
defines the points on a line. So it means that we have found the #» v and c
that characterize lint (l1 , l2 ):
#» v# l» v# l»
v = # »1 + # »2
kvl1 k kvl2 k
cl cl
c = # »1 + # »2
kvl1 k kvl2 k
The reasoning is very similar for lext (l1 , l2 ), the only difference being
signs. Both can be implemented as follows.
line bisector(line l1, line l2, bool interior) {
assert(cross(l1.v, l2.v) != 0); // l1 and l2 cannot be parallel!
double sign = interior ? 1 : -1;
return {l2.v/abs(l2.v) + l1.v/abs(l1.v) * sign,
l2.c/abs(l2.v) + l1.c/abs(l1.v) * sign};
}
2.5 Segments
53
P
A B A B
# » # » # » # »
PA·PB ≤ 0 PA·PB > 0
in disk out of disk
bool inDisk(pt a, pt b, pt p) {
return dot(a-p, b-p) <= 0;
}
Math insight
# » # »
In fact, we can notice that P A · P B is equal to the power of point P
with respect to the circle of diameter [AB]: if O is the center of that
# » # »
circle and r its radius, then P A · P B = |OP |2 − r2 . This makes it
perfect for our purpose.
B
A
bool onSegment(pt a, pt b, pt p) {
return orient(a,b,p) == 0 && inDisk(a,b,p);
}
Finding the precise intersection between two segments [AB] and [CD] is
quite tricky: many configurations are possible and the intersection itself
might be empty, a single point or a whole segment.
54
To simplify things, we will separate the problem in two distinct cases:
1. Segments [AB] and [CD] intersect properly, that is, their intersection
is one single point which is not an endpoint of either segment. This is
easy to test with orient().
2. In all other cases, the intersection, if it exists, is determined by the
endpoints. If it is a single point, it must be one of A, B, C, D, and if
it is a whole segment, it will necessarily start and end with points in
A, B, C, D.
Let’s deal with the first case: there is a single proper intersection point
I. To test this, it suffices to test that A and B are on either side of line CD,
and that C and D are on either side of line AB. If the test is positive, we
find I as a weighted average of A and B.4
C
B
A D
proper intersection
Then to deal with the second case, we will test for every point among
A, B, C, D if it is on the other segment. If it is, we add it to a set S. Clearly,
an endpoint cannot be in the middle of the intersection segment, so S will
always contain 0, 1 or 2 distinct points, describing an empty intersection, a
single intersection point or an intersection segment.
4
We can understand the formula as the center of gravity of point A with weight |oB | and
point B with weight |oA |, which gives us a point I on [AB] such that |IA|/|IB| = oA /oB .
55
A
B
B
A B
C D
D A
C D C
S=∅ S = {A} S = {B, C}
no intersection intersection point intersection segment
set<pt,cmpX> inters(pt a, pt b, pt c, pt d) {
pt out;
if (properInter(a,b,c,d,out)) return {out};
set<pt,cmpX> s;
if (onSegment(c,d,a)) s.insert(a);
if (onSegment(c,d,b)) s.insert(b);
if (onSegment(a,b,c)) s.insert(c);
if (onSegment(a,b,d)) s.insert(d);
return s;
}
To find the distance between segment [AB] and point P , there are two cases:
either the closest point to P on [AB] is strictly between A and B, or it is
one of the endpoints (A or B). The first case happens when the orthogonal
projection of P onto AB is between A and B.
P P
P
B B B
A A A
closest to A closest to projection closest to B
56
To check this, we can use the cmpProj() method in line.
double segPoint(pt a, pt b, pt p) {
if (a != b) {
line l(a,b);
if (l.cmpProj(a,p) && l.cmpProj(p,b)) // if closest to
projection
return l.dist(p); // output distance to
line
}
return min(abs(p-a), abs(p-b)); // otherwise distance to A or B
}
We can find the distance between two segments [AB] and [CD] based on
the segment-point distance if we separate into the same two cases as for
segment-segment intersection:
1. Segments [AB] and [CD] intersect properly, in which case the distance
is of course 0.
2. In all other cases, the shortest distance between the segments is at-
tained in at least one of the endpoints, so we only need to test the four
endpoints and report the minimum.
This can be readily implemented with the functions at our disposal.
double segSeg(pt a, pt b, pt c, pt d) {
pt dummy;
if (properInter(a,b,c,d,dummy))
return 0;
return min({segPoint(a,b,c), segPoint(a,b,d),
segPoint(c,d,a), segPoint(c,d,b)});
}
C A A
B B B D
D
A D C C
proper intersection between points between point and
line
57
2.6 Polygons
In this section we will discuss basic tasks on polygons: how to find their
area and two ways to detect if a point is inside or outside them.
B
θ
We know that the area of this triangle is 12 |AB||AC| sin θ, because |AC| sin θ
is the length of the height coming down from C. This looks a lot like the
definition of cross product: in fact,
1 1 # » # »
|AB||AC| sin θ = AC × AC
2 2
Since O is the origin, it can be implemented simply like this:
double areaTriangle(pt a, pt b, pt c) {
return abs(cross(b-a, c-a)) / 2.0;
}
Now that we can compute the area of a triangle, the intuitive way to
find the area of a polygon would be to
1. divide the polygon into triangles;
2. add up all the areas.
However, it turns out that reliably dividing a polygon into triangles is a
difficult problem in itself. So instead we’ll add and subtract triangle areas
in a clever way. Let’s take this quadrilateral as an example:
D C
A B
58
Let’s take an arbitrary reference point O. Let’s consider the vertices of
ABCD in order, and for every pair of consecutive points P1 , P2 , we’ll add
# »
the area of OP1 P2 to the total if P1 P2 goes counter-clockwise around O, and
subtract it otherwise. Additions are marked in blue and subtractions in red.
We can see that this will indeed compute the area of quadrilateral
ABCD. In fact, it works for any polygon (draw a few more examples to
convince yourself).
Note that the sign (add or subtract) that we take for the area of OP1 P2
is exactly the sign that the cross product takes. If we take the origin as
reference point O, it gives this simple implementation:
double areaPolygon(vector<pt> p) {
double area = 0.0;
for (int i = 0, n = p.size(); i < n; i++) {
area += cross(p[i], p[(i+1)%n]); // wrap back to 0 if i == n
-1
}
return abs(area) / 2.0;
}
We have to take the absolute value in case the vertices are given in clockwise
order. In fact, testing the sign of area is a good way to know whether the
vertices are in counter-clockwise (positive) or clockwise (negative) order. It
is good practice to always put your polygons in counter-clockwise order,
by reversing the array of vertices if necessary, because some algorithms on
polygons use this property.
59
A2
A1
However, sometimes this can go wrong if the ray touches a vertex of the
polygon, as below. The ray from A3 intersects the polygon twice, but A3
is inside. We can try to solve the issue by counting one intersection per
segment touched, which would give three intersections for A3 , but then the
ray from A4 will intersect the polygon twice even though A4 is inside.
A4
A3
Below we show for some segments whether they are considered to cross
the ray or not. We can see in the last two examples that the behavior is
different if the segment touches the ray from below or from above.
60
Pi+1 Pi+1
A A
Pi Pi
• touches ray: OK • touches ray: KO
• halves 6=: OK • halves 6=: OK
⇒ crossing ⇒ no crossing
Pi+1
Pi+1
A A Pi
Pi
• touches ray: OK • touches ray: OK
• halves 6=: OK • halves 6=: KO
⇒ crossing ⇒ no crossing
Exercise 2
Verify that, with this new definition of crossing, A3 and A4 are cor-
rectly detected to be inside the polygon.
Checking the halves to which the points belong is easy, but checking that
the segment touches the ray is a bit more tricky. We could check whether
the segments [Pi , Pi+1 ] and [AB] intersect for B very far on the ray, but
it actually we can do it more simply using orient: if Pi is below and Pi+1
above, then orient(A, Pi , Pi+1 ) should be positive, and otherwise it should
be negative. We can then implement this with the code below:
// true if P at least as high as A (blue part)
bool above(pt a, pt p) {
return p.y >= a.y;
}
// check if [PQ] crosses ray from A
bool crossesRay(pt a, pt p, pt q) {
return (above(a,q) - above(a,p)) * orient(a,p,q) > 0;
}
61
int numCrossings = 0;
for (int i = 0, n = p.size(); i < n; i++) {
if (onSegment(p[i], p[(i+1)%n], a))
return !strict;
numCrossings += crossesRay(a, p[i], p[(i+1)%n]);
}
return numCrossings & 1; // inside if odd number of crossings
}
A2
A3
A5
A4
A6
A1
Points A1 and A2 are completely out of the curve so the winding number
around them is 0 (no turn). Points A3 and A4 are inside the main loop, which
goes counterclockwise, so the winding number around them is 1. The curve
turns twice counterclockwise around A5 , so the winding number is 2. Finally
the curve goes clockwise around A6 , for a winding number of −1.
62
Math insight
In fact, we can move the curve continuously without changing the
winding number as long as we don’t touch the reference point. There-
fore we can “untie” loops which don’t contain the point. That’s why,
when looking at A3 or A4 , we can completely ignore the loops that
contain A5 and A6 .
Math insight
If we move the reference point while keeping the curve unchanged,
the value of the winding number will only change when it crosses
the curve. If it crosses the curve from the right (according to its
orientation), the winding number increases by 1, and if it crosses it
from the left, the winding number decreases by 1.
+1 −1
Exercise 3
What value will areaPolygon() (section 2.6.1) give when applied to a
closed polyline that crosses itself, like the curve above, instead of a
simple polygon? Assume we don’t take the absolute value.
[Go to solution]
63
Then we simply sum it all up and figure out how many turns were made:
int windingNumber(vector<pt> p, pt a) {
double ampli = 0;
for (int i = 0, n = p.size(); i < n; i++)
ampli += angleTravelled(a, p[i], p[(i+1)%n]);
return round(ampli / (2*M_PI));
}
Warning
The winding number is not defined if the reference point is on the
curve/polyline. If it is the case, this code will give arbitrary results,
and potentially (int)NAN.
While the code above works, its use of floating-point numbers makes it non
ideal, and when coordinates are integers, we can do better. We will define a
new way to work with angles, as a type angle. This type will also be useful
for other tasks, such as for sweep angle algorithms.
Instead of working with amplitudes directly, we will represent angles by a
point and a certain number of full turns.5 More precisely, in this case, we will
use point (x, y) and number of turns t to represent angle atan2(y, x) + 2πt.
We start by defining the new type angle. We also define a utility function
t360() which turns an angle by a full turn.
struct angle {
pt d; int t = 0; // direction and number of full turns
angle t180(); // to be defined later
angle t360() {return {d, t+1};}
};
The range of angles which have the same value for t is (−π + 2πt, π + 2πt].
We will now define a comparator between angles. The approach is the
same as what we did for the polar sort in section 2.3.2, so we will reuse the
function half() which separates the plane into two halves so that angles
within one half are easily comparable:
bool half(pt p) {
return p.y > 0 || (p.y == 0 && p.x < 0);
}
5
This approach is based on an original idea in [2], see “Angle.h”.
64
It returns true for the part highlighted in blue and false otherwise.
Thus, in practice, it allows us to separate each range (−π + 2πt, π + 2πt]
into the subranges (−π + 2πt, 2πt], for which half() returns false, and
(2πt, π + 2πt], for which half() returns true.
(2πt, π + 2πt]
We can now write the comparator between angles, which is nearly iden-
tical to the one we used for polar sort, except that we first check the number
of full turns t.
bool operator<(angle a, angle b) {
return make_tuple(a.t, half(a.d), 0) <
make_tuple(b.t, half(b.d), cross(a.d,b.d));
}
We also define the function t180() which turns an angle by half a turn
counterclockwise. The resulting angle has an opposite direction. To find the
number of full turns t, there are two cases:
• if half(d) is false, we are in the lower half (−π + 2πt, 2πt], and we
will move to the upper half (2πt, π + 2πt], without changing t;
• if half(d) is true, we are in the upper half (2πt, π + 2πt], and we will
move to (−π + 2π(t + 1), 2π(t + 1)], the lower half for t + 1.
angle t180() {return {d*(-1), t + half(d)};}
We will now implement the function that will allow us to compute the
winding number. Consider an angle with direction point D. Given a new
direction D0 , we would like to move the angle in such a way that if direction
D0 is to the left of D, the angle increases, and if D0 is to the right of D, the
angle decreases.
D
D
O
D0
O
D0
D0left of D
D0 right of D
⇒ increase
⇒ decrease
65
angle moveTo(angle a, pt newD) {
// check that segment [DD’] doesn’t go through the origin
assert(!onSegment(a.d, newD, {0,0}));
We know that b as it is first defined is less than a full turn away from a, so
the two conditions are enough to bring it within half a turn of a.
We can use this to implement a new version of windingNumber() very
simply. We start at some vertex of the polygon, move vertex to vertex while
maintaining the angle, then read the number of full turns once we come back
to it.
int windingNumber(vector<pt> p, pt a) {
angle a{p.back()}; // start at last vertex
for (pt d : p)
a = moveTo(a, d); // move to first vertex, second, etc.
return a.t;
}
2.7 Circles
(x − x0 )2 + (y − y0 )2 = r2
66
r sin θ
θ
r cos θ
2.7.2 Circumcircle
The circumcircle of a triangle ABC is the circle that passes through all
three points A, B and C.
C
O
B
The radius can then be found by taking the distance to any of the three
points, or directly taking the length of perp(b*sq(c)- c*sq(b))/cross(b,c)
# »
/2, which represents vector AO.
67
Math insight
The formula can be easily interpreted as the intersection point of the
line segment bisectors of segments [AB] and [AC], when considering
coordinates relative to A. Consider the bisector of [AB]. It is per-
# »
pendicular to [AB], so its direction vector is perp(AB), and it passes
# »
through the middle point 21 AB, so its constant term (variable c in
# » # »
the line structure) is perp(AB) × 21 AB = − 12 |AB|2 . Similarly, the
# »
bisector of [AC] is defined by direction vector perp(AC) and constant
term − 21 |AC|2 . We then just plug those into the formula for line
intersection found in section 2.4.6:
# » # »
# » (− 21 |AB|2 ) perp(AC) − (− 21 |AC|2 ) perp(AB)
AO = # » # »
perp(AB) × perp(AC)
# » # »
perp(|AC|2 AB − |AB|2 AC)
= # » # »
2 AB × AC
Let’s assume there are two intersection points I and J. We first find the
midpoint of [IJ]. This happens to be the projection of O onto the line l,
which we will call P .
r
h
O P
d
68
Once we have found P , to find I and J we need to move√along the line
by a certain distance h. By the Pythagorean theorem, h = r2 − d2 where
d is the distance from O to l.
This gives the following implementation (note that we have to divide
by k v#»l k so that we move by the correct distance). It returns the number
of intersections, and places them in out. If there is only one intersection,
out.first and out.second are equal.
Similarly to the previous section, two circles (O1 , r1 ) and (O2 , r2 ) can have
either 0, 1, 2 or an infinity of intersection points (in case the circles are
identical).
0 intersections 1 intersection
2 intersections ∞ intersections
69
J
r1 r2
h
O1 O2
P
d2 + r12 − r22
cos(]O2 O1 J) =
2dr1
and since O1 P J is a right triangle,
d2 + r12 − r22
|O1 P | = r1 cos(]O2 O1 J) =
2d
which allows us to find P .
Now to find h = |P I| = |P J|,
p we apply the Pythagorean theorem on
triangle O1 P J, which gives h = r12 − |O1 P |2 .
This gives the following implementation, which works in a very similar
way to the code in the previous section. It aborts if the circles are identical.
int circleCircle(pt o1, double r1, pt o2, double r2, pair<pt,pt> &
out) {
pt d=o2-o1; double d2=sq(d);
if (d2 == 0) {assert(r1 != r2); return 0;} // concentric circles
double pd = (d2 + r1*r1 - r2*r2)/2; // = |O_1P| * d
double h2 = r1*r1 - pd*pd/d2; // = hˆ2
if (h2 >= 0) {
pt p = o1 + d*pd/d2, h = perp(d)*sqrt(h2/d2);
out = {p-h, p+h};
}
return 1 + sgn(h2);
}
70
Math insight
Let’s check that if d 6= 0 and variable h2 in the code is nonnegative,
there are indeed 1 or 2 intersections (the opposite is clearly true: if
h2 is negative, the length h cannot exist). The value of h2 is
d + r1 ≥ r2 d + r2 ≥ r1 r1 + r2 ≥ d
Here we will try and find a line which is tangent to two circles (O1 , r1 )
and (O2 , r2 ). There are two types of such tangents: outer tangents, for
which both circles are on the same side of the line, and inner tangents, for
which the circles are on either side.
71
outer tangents inner tangents
We will study the case of outer tangents. Our first goal is to find a unit
vector parallel to the rays [O1 P1 ] and [O2 P2 ], in other words, we want to
# »
find #»
v = O1 P1 /r1 . To do this we will try to find angle α marked on the
figure.
P1
v#»
P2
r2
h
r1 −
α
O1 d O2
If we project O2 onto the ray from O1 , this forms a right triangle with
hypothenuse d = |O1 O2 | and adjacent side r1 − r2 , whichp means cos α =
r1 −r2
d . By the Pythagorean theorem, the third side is h = d2 − (r1 − r2 )2 ,
h
and we can compute sin α = d .
# » # »
From this we find #»
v in terms of O1 O2 and perp O1 O2 as
# » # »
#»
v = cos α O1 O2 /d ± sin α perp O1 O2 /d
# » # »
(r1 − r2 ) O1 O2 ± h perp O1 O2
=
d2
where the ± depends on which of the two outer tangents we want to find.
We can then compute P1 and P2 as
P1 = O1 + r1 #»
v and P2 = O2 + r2 #»
v
72
Exercise 4
Study the case of the inner tangents and show that it corresponds
exactly to the case of the outer tangents if r2 is replaced by −r2 .
This will allow us to write a function that handles both cases at once
with an additional argument bool inner and this line:
if (inner) r2 = -r2;
This gives the following code. It returns the number of tangents of the
specified type. Besides,
• if there are 2 tangents, it fills out with two pairs of points: the pairs
of tangency points on each circle (P1 , P2 ), for each of the tangents;
• if there is 1 tangent, the circles are tangent to each other at some point
P , out just contains P 4 times, and the tangent line can be found as
line(o1,p).perpThrough(p) (see 2.4.3);
• if there are 0 tangents, it does nothing;
• if the circles are identical, it aborts.
int tangents(pt o1, double r1, pt o2, double r2, bool inner, vector<
pair<pt,pt>> &out) {
if (inner) r2 = -r2;
pt d = o2-o1;
double dr = r1-r2, d2 = sq(d), h2 = d2-dr*dr;
if (d2 == 0 || h2 < 0) {assert(h2 != 0); return 0;}
for (double sign : {-1,1}) {
pt v = (d*dr + perp(d)*sqrt(h2)*sign)/d2;
out.push_back({o1 + v*r1, o2 + v*r2});
}
return 1 + (h2 > 0);
}
Conveniently, the same code can be used to find the tangent to a circle
passing through a point by setting r2 to 0 (in which case the value of inner
doesn’t matter).
73
Chapter 3
3D geometry
In this chapter we cover the basic objects and concepts of 3D geometry, and
the operations we can do on them. Although the additional dimension adds
some complexity which can make things more difficult and interesting, a
lot of things work the same way as they do in 2D. Therefore, we will make
frequent references to chapter 2, which we consider as a prerequisite.
In this first section, we define our point representation, we explain how the
dot and cross products we saw in 2D translate into 3D, and we show how a
combination of the two, the mixed product, can help us define a 3D analog
of the orient() function.
We will define points and vectors in space by their coordinates (x, y, z): their
positions along three perpendicular axes.
P = (1, 3, 2)
y
x
74
struct p3 {
T x,y,z;
// Some comparators
bool operator==(p3 p) {return tie(x,y,z) == tie(p.x,p.y,p.z);}
bool operator!=(p3 p) {return !operator==(p);}
};
Choosing the scalar type T is done the same way as in 2D, see 2.1.2 for our
remarks on that.
#»
For convenience, we also define a zero vector 0 = (0, 0, 0):
p3 zero{0,0,0};
#» #»
w
w θ θ
θ
#»
v #»
v #»
w #»
v
θ < π/2 θ = π/2 θ > π/2
#»
v ·w#» = 5 #»
v ·w#» = 0 #»
v ·w#» = −5
75
• linearity: (v#»1 + v#»2 ) · w
#» = v#» · w
1
#» + v#» · w
2
#» (also on the right);
• #» #» = 0 iff #»
v ·w v and w #» are perpendicular;
• #» #» stays constant if w
v ·w #» moves perpendicular to #» v.
Like in 2D, dot product is very simple to implement:
T operator|(p3 v, p3 w) {return v.x*w.x + v.y*w.y + v.z*w.z;}
Since 3D geometry uses dot and cross product a lot, to shorten notations
we will be using operator | for dot product and operator * for cross product.
We chose them for these mostly arbitrary reasons:
• | has a lower precedence than * which is desirable, it kind of looks like
a “parallel” operator, and since it most often has to be parenthesized
(e.g. (v|w) == 0), it can be a bit reminiscent of the inner product
notation hv, wi;
• * is the closest thing to a cross in overloadable R++ operators, and
the compiler doesn’t produce a warning if you don’t paranthesize b*c
in expressions such as a|b*c, which we will use a lot.
We define the usual sq() and abs() based on dot product and add a
unit() function that makes the norm of a vector equal to 1 while preserving
its direction.
T sq(p3 v) {return v|v;}
double abs(p3 v) {return sqrt(sq(v));}
p3 unit(p3 v) {return v/abs(v);}
Like in 2D, we can also use dot product to find the amplitude in [0, π]
of the angle between vectors #» #» see secton 2.3.1 for more details.
v and w,
double angle(p3 v, p3 w) {
double cosTheta = (v|w) / abs(v) / abs(w);
return acos(max(-1.0, min(1.0, cosTheta)));
}
76
#» #»
v ×w
#»
w
#»
v
The right-hand rule says this: if you take your right hand, align your
thumb with #» v and your extended index w, #» and fold your middle finger at
a 90◦ angle, then it will point in the direction of #» #» Another way to
v × w.
#» #»
express it is to say that if you draw v and w on a sheet of paper and look at
the sign of their 2D cross product, if it is positive #»
v ×w#» will point up from
#» #»
the sheet, and if it is negative v × w will point down through the sheet.
#» #»
v
w #» #»
v ×w
#»
v
#» #»
v ×w #»
w
points towards you points away from you
77
Note that the z-coordinate is the same expression as the 2D cross prod-
uct: indeed, it is the value of the 2D cross product if #» #» are projected
v and w
onto plane z = 0.
A very useful combination of dot product and cross product is the mixed
product. We define the mixed product of three vectors #»
u , #» #» as
v and w
( #»
u × #» #»
v)·w
#» #»
w #»
n n
#»v #»v
d d
#»u an #» an
f f u
neo ne
o #»
w
pla pla
( #»
u × #» #» > 0
v)·w ( #»
u × #» #» < 0
v)·w
R
S R R
S P Q
P Q P Q S
S above P QR S on P QR S below P QR
orient(P, Q, R, S) > 0 orient(P, Q, R, S) = 0 orient(P, Q, R, S) < 0
78
This orient() function has several very nice properties. First, it stays
the same if we swap any three arguments in a circular way: for exam-
ple let’s take P, Q, S, then orient(P, Q, R, S) = orient(Q, S, R, P ). On the
other hand, swapping any two arguments changes its sign: for example,
orient(P, Q, R, S) = − orient(P, S, R, Q).
Exercise 5
Show that those properties also apply to the mixed product. For
example, ( #»
u × #» #» = ( #»
v)·w #» · #»
v × w) u and ( #»
u × #» #» = −( #»
v)·w #» · #»
u × w) v.
[Go to solution]
Exercise 6
A convenient way to check whether two lines P Q and RS are skew is
# » # » # »
to check whether P Q × RS · P R 6= 0: in fact you can replace P R
by any vector going from P Q to RS.
Using the properties of dot product, cross product and orient(), prove
that # » # » # »
P Q × RS · P R = − orient(P, Q, R, S)
[Go to solution]
79
#» R
n #» #»
n R n
P P R
P
Q Q Q
3.2 Planes
In this section we will discuss how to represent planes, and many ways we
can use them. We will see that they play a very similar role to the role lines
play in 2D, and many of the operations we defined in section 2.4 have a
direct equivalent here. Since the explanations are nearly identical, we make
them a bit shorter here; please refer to section 2.4 if you want more details.
Planes are sets of points (x, y, z) which obey an equation of the form ax +
by + cz = d. Here, a, b, c determine the orientation of the plane, while d
determines its position relative to the origin.
ax + by + cz = d + 1
ax + by + cz = d
ax + by + cz = d − 1
Vector #»
n = (a, b, c) is perpendicular to the plane and is called a normal
of the plane. The equation can be rewritten as #» n · (x, y, z) = d: that is, the
plane is formed by all points whose dot product with #» n is equal to constant
d. This makes sense, because we know that the dot product doesn’t change
when one vector (here, the point (x, y, z)) moves perpendicularly to the
other (here, the normal #» n ).
80
#»
n
Here are some other ways to define a plane, and how to find #»n and d in
each case:
• if we know the normal #» n and a point P belonging to the plane: we
#»
can find d as n · P ;
• if we know a point P and two (non-parallel) vectors #»
v and w#» that are
parallel to the plane: we can define #»
n = #»
v ×w#» then find d as above;
• if we know 3 non-collinear points P, Q, R on the plane: we first find
# » # »
two vectors #»
v = P Q and w#» = P R that are parallel from the plane,
#»
then find n and d as above.
We implement this with the following structure and constructors:
struct plane {
p3 n; T d;
// From normal n and offset d
plane(p3 n, T d) : n(n), d(d) {}
// From normal n and point P
plane(p3 n, p3 p) : n(n), d(n|p) {}
// From three non-collinear points P,Q,R
plane(p3 p, p3 q, p3 r) : plane((q-p)*(r-p), p) {}
The first thing we are interested for a plane Π in is the value of ax+by+cz−d.
If it’s zero, the point is on Π, and otherwise it tells us which side of Π the
point lies. We name it side(), like its 2D line equivalent:
81
T side(p3 p) {return (n|p)-d;}
#»
n P #»
n #»
n
P
P
P above P on plane P below
sideΠ (P ) = 3 sideΠ (P ) = 0 sideΠ (P ) = −2
And just like the side() for 2D lines, we can get the distance from it, if
we compensate for the norm of #»n:
double dist(p3 p) {return abs(side(p))/abs(n);}
#»
t
To find the new value d0 , we use the same argument as for 2D lines:
#»
suppose point P is on the old plane, that is, #»
n · P = d. Then P + t should
be in the new plane:
#» #» #»
d0 = #»
n · (P + t ) = #»
n · P + #»
n · t = d + #»
n· t
82
3.2.4 Orthogonal projection and reflection
P
P
P0
P0
projection reflection
When we have a plane Π, sometimes we will want to know what are the
coordinates of a point in Π. That is, suppose we have a few points that we
know are coplanar, and we want to use some 2D algorithm on them. How
do we get the 2D coordinates of those points?
To do this, we need to chose an origin point O on Π and two vectors
#» #» #»
dx and dy indicating the desired x and y directions. Let’s first assume dx
#»
and dy are perpendicular and their norms are 1. Then, to find the x- and
83
y-coordinate of a point P on Π, we simply need to compute
# » #»
x = OP · dx
# » #»
y = OP · dy
#» #» #»
and if we also have vector dz = dx × dy perpendicular to the first two, we
can find the “height” of P respective to Π:
# » #»
z = OP · dz
#»
dz
P
#»
#» dy
dx
84
coords(p3 p, p3 q, p3 r, p3 s) :
o(p), dx(q-p), dy(r-p), dz(s-p) {}
Math insight
The second constructor allows us specify points P, Q, R, S and choose
#» # » #» # » #» # »
dx = P Q, dy = P R and dz = P S directly. This can be useful if we
#» #»
don’t care that the 2D coordinate system (dx , dy ) is not orthonormal
(perpendicular and of norm 1), because it allows us to keep using
integer coordinates.
#» #»
If dx and dy are not perpendicular or do not have a norm of 1, the
computed distances and angles will not be correct. But if we are only
interested in the relative positions of lines and points, and finding
the intersections of lines or segments, then everything works fine.
Computing the 2D convex hull of a set of points is an example of
such a problem, because it only requires that the sign of orient() is
correct.
3.3 Lines
In this section we will discuss how to represent lines in 3D, and some related
problems, including finding the line at the intersection of two planes, and
finding the intersection point of a plane and a line.
a1 x + b1 y + c1 z = d1
a2 x + b2 y + c2 z = d2
but this representation is not very convenient to work with, and there are
many pairs planes we could choose.
85
a1 x + b1 y + c1 z = d1
a2 x + b2 y + c2 z = d2
O
#»
d
#»
Note that here, d plays the same role as #»
v did for 2D lines (see sec-
tion 2.4.1).
If we are given two points P, Q on the line, we can set O = P and
#» # »
d = P Q, so we implement the structure like this:
struct line3d {
p3 d, o;
// From two points P, Q
line3d(p3 p, p3 q) : d(q-p), o(p) {}
// From two planes p1, p2 (requires T = double)
line3d(plane p1, plane p2); // will be implemented later
86
3.3.2 Distance from a line
#» # »
A point P is on a line l described by O, d if and only if OP is parallel to
#» #» # » #»
d . This is the case when d × OP = 0 .
#» # »
d × OP
# »
OP
P O #»
d
#» # »
More generally, this product d × OP also gives us information about the
distance to the line: the distance from l to P is equal to
#» # »
d × OP
#»
d
Just like we did with 2D lines in section 2.4.4, we can sort points according
to their position along a line l. To find out if a point P should come before
another point Q, we simply need to check whether
#» #»
d ·P < d ·Q
Let’s say we want to project P on a line l, that is find the closest point
to P on l. Our usual approach to projecting things, which is to move P
perpendicularly until it touches l, doesn’t work as well here: indeed, there
are many possible directions that are perpendicular to l.
87
P P
P0
P0
O O
#» #»
d d
projection reflection
Instead we will start from O and move along the line until we reach the
projection of P . We have seen in the previous section that taking the dot
#»
product with d tells us how far some point is along l. In fact, if we compute
#» # »
d · OP , this tells us the (signed) distance from O to the projection of P ,
#»
multiplied by k d k. So we can find the projection this way:
p3 proj(p3 p) {return o + d*(d|(p-o))/sq(d);}
Once we’ve found the projection P 0 , we can find the reflection P 00 easily,
# »
since it is twice as far in the same direction: we have P 00 = P 0 + P P 0 , which
becomes 2P 0 − P if we allow vector operations on points.
p3 refl(p3 p) {return proj(p)*2 - p;}
Let’s say we have a plane Π, represented by vector #»n and real d, and a line
#»
l, represented by point O and d . To find an intersection between them, we
#»
need to find a point O + k d that lies on Π, that is, such that
#»
#»
n · O +kd = d
#»
n
O
#»
d
88
which can be implemented directly:
p3 inter(plane p) {return o - d*p.side(o)/(p.n|d);}
#»
Note that this is undefined when #»
n · d = 0, that is, when Π and l are
parallel.2
Π2
Π1
O
v#»1 v#»2
2
We take a closer look at criteria for parallelism and perpendicularity in section 3.4.
89
Let’s verify that O is on Π1 , that is, n# »1 · O = n# »1 · (v#»1 + v#»2 ) = d1 . Since
v#»2 is perpendicular to n# »1 , clearly n# »1 · v#»2 = 0. What remains is
d1 #»
n# »1 · v#»1 =
#»
2 n# »2 × d · n# »1
d
d1 #»
=
#»
2 n# »1 × n# »2 · d
d
d1 #» #»
=
#»
2 d · d
d
= d1
where between the first two lines, we used the fact that the mixed product
is conserved if we swap its arguments in a circular way (see exercise 5).
We prove similarly that O is on Π2 . Finally we note that both v#»1 and
#» #»
v2 are perpendicular to d , so the vector from the origin to O is perpencular
#»
to d , the direction of the line. Therefore it must necessarily arrive on the
line at the closest point to the origin.
We can implement this as the following constructor:
line3d(plane p1, plane p2) {
d = p1.n*p2.n;
o = (p2.n*p1.d - p1.n*p2.d)*d/sq(d);
}
90
Let’s call C1 the point of l1 that is closest to l2 , and C2 the point on l2
# »
that is closest to l1 . Direction C1 C2 should be perpendicular to both l1 and
l2 . Indeed, if it were not, it would be possible to get a smaller distance by
# » #» #»
moving either C1 or C2 . So C1 C2 is parallel to #» n = d1 × d2 .
l2
#»
d2
l1
C2
#»
n
C1
#»
d1
We don’t know C1 or C2 yet, but since the dot product doesn’t change when
one of the vectors moves perpendicular to the other, we can move C1 to O1
and C2 to O2 , so that we get
# » #» # » #» # » #»
C1 C2 · n O1 C2 · n O1 O2 · n
#» = #» =
knk knk k #»
nk
Now, how can we find C1 and C2 themselves? Let’s call Π2 the plane
#»
that contains l2 and is parallel to n2 ; its normal is n# »2 = d2 × #»
n . Then C1 is
the intersection between that plane and line l1 , so we can use the formula
obtained in section 3.3.5 to find it:
# »
O1 O2 · n# »2
C1 = O1 + #» # »
d1 · n2
91
l2
#»
d2
Π2
C2
#»
n
C1
l1 n# »2
The angle between two planes is equal to the angle between their normals.
Since usually two angles of distinct amplitudes θ and π − θ are formed, we
take the smaller of the two, in [0, π2 ].
n# »2
n# »1 # »
n2
θ
π−θ
We can find it with the following code. We take the minimum with 1 to
avoid nan in case of imprecisions.
92
double smallAngle(p3 v, p3 w) {
return acos(min(abs(v|w)/abs(v)/abs(w), 1.0));
}
double angle(plane p1, plane p2) {
return smallAngle(p1.n, p2.n);
}
The situation with lines is exactly the same: their angle is equal to the angle
between their direction vectors. Note that the lines aren’t necessarily in the
same plane, so the angle is taken as if they were moved until they touch.
θ
#» #»
d1 d2
π−θ
93
3.4.3 Between a plane and a line
The situation when considering a plane and a line is a bit different. Let’s
#»
consider a plane Π of normal #» n and a line l of direction vector d . When
#»
they are perpendicular, #» n is parallel to d , and inversely when they are
#» #»
parallel, n is perpendicular to d . In general, if the angle between #»
#» n and d
is θ ∈ [0, π2 ], then the angle between the plane and the line is π2 − θ.
θ
#»
#»
n d
O + #»
n
#»
n O
#»
The plane perpendicular to a line l of direction vector d and going
#»
through a point O is simply the plane containing O and whose normal is d .
94
#»
d
O
3.5 Polyhedrons
3.5.1 Definition
y
x
95
Note that those conditions don’t exclude nonconvex polyhedrons (like
the example above), but they do exclude self-crossing polyhedrons.
z
P4
P1 P3
y
x P2
How do we find the area of a face P1 · · · Pn ? First let’s take the most
# » # »
simple case: a triangle ABC. If we compute cross product AB × AC, its
direction is perpendicular to the triangle and its norm is
|AB| |AC| sin θ
# » # »
where θ is the amplitude of angle between AB and AC. This is twice the
area of triangle ABC,
# »as was already noted in section 2.6.1. So the area of
# »
triangle ABC is 21
AB × AC
.
# » # »
AB × AC
C
A
B
96
• Triangles ABC and ACD are oriented in different ways (the angle at
# » # » # » # »
either B or D is concave). In this case, vectors AB × AC and AC × AD
point in opposite directions, and we should take the difference of the
# » # » # » # »
areas. So again, taking the vector sum AB × AC + AC × AD will
produce the desired effect.
D
D
C
A A B
B C
# » # » # » # » # » # » # » # »
AB × AC AC × AD AB × AC AC × AD
If we continue with more and more vertices, we can make the following
general conclusion. Given a polygon P1 · · · Pn , we can compute the vector
area
#» P1 × P2 + P2 × P3 + · · · + Pn × P1
S =
2
#»
This vector is perpendicular to the polygon, and
S
gives its area. It
is oriented such that when it points towards the observer, the vertices are
numbered in counter-clockwise order.
#»
S
P5
P4
P1 P3
P2
Note that this is very similar to the 2D formula for area we obtained in
section 2.6.1. The only thing that changes is that the result is a vector. We
can implement this straightforwardly:
97
p3 vectorArea2(vector<p3> p) { // vector area * 2 (to avoid
divisions)
p3 S = zero;
for (int i = 0, n = p.size(); i < n; i++)
S = S + p[i]*p[(i+1)%n];
return S;
}
double area(vector<p3> p) {
return abs(vectorArea2(p)) / 2.0;
}
#»
S pointing outside (lengths not to scale)
While it’s not clear when looking at individual faces which side is inside
the polyhedron, it’s easy to deduce the correct orientation by looking at the
orientation of an adjacent face. If two faces share an edge [P Q], and the
first face lists P then Q in this order, then the other face should list them
in the other order, so that they “rotate” in the same direction. Note that
because of circularity, in P1 · · · Pn , Pn is considered to come before P1 , not
after.
98
Q Q
R R
S S
P P
99
int v = es[{b,a}];
g[u].push_back({v,false});
g[v].push_back({u,false});
} else { // not seen yet
es[{a,b}] = u;
}
}
}
3.5.4 Volume
Suppose that the faces are oriented correctly. We will show how to compute
the volume of the polyhedron by taking the same approach as in section 2.6.1
for the 2D polygon area.
Let’s choose an arbitrary reference point O. We will compute the volume
of the polyhedron face by face: for a face P1 · · · Pn , if the side of the face
seen from O is inside the polygon, add the volume of pyramid OP1 · · · Pn ,
100
and otherwise subtract it. That way, by inclusion and exclusion, the final
result will be the volume inside the polyhedron.
#»
Since S always points towards the outside of the polygon, we just have
#»
to check whether S points away from O (add) or towards O (subtract).
#»
S
P2 P3
P1 P1
P3 P2
#»
S
O O
#» #»
S away from O S towards O
+ volume of OP1 P2 P3 − volume of OP1 P2 P3
#»
S
P2
P1
P3
|OP1 | cos θ
θ
O
#» # »
So the absolute value of S · OP1 is equal to the volume of pyramid
#»
OP1 · · · Pn , and the dot product is positive if S points away from O, and
#»
negative if S points towards O. This is exactly the sign we want.
For the implementation, we take O to be the origin for convenience. We
divide by 6 at the end because we need to divide by 2 to get the correct area
and then by 3 because of the formula for the volume of a pyramid.
101
double volume(vector<vector<p3>> fs) {
double vol6 = 0.0;
for (vector<p3> f : fs)
vol6 += (vectorArea2(f)|f[0]);
return abs(vol6) / 6.0;
}
#»
In case the vector areas S all point towards the inside of the polyhedron
(which may happen after applying the face orientation procedure in the
previous section), vol6 will be correct but will have a negative sign. If this
#»
happens, flip all the faces so that all S now point outside.
(x − x0 )2 + (y − y0 )2 + (z − z0 )2 = r2
r
r
y
x
102
us how far North the point is (or South, if ϕ < 0), and an angle λ, the
longitude, which tells us how far East the point is (or West, if λ < 0). We
usually take ϕ ∈ [− π2 , π2 ] and λ ∈ (−π, π].
North
West
ϕ East
λ
x
South
103
r
h
O P
d
This function returns the number of intersection points and places them
in pair out if they exist.
int sphereLine(p3 o, double r, line3d l, pair<p3,p3> &out) {
double h2 = r*r - l.sqDist(o);
if (h2 < 0) return 0; // the line doesn’t touch the sphere
p3 p = l.proj(o); // point P
p3 h = l.d*sqrt(h2)/abs(l.d); // vector parallel to l, of length
h
out = {p-h, p+h};
return 1 + (h2 > 0);
}
B
O θ
104
}
This code also works if A and B are not actually on the sphere, in which
case it will give the distance between their projections on the sphere:
Greenland
London
Los Angeles
105
Given two spherical segments [AB] and [CD], we would like to figure
out if they intersect, and what their intersection is. This is part of a more
general problem: given two segments [AB] and [CD] in space, if we view
them from an observation point O, does one of them hide part of the other,
that is, is there a ray from O that touches them both?
B
B C
C D
D O
A
A
common point on sphere common ray from O
We will solve the general problem, and make sure that our answer is
always exact when points A, B, C, D are integer points. To do this, we will
separate cases just like we did for 2D segment intersection in section 2.5.2:
1. Segments [AB] and [CD] intersect properly, that is, their intersection
is one single point which is not an endpoint of either segment. For
the general problem this means that there is a single ray from O that
touches both [AB] and [CD], and it doesn’t touch any of A, B, C, D.
2. In all other cases, the intersection, if it exists, is determined by the
endpoints. If it is a single point, it must be one of A, B, C, D, and if
it is a whole segment, it will necessarily start and end with points in
A, B, C, D.
For the rest of explanation, we will consider the case of spherical segment
intersection, because it’s easier to visualize, but it should be easy to verify
that this also applise to the general problem.
Proper intersection
Let’s deal with the first case: there is a single proper intersection point I.
For this to be the case, A and B must be on either side of plane OCD, and
C and D must be on either side of plane OAB. Put another way, A and
B must be on either side of the great circle containing C and D, and vice
versa.
106
B
C D
D
C
A
To eliminate this kind of case, we also have to check that oA and oC have
opposite signs (this is clearly not the case here).
Exercise 7
Consider a few more examples and verify that these criteria correctly
detect proper intersections in all cases.
107
ob = sgn(cd|b),
oc = sgn(ab|c),
od = sgn(ab|d);
out = ab*cd*od; // four multiplications => careful with overflow
!
return (oa != ob && oc != od && oa != oc);
}
Improper intersections
To deal with the second case, we will do as for 2D segments and test for
every point among A, B, C, D if it is on the other segment. If it is, we add
it to a set S. S will contain 0, 1, or 2 distinct points, describing an empty
intersection, a single intersection point or an intersection segment.
C
D D
D
C C
A B A B A B
108
#»
n #»
n #»
n
P O O O
P
A B A P B A B
There remains only one special case: if A and B are the same point on
#»
the sphere, then #»
n = 0 , and then we should just check that P is also that
same point.
We arrive at the following implementation. To handle the general prob-
lem, instead of directly checking for equality between P and A or B, we
check that they are in the same direction with the cross product.
bool onSphSegment(p3 a, p3 b, p3 p) {
p3 n = a*b;
if (n == zero)
return a*p == zero && (a|p) > 0;
return (n|p) == 0 && (n|a*p) >= 0 && (n|b*p) <= 0;
}
Now we just have to put all of this together in one function. First we
check for a proper intersection, then if there is none we check segment by
segment and add the points to set S. Since (as mentioned) we can’t check
for equality directly, we use a custom set structure that checks if the cross
product is zero for every point already in the set.
struct directionSet : vector<p3> {
using vector::vector; // import constructors
void insert(p3 p) {
for (p3 q : *this) if (p*q == zero) return;
push_back(p);
}
};
directionSet intersSph(p3 a, p3 b, p3 c, p3 d) {
assert(validSegment(a, b) && validSegment(c, d));
p3 out;
if (properInter(a, b, c, d, out)) return {out};
directionSet s;
if (onSphSegment(c, d, a)) s.insert(a);
109
if (onSphSegment(c, d, b)) s.insert(b);
if (onSphSegment(a, b, c)) s.insert(c);
if (onSphSegment(a, b, d)) s.insert(d);
return s;
}
Given two spherical segments [AB] and [AC] on a sphere around the origin
O, how do we compute the amplitude of the angle they form on the surface
of the sphere at A? This angle is equal to the angle between planes OAB
and OAC, or more precisely between their normals A × B and A × C. Thus
we can find the angle using angle(), which gives values in [0, π].
A×B A×C
O
C C
A θ B A θ B
double angleSph(p3 a, p3 b, p3 c) {
return angle(a*b, a*c);
}
If instead of values in [0, π], we want to know the oriented angle3 between
[AB] and [AC], that is, how much we rotate if we go from B to C around A
counterclockwise, then we need to know on which side of plane OAB point
C lies. If C lies “to the left” of plane OAB, the angle given by angle() is
correct, but if C lies “to the right”, we should subtract it from 2π.
C B
A θ B A θ C
≈ 30◦
angleSph(a,b,c) angleSph(a,b,c)≈ 30◦
orientedAngleSph(a,b,c) ≈ 30◦ orientedAngleSph(a,b,c) ≈ 330◦
3
We defined a similar notion in 2D in section 2.3.2.
110
double orientedAngleSph(p3 a, p3 b, p3 c) {
if ((a*b|c) >= 0)
return angleSph(a, b, c);
else
return 2*M_PI - angleSph(a, b, c);
}
D D
C C
A A
B B
111
SA
C
SC S
B
A
SB
Since those four triangles cover half of the sphere together, we have
S + SA + SB + SC = 2πr2 . Besides if we take triangle ABC and the triangle
on the other side of [BC], together they form a whole spherical wedge 4 of
angle α. So their combined area should be α/2π of the total area, that is
S + SA = 2αr2 . Similarly, we obtain S + SB = 2βr2 and S + SC = 2γr2 .
Combining those four equations, we obtain
2S = (S + SA ) + (S + SB ) + (S + SC ) − (S + SA + SB + SC )
= r2 (2α + 2β + 2γ − 2π)
112
Earth, even though their actual sizes are very different. So we say that they
subtend the same solid angle at the observation point that is Earth.
Let’s define it more precisely. Just like the planar angle subtended by
an object is the length of the object once it is projected onto a unit circle
around the point, the solid angle subtended by an object is the area of the
object once it is projected onto a unit sphere around the point.
O O
The unit for solid angles is the steradian (sr), and because the area of
a unit sphere is 4π, a solid angle of 4π means that the observation point is
completely surrounded, while a solid angle of 2π means that half of the view
is covered.
We can easily find the solid angle subtended by a polygon by using the
function areaOnSphere() we just defined and setting r = 1. Indeed, all it
does is compute angles, so the distance from the origin O doesn’t matter.
113
Math insight
#»
The solid angle subtended by a small surface d S at a position #» r is
inversely proportional to the square of the distance from the origin,
#»
k #»
r k2 , and proportional to the cosine of the angle between #»
r and d S ,
because a surface seen from sideways occupies less of the view. So we
can also find solid angles with the following integral:
Z #» #»
r · dS
Ω=
k #»
r k3
#»
#» S
S
114
To find what quantity to add or subtract for each face, we will use
#»
function areaOnSphere() directly. If S points away from O, it will return a
value in (0, 2π), the area of the projection of the face a unit sphere, which
#»
we should keep. If S points towards O, it will return a value in (2π, 4π),
the area of the rest of the unit sphere, so we should remove 4π to get the
subtraction we want.
Since we always want the value to be in (−2π, 2π) in the end, we can
use function remainder(), giving the simple implementation below.
int windingNumber3D(vector<vector<p3>> fs) {
double sum = 0;
for (vector<p3> f : fs)
sum += remainder(areaOnSphere(1, f), 4*M_PI);
return round(sum / (4*M_PI));
}
115
Appendix A
Exercise 1
Prove that (r1 cis ϕ1 ) ∗ (r2 cis ϕ2 ) = (r1 r2 ) cis(ϕ1 + ϕ2 ) using this new
definition of product.
[Back to exercise]
Exercise 3
What value will areaPolygon() (section 2.6.1) give when applied to a
closed polyline that crosses itself, like the curve above, instead of a
simple polygon? Assume we don’t take the absolute value.
It will give the sum of the areas of the parts delimited by the curve, multi-
plied by their corresponding winding numbers. For the curve below, it will
give the sum of:
• 1 × the area of the part containing A3 and A4 ;
• 2 × the area of the part containing A5 ;
• −1 × the area of the part containing A6 .
116
A2
A3
A5
A4
A6
A1
[Back to exercise]
Exercise 5
Show that those properties also apply to the mixed product. For
example, ( #»
u × #» #» = ( #»
v)·w #» · #»
v × w) u and ( #»
u × #» #» = −( #»
v)·w #» · #»
u × w) v.
[Back to exercise]
Exercise 6
A convenient way to check whether two lines P Q and RS are skew is
# » # » # »
to check whether P Q × RS · P R 6= 0: in fact you can replace P R
by any vector going from P Q to RS.
Using the properties of dot product, cross product and orient(), prove
that # » # » # »
P Q × RS · P R = − orient(P, Q, R, S)
117
We develop:
# » # » # » # » # » # » # »
P Q × RS · P R = P Q × P S − P R ·PR
# » # » # » # » # »
= PQ × PS − PQ × PR ·PR
# » # » # » # » # » # »
= PQ × PS ·PR − PQ × PR ·PR
# » # » # »
= PQ × PS ·PR
= orient(P, Q, S, R)
= − orient(P, Q, R, S)
[Back to exercise]
118
Appendix B
Omitted proofs
119
Lemma 1. Let f (x) = (1 + )x − 1. If a, b > 0, then f (a) + f (b) ≤ f (a + b).
We find
| round(a0 ∗ b0 ) − (a ∗ b)|
= |(round(a0 ∗ b0 ) − (a0 ∗ b0 )) + ((a0 ∗ b0 ) − (a ∗ b))|
≤ | round(a0 ∗ b0 ) − (a0 ∗ b0 )| + |(a0 − a) ∗ (b0 − b)|
≤ M d + |a0 − a| + |b0 − b|
≤ M d + M d ((1 + )na − 1) + M d ((1 + )nb − 1)
= M d [f (1) + f (na ) + f (nb )]
≤ M d f (na + nb + 1)
= M d (1 + )na +nb +1 − 1
where the step before last follows from two applications of Lemma 1.
120
Proof. From the hypotheses we know that
|a| ≤ M d
|b| ≤ M d
|a0 − a| ≤ M da ((1 + )na − 1)
|b0 − b| ≤ M da ((1 + )nb − 1)
|a0 b0 | ≤ M da +db .
We find
| round(a0 b0 ) − ab)|
= |(round(a0 b0 ) − a0 b0 ) + (a0 b0 − ab)|
≤ | round(a0 b0 ) − a0 b0 | + |(a0 − a)b + (b0 − b)a + (a0 − a)(b0 − b)|
≤ M da +db + |a0 − a||b| + |b0 − b||a| + |a0 − a||b0 − b|
≤ M da +db + M da ((1 + )na − 1) M db + M db ((1 + )nb − 1) M da
+ M da ((1 + )na − 1) M db ((1 + )nb − 1)
= M da +db + ((1 + )na − 1) + ((1 + )nb − 1)
+ ((1 + )na − 1) ((1 + )nb − 1)
= M da +db + (1 + )na +nb − 1
= M da +db [f (1) + f (na + nb )]
≤ M da +db f (na + nb + 1)
= M da +db (1 + )na +nb +1 − 1
121
Bibliography
[2] Simon Lindholm et al. KTH ACM Contest Template Library. 2017.
url: https://github.com/kth-competitive-programming/kactl.
122