C++

Points, Vectors and the Orientation Predicate

The one question almost every geometry algorithm asks, "which side of the line is this point on?", answered with a sign instead of an angle.

This lesson is a free sample of the course.

One question under every geometry algorithm

A convex hull wraps a set of points in the smallest convex shape. A triangulation cuts a polygon into triangles. A sweep line finds out whether any two segments cross. These look like different problems, and their textbook descriptions are full of different ideas. Underneath, they keep asking the same small question:

Given a directed line from a to b and a third point c, does c lie to the left of the line, to the right, or on it?

Call that the orientation of the triple (a, b, c). It is the same as asking whether walking a, then b, then c makes a left turn, a right turn, or goes straight on. A convex hull algorithm keeps a point only if it makes a turn in the right direction. A triangulation keeps an ear only if its corner turns the right way. The sweep line asks whether the endpoints of one segment fall on opposite sides of another.

If that one function is right, a large part of geometry is right. If it is subtly wrong, algorithms built on it fail in ways that are hard to trace: a hull that has a dent in it, a triangulation with a sliver of negative area, an infinite loop on input that "should" be fine. That is why this course begins with a single function and spends a whole lesson on it.

This lesson builds the function, shows how to read it, and shows the two ways people usually get it wrong. The next lesson, Floating-Point Robustness, deals with what happens when the coordinates are not whole numbers.

Points and vectors in a dozen lines

A point is a position. A vector is a displacement: the difference between two positions. In code they are both two numbers, and it helps to keep them as separate types in your head even when you store them in one struct.

The two operations on vectors that matter here are the dot product and the 2D cross product. For vectors u = (ux, uy) and v = (vx, vy):

  • dot(u, v) = ux*vx + uy*vy. It is zero when the vectors are perpendicular, positive when they point roughly the same way, negative when they point roughly opposite ways.
  • cross(u, v) = ux*vy - uy*vx. It is the z-component of the 3D cross product of the two vectors laid in the xy-plane. It is zero when they are parallel, positive when v is counter-clockwise from u, negative when it is clockwise.
main.cppc++23
#include <print>

struct Vec { double x, y; };

Vec operator-(Vec a, Vec b) { return {a.x - b.x, a.y - b.y}; }
double dot(Vec a, Vec b) { return a.x * b.x + a.y * b.y; }
double cross(Vec a, Vec b) { return a.x * b.y - a.y * b.x; }

int main() {
    Vec a{1, 2}, b{4, 6};
    Vec ab = b - a;  // the vector that takes you from a to b
    std::println("ab = ({}, {})", ab.x, ab.y);
    std::println("dot(ab, ab) = {}", dot(ab, ab));  // 3*3 + 4*4

    Vec east{1, 0}, north{0, 1};
    std::println("dot(east, north)   = {}", dot(east, north));    // perpendicular
    // north is counter-clockwise from east; swap the arguments and the sign flips
    std::println("cross(east, north) = {}", cross(east, north));
    std::println("cross(north, east) = {}", cross(north, east));
}
Output

Press Run to compile and run this example.

ab = (3, 4)
dot(ab, ab) = 25
dot(east, north)   = 0
cross(east, north) = 1
cross(north, east) = -1

Notice the last two lines. The cross product is antisymmetric: swapping its arguments flips the sign. That single fact is why one number can answer a question with three possible answers.

The sign of a cross product is the turn

Take the three points a, b, c. The vector b - a is the direction of travel. The vector c - a points from the start to the candidate. Their cross product,

cross(b - a, c - a) = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x),

is positive when c - a is counter-clockwise from b - a, which means c is on the left of the directed line. It is negative when c is on the right, and zero when the three points lie on one line.

The magnitude means something too. |cross(b - a, c - a)| is the area of the parallelogram spanned by the two vectors, so half of it is the area of the triangle abc. The signed value is twice the signed area of the triangle, positive when a, b, c run counter-clockwise. You will use this again for polygon area in lesson 5.

A function that returns the raw number is useful, but most callers want only the sign. A small enum makes the call sites read like the question being asked:

main.cppc++23
#include <print>

struct Point { double x, y; };

// Twice the signed area of triangle abc.
double cross3(Point a, Point b, Point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

enum class Turn { Right = -1, Collinear = 0, Left = 1 };

Turn orientation(Point a, Point b, Point c) {
    double v = cross3(a, b, c);
    return v > 0 ? Turn::Left : v < 0 ? Turn::Right : Turn::Collinear;
}

const char* name(Turn t) {
    switch (t) {
        case Turn::Left: return "left";
        case Turn::Right: return "right";
        case Turn::Collinear: return "collinear";
    }
    return "?";
}

int main() {
    Point a{0, 0}, b{4, 0};  // the line points east
    for (Point c : {Point{2, 3}, Point{2, -3}, Point{9, 0}, Point{-1, 0}}) {
        Turn t = orientation(a, b, c);
        double v = cross3(a, b, c);
        std::println("c = ({}, {}): {:<9} (cross = {})", c.x, c.y, name(t), v);
    }
}
Output

Press Run to compile and run this example.

c = (2, 3): left      (cross = 12)
c = (2, -3): right     (cross = -12)
c = (9, 0): collinear (cross = 0)
c = (-1, 0): collinear (cross = 0)

Two of those answers deserve a second look. Both (9, 0) and (-1, 0) are "collinear", even though one is beyond b and the other is behind a. Orientation says only that the three points are on one line, not that c is on the segment from a to b. The segment question needs a second test, which comes later in this lesson.

Seeing it

Numbers about left and right are easier to trust once you have watched them sort a field of points. This program places a grid of points around a directed line and colors each one by its orientation: blue for left, red for right, green for exactly on the line.

SVG's y-axis points down, while the mathematics above has it pointing up. The program computes everything in the usual "y up" coordinates and flips only when it writes the picture.

main.cppc++23
#include <print>

struct Point { double x, y; };

double cross3(Point a, Point b, Point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

constexpr double kHeight = 120;
double flip(double y) { return kHeight - y; }  // SVG's y axis points down

void label(const char* text, double x, double y) {
    std::println(R"(<text x="{}" y="{}" font-size="9">{}</text>)", x, flip(y), text);
}

int main() {
    Point a{10, 90}, b{190, 30};
    std::println(R"(<svg viewBox="0 0 200 120">)");
    std::println(R"(<line x1="{}" y1="{}" x2="{}" y2="{}" )"
                 R"(stroke="#1f2937" stroke-width="2"/>)",
                 a.x, flip(a.y), b.x, flip(b.y));
    for (double x = 10; x < 200; x += 20) {
        for (double y = 10; y < 120; y += 20) {
            double v = cross3(a, b, {x, y});
            const char* color = v > 0 ? "#2563eb" : v < 0 ? "#dc2626" : "#16a34a";
            std::println(R"(<circle cx="{}" cy="{}" r="3.5" fill="{}"/>)",
                         x, flip(y), color);
        }
    }
    label("a", a.x - 2, a.y - 12);
    label("b", b.x - 2, b.y + 6);
    std::println("</svg>");
}
Output

Press Run to compile and run this example.

The line runs from a at the upper left down to b at the lower right. Walking from a to b, the points above the line are on your left-hand side, which is why they are blue. The green points lie exactly on the line; because a and b were chosen on the grid, four grid points land on it exactly.

Why not just use angles

Most people's first idea is trigonometry: compute the angle of a→b, compute the angle of a→c, and compare them. It feels natural, and it is wrong in a way that the cross product is not.

Angles live on a circle, so their difference wraps around. If c is straight behind a, the two directions are exactly opposite, the difference is ±π, and a comparison against zero reports a turn even though the three points are collinear. You can patch that with special cases, and each patch adds another place to be wrong. The angle version also calls atan2 twice, which costs far more than the two multiplications and five subtractions of the cross product.

This program counts how often the angle version disagrees with the exact answer over every ordered triple of distinct-direction points in a 5 by 5 grid of integers.

main.cppc++23
#include <cmath>
#include <numbers>
#include <print>
#include <vector>

struct Point { int x, y; };

// Exact: the sign of a cross product in 64-bit integers.
int exactTurn(Point a, Point b, Point c) {
    long long v = 1LL * (b.x - a.x) * (c.y - a.y) - 1LL * (b.y - a.y) * (c.x - a.x);
    return (v > 0) - (v < 0);
}

// The tempting version: compare the directions a->b and a->c as angles.
int angleTurn(Point a, Point b, Point c) {
    double ab = std::atan2(b.y - a.y, b.x - a.x);
    double ac = std::atan2(c.y - a.y, c.x - a.x);
    double d = ac - ab;
    const double pi = std::numbers::pi;
    if (d > pi) d -= 2 * pi;
    if (d <= -pi) d += 2 * pi;
    const double eps = 1e-9;
    return (d > eps) - (d < -eps);
}

bool same(Point p, Point q) { return p.x == q.x && p.y == q.y; }

int main() {
    std::vector<Point> grid;
    for (int x = 0; x < 5; ++x)
        for (int y = 0; y < 5; ++y) grid.push_back({x, y});

    long triples = 0, wrong = 0, behind = 0;
    for (Point a : grid)
        for (Point b : grid)
            for (Point c : grid) {
                if (same(a, b) || same(a, c)) continue;  // a direction needs 2 points
                ++triples;
                bool bad = exactTurn(a, b, c) != angleTurn(a, b, c);
                wrong += bad;
                long long dot = 1LL * (b.x - a.x) * (c.x - a.x)
                              + 1LL * (b.y - a.y) * (c.y - a.y);
                behind += bad && exactTurn(a, b, c) == 0 && dot < 0;
            }
    std::println("triples checked:            {}", triples);
    std::println("misclassified by angles:    {}", wrong);
    std::println("  of which c is behind a:   {}", behind);
}
Output

Press Run to compile and run this example.

triples checked:            14400
misclassified by angles:    304
  of which c is behind a:   304

Every one of the 304 failures is a collinear triple with c behind a. The cross product has no such case: it never asks for an angle, so there is nothing to wrap.

Which side of a segment

The orientation test answers questions about a line. Many algorithms need the answer for a segment, where "on the line" splits into three cases: between the endpoints, beyond b, or behind a. When the cross product is zero, the dot product separates them.

This version uses integer coordinates and long long. That is exact as long as the coordinates stay within a bound that you can state: if every coordinate satisfies |x| <= 10^9, each difference is at most 2*10^9, each product at most 4*10^18, and the difference of two products at most 8*10^18, which is below the largest long long (about 9.22*10^18). The static_assert records that reasoning in the code, where the next reader will see it.

main.cppc++23
#include <limits>
#include <print>

struct Point { long long x, y; };

constexpr long long kMaxCoord = 1'000'000'000;
static_assert(8 * kMaxCoord * kMaxCoord < std::numeric_limits<long long>::max(),
              "cross product of coordinates up to kMaxCoord must not overflow");

long long cross3(Point a, Point b, Point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
long long dot3(Point a, Point b, Point c) {  // (b - a) . (c - a)
    return (b.x - a.x) * (c.x - a.x) + (b.y - a.y) * (c.y - a.y);
}

enum class Place { Left, Right, Behind, OnSegment, Beyond };

Place classify(Point a, Point b, Point c) {
    long long turn = cross3(a, b, c);
    if (turn > 0) return Place::Left;
    if (turn < 0) return Place::Right;
    if (dot3(a, b, c) < 0) return Place::Behind;   // on the line, before a
    if (dot3(b, a, c) < 0) return Place::Beyond;   // on the line, past b
    return Place::OnSegment;
}

const char* name(Place p) {
    switch (p) {
        case Place::Left: return "left of the line";
        case Place::Right: return "right of the line";
        case Place::Behind: return "on the line, behind a";
        case Place::OnSegment: return "on the segment";
        case Place::Beyond: return "on the line, beyond b";
    }
    return "?";
}

int main() {
    Point a{0, 0}, b{6, 3};
    const Point candidates[] = {{2, 5}, {4, -1}, {-2, -1}, {4, 2},
                                {8, 4}, {6, 3}, {0, 0}};
    for (Point c : candidates) {
        std::println("({:>2}, {:>2}) is {}", c.x, c.y, name(classify(a, b, c)));
    }
}
Output

Press Run to compile and run this example.

( 2,  5) is left of the line
( 4, -1) is right of the line
(-2, -1) is on the line, behind a
( 4,  2) is on the segment
( 8,  4) is on the line, beyond b
( 6,  3) is on the segment
( 0,  0) is on the segment

The endpoints themselves count as on the segment. Whether that is the right choice depends on the algorithm: a point-in-polygon test, for example, has to decide what to do when a point is exactly on an edge. Make the decision once, in one function, and write it down.

A first look at robustness

All of the above is exact because the coordinates are integers. With double coordinates the same formula can return the wrong sign, and the wrong sign is exactly what breaks hulls and triangulations.

Here is the sign function written once for any coordinate type, applied to three points that are exactly collinear on paper (they lie on the line y = 3x):

main.cppc++23
#include <print>

template <class T>
struct Point2 { T x, y; };

template <class T>
int orientation(Point2<T> a, Point2<T> b, Point2<T> c) {
    T v = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
    return (v > T{}) - (v < T{});  // -1, 0 or +1
}

int main() {
    // The same three points, written as integers and as decimals.
    int exact = orientation<long long>({1, 3}, {2, 6}, {3, 9});
    int approx = orientation<double>({0.1, 0.3}, {0.2, 0.6}, {0.3, 0.9});
    std::println("integers: {}", exact);
    std::println("doubles:  {}", approx);
}
Output

Press Run to compile and run this example.

integers: 0
doubles:  1

The doubles disagree with arithmetic. The value 0.1 has no exact binary representation, so the subtractions leave tiny errors, and the cross product comes out around 2e-17 instead of 0. The sign function faithfully reports "left". The next lesson is about what to do about that: when an epsilon helps, when it hurts, and how to get a predicate whose sign is always right.

For now the practical rule is simple. If your coordinates can be integers, make them integers. Scale decimal data to a fixed grid (millimeters instead of meters) and use long long. You get exact predicates and no tolerance to tune.

A screen is upside down

Everything in this lesson assumes the y-axis points up, as in a math textbook. Most screens and image formats put the origin at the top-left with y growing downward. Flip the axis and every orientation flips with it: a triple that turns counter-clockwise on paper turns clockwise on the screen, and Left becomes what a user sees as "right".

This is a frequent source of "my convex hull is inside out" bugs, because nothing fails. The function is correct; the picture is mirrored.

main.cppc++23
#include <print>

struct Point { double x, y; };

double cross3(Point a, Point b, Point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

int main() {
    // The corners of a triangle, listed counter-clockwise in maths coordinates.
    Point p{0, 0}, q{4, 0}, r{0, 3};
    auto turn = [](double v) { return v > 0 ? "counter-clockwise" : "clockwise"; };
    std::println("maths coordinates (y up):    {}", turn(cross3(p, q, r)));

    // The same triangle in screen coordinates: only the y values change sign.
    Point ps{0, 0}, qs{4, -0.0}, rs{0, -3};
    std::println("screen coordinates (y down): {}", turn(cross3(ps, qs, rs)));
}
Output

Press Run to compile and run this example.

maths coordinates (y up):    counter-clockwise
screen coordinates (y down): clockwise
Pitfall

Three mistakes account for most orientation bugs, and none of them is about the formula.

  • Swapping the argument order. orientation(a, b, c) is not orientation(a, c, b): the two have opposite signs. Write the call so that it reads as "walk from a to b, then look at c".
  • Comparing a floating-point cross product to 0 with ==, or to a fixed epsilon. == treats 1e-17 as a turn. A fixed epsilon such as 1e-9 is wrong for coordinates near 1e6 (where it is far too small) and for coordinates near 1e-6 (where it is far too large). Lesson 2 covers the correct approach.
  • Mixing coordinate systems. Computing in y-up and reading the result in y-down, or the reverse, mirrors your answer. Convert at the boundary and say which system each function expects.
Best Practice

Write the orientation test once, name it, and call it everywhere. Do not inline (b.x - a.x) * (c.y - a.y) - ... in a dozen algorithms.

  • A single function is the only place to fix when you decide to change number types, add exact arithmetic or add a tolerance. Dozens of inlined copies will not all be fixed.
  • Return a sign (-1, 0, +1) or an enum, not the raw value. Callers then cannot accidentally depend on the magnitude, and swapping in an exact implementation later does not change their code.
  • Template the coordinate type, as in the last-but-one program. Algorithms written against orientation<T> work with long long today and a better type tomorrow.
  • Test it with the awkward inputs first: all three points equal, two equal, collinear in each of the three possible orders, and coordinates at the largest magnitude you promise to support.

Try it yourself

  1. In the program that classifies points around a segment, add the point (3, 1). Predict first which of the five places it falls in, then run it. Then compute the cross product by hand.
  2. In the drawing program, change b to {190, 31}. Predict how many green points remain. Run it. What does that tell you about how rare exact collinearity is in real data?
  3. In the angle comparison, change eps from 1e-9 to 0. Which of the three counters changes, and why?
  4. In the last program, write the screen triangle with rs{0, -3} replaced by rs{0, 3}. Is it clockwise or counter-clockwise now, and what does that say about which corner order is "the same triangle"?
  5. Make the grid in the angle comparison 7 by 7. Predict what happens to the fraction of misclassified triples before you run it.
Note

CGAL, the Computational Geometry Algorithms Library, exposes this exact test as CGAL::orientation(p, q, r), and its algorithms (convex hull, triangulation, Boolean operations on polygons) are written in terms of it. Lesson 2 shows why CGAL takes so much care over the implementation of that one function that it lets you choose between exact and filtered arithmetic.

Computer graphics rasterizers decide whether a pixel is inside a triangle with the same sign test, applied to each of the triangle's three edges. A pixel is inside when all three agree.

Note
  • Mark de Berg, Otfried Cheong, Marc van Kreveld and Mark Overmars, "Computational Geometry: Algorithms and Applications", 3rd edition, Springer, 2008. The standard textbook; its first chapter builds the convex hull on this test.
  • Jonathan Richard Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", Discrete & Computational Geometry 18(3), 1997. The paper behind the robust orientation test used by many mesh generators.
  • CGAL documentation, for the library's own description of its predicates and kernels.
Key Takeaways
  • Orientation asks whether c is left of, right of, or on the directed line a→b. Convex hulls, triangulations and sweeps are all built from it.
  • cross(b - a, c - a) gives the answer in its sign and twice the triangle's area in its magnitude. Swapping two arguments flips the sign.
  • Angles are the wrong tool: differences wrap around, so opposite directions are reported as a turn, and atan2 is slow.
  • Zero cross product means "same line", not "on the segment". Use the dot product to place a collinear point.
  • With integer coordinates up to 10^9 in absolute value, long long gives exact, overflow-free orientation. State such bounds in a static_assert.
  • With double, exactly collinear points can produce a tiny non-zero cross product. Prefer integer coordinates, and see the next lesson for the rest.
  • Write the predicate once, return a sign, and template the coordinate type.

Quiz — Test Your Knowledge

(40 XP)

1. What does the sign of cross(b - a, c - a) tell you?

2. You compute orientation on screen coordinates, where y grows downward, using the formula from this lesson. A triangle you listed counter-clockwise on paper gives a negative cross product. What happened?

3. What is |cross(b - a, c - a)| / 2?

4. Why does comparing atan2(c - a) with atan2(b - a) fail to detect collinear points?

5. All coordinates are integers with absolute value at most 10^9. Which type computes the orientation determinant exactly with plain arithmetic?

Next
Floating-Point Robustness: Epsilons, Exact Predicates and Filters