Geometric Primitives & Orientation

This lesson opens the computational geometry module, where the objects are points, segments, and polygons in the plane rather than numbers or graphs. The algorithms ahead, including convex hulls, sweep-line intersection, and closest pairs, look involved, but they nearly all rest on one small operation applied repeatedly: given three points, does the path through them turn left or right? If that primitive is exact, the rest is bookkeeping; if it is subtly wrong, every structure built on top inherits the error.

The central design decision of this lesson is therefore exact arithmetic. The natural geometric quantities (angles, lengths, slopes) are irrational and force floating point, where a quantity that should be zero comes out as $10^{-16}$ and a collinearity test flips the wrong way. We avoid them. With integer input coordinates, the orientation and area primitives below are polynomials in the coordinates, so they evaluate to exact integers and their signs are never in doubt.¹ Slopes, square roots, and $\arctan$ never appear.

Points as vectors

We identify a point $P=(x,y)$ with the vector from the origin to it, which lets us do arithmetic on geometry. For points $A=(a_x,a_y)$ and $B=(b_x,b_y)$ and a scalar $t$:

The subtraction $B-A$ is the most important of these: it is the displacement vector pointing from $A$ to $B$, and almost every primitive below is phrased in terms of such difference vectors anchored at a common point. If $A$, $B$ have integer coordinates, so does $B-A$, and exactness is preserved by every one of these operations.

from __future__ import annotations

from typing import NamedTuple

class Point(NamedTuple):
  """
    A planar point / vector with x and y components.\n
    Integer coordinates keep all derived primitives exact.\n
  """
  x: int
  y: int

  def add(self, other: Point) -> Point:
    """
      Componentwise vector sum.\n
      A named method rather than `+` so the type does not have to override\n
      tuple's concatenating `__add__`, whose return type is incompatible\n
      with a fixed two-component Point.\n
    """
    return Point(self.x + other.x, self.y + other.y)

  def __sub__(self, other: Point) -> Point:
    """
      Displacement from `other` to this point — the vector other -> self.\n
    """
    return Point(self.x - other.x, self.y - other.y)

  def scale(self, factor: int) -> Point:
    """
      Scalar multiple of the vector.\n
    """
    return Point(self.x * factor, self.y * factor)

The dot product: angle and projection

The dot product of two vectors measures how much they point the same way:

where $\theta$ is the angle between them. The second equality is the useful one in reverse: because $|\vec{a}|,|\vec{b}|>0$, the sign of the dot product is the sign of $\cos \theta$, so it classifies the angle without ever computing it.

Three uses recur. Projection: the scalar projection of $\vec{b}$ onto $\vec{a}$ is $(\vec{a}\cdot \vec{b})/|\vec{a}|$, the signed length of $\vec{b}$'s shadow along $\vec{a}$. Perpendicularity: $\vec{a}\perp \vec{b}\leftrightarrow \vec{a}\cdot \vec{b}=0$, an exact integer test. Angle: $\cos \theta =(\vec{a}\cdot \vec{b})/(|\vec{a}||\vec{b}|)$ when the actual angle is genuinely needed (the one place a square root sneaks in). The dot product is symmetric and says nothing about which side one vector lies on; for that we need the cross product.

import math

from vector import Point

def dot(first: Point, second: Point) -> int:
  """
    The dot product a . b = a_x b_x + a_y b_y, exact for integer vectors.\n
  """
  return first.x * second.x + first.y * second.y

def angle_sign(first: Point, second: Point) -> int:
  """
    Sign of the angle between the vectors: +1 acute, 0 right, -1 obtuse.\n
    This is the sign of cos(theta), read straight off the dot product.\n
  """
  product: int = dot(first, second)
  if product > 0:
    return 1
  if product < 0:
    return -1
  return 0

def is_perpendicular(first: Point, second: Point) -> bool:
  """
    Whether the two vectors meet at a right angle (a . b == 0).\n
  """
  return dot(first, second) == 0

def projection_length(vector: Point, onto: Point) -> float:
  """
    Signed length of `vector`'s shadow along `onto`: (a . b) / |a|.\n
    Negative when `vector` points opposite `onto`. The lone square root in\n
    the lesson, so this returns a float rather than an exact integer.\n
  """
  magnitude: float = math.hypot(onto.x, onto.y)
  if magnitude == 0.0:
    raise ValueError("cannot project onto the zero vector")
  return dot(vector, onto) / magnitude

The cross product: signed area and orientation

In the plane the cross product of two vectors is a single scalar:

Its magnitude $|\vec{a}\times \vec{b}|$ equals the area of the parallelogram spanned by $\vec{a}$ and $\vec{b}$ (and twice the area of the triangle they form). Its sign is the sign of $\sin \theta$, which encodes orientation: positive when $\vec{b}$ lies counterclockwise from $\vec{a}$, negative when clockwise, zero when the two are parallel (collinear). Unlike the dot product, the cross product is antisymmetric: $\vec{a}\times \vec{b}=-(\vec{b}\times \vec{a})$. This single quantity, an exact integer for integer inputs, underlies nearly everything that follows.

For example, take $\vec{a}=(4,1)$ and $\vec{b}=(1,3)$. Then

The value $11$ is the area of the parallelogram those two vectors span, so the triangle $O$–$\vec{a}$–$\vec{b}$ has area $\frac{11}{2}$; the positive sign says $\vec{b}$ sits counterclockwise from $\vec{a}$, and swapping the arguments flips the sign but not the magnitude. Contrast the dot product on the same vectors, $\vec{a}\cdot \vec{b}=4\cdot 1+1\cdot 3=7>0$, which only reports that the angle between them is acute — it cannot tell counterclockwise from clockwise. The two products are complementary: dot for how aligned, cross for which side and how much area.

from fractions import Fraction

from vector import Point

def cross(first: Point, second: Point) -> int:
  """
    The scalar cross product a x b = a_x b_y - a_y b_x, exact for integers.\n
  """
  return first.x * second.y - first.y * second.x

def parallelogram_area(first: Point, second: Point) -> int:
  """
    Area of the parallelogram spanned by the two vectors: |a x b|.\n
  """
  return abs(cross(first, second))

def triangle_area(corner: Point, first: Point, second: Point) -> Fraction:
  """
    Unsigned area of the triangle with vertices `corner`, `first`, `second`.\n
    Half the cross product of the two edges anchored at `corner`; an exact\n
    half-integer, so a Fraction keeps it exact.\n
  """
  signed: int = cross(first - corner, second - corner)
  return Fraction(abs(signed), 2)

The orientation test

Anchor two difference vectors at a common point $A$ and take their cross product. This is the orientation (or ccw) test, the fundamental primitive of planar computational geometry:

It reports the sense of the turn made by the directed path $A\to B\to C$:

The three outcomes are best seen on concrete integer coordinates. Fix $A=(0,0)$ and $B=(4,0)$ and slide the third point: $C$ above the line $AB$ gives a left turn, $C$ below gives a right turn, and $C$ on the line gives $0$. Each value below is a single subtraction of products, exact to the last bit.

Three properties matter. The test uses only additions and multiplications of the input coordinates, so for integer inputs it is exact. It is antisymmetric in a way that respects the geometry: swapping $B$ and $C$ flips the sign, matching the reversal of the turn. And it answers, in $O(1)$, the question every higher geometric algorithm reduces to — which side of the line $AB$ does $C$ lie on? The collinearity test are $A$, $B$, $C$ on one line? is just $\mathrm{ccw}(A,B,C)=0$, with no division by a slope and hence no vertical-line special case.

To see the arithmetic end to end, fix $A=(1,1)$ and $B=(5,3)$ and test three different third points, evaluating $d=(b_x-a_x)(c_y-a_y)-(b_y-a_y)(c_x-a_x)$ with $(b_x-a_x,b_y-a_y)=(4,2)$ throughout:

The middle row is the important one: $(3,2)$ is the exact midpoint of $A$ and $B$, so it lies on the line, and the primitive returns a clean integer $0$ rather than a floating-point near-zero that a threshold test might misclassify. That exactness is why the collinearity test never needs an epsilon.

from cross_product import cross
from vector import Point

def orientation(first: Point, second: Point, third: Point) -> int:
  """
    Sign of the turn at first -> second -> third.\n
    +1 left turn (counterclockwise), -1 right turn (clockwise),\n
    0 collinear. Exact integer arithmetic.\n
  """
  turn: int = cross(second - first, third - first)
  if turn > 0:
    return 1
  if turn < 0:
    return -1
  return 0

def is_collinear(first: Point, second: Point, third: Point) -> bool:
  """
    Whether the three points lie on a single line — ccw == 0.\n
    No slope and hence no vertical-line special case.\n
  """
  return orientation(first, second, third) == 0

def is_left_turn(first: Point, second: Point, third: Point) -> bool:
  """
    Whether the path turns counterclockwise (strictly left) at `second`.\n
  """
  return orientation(first, second, third) > 0

def is_right_turn(first: Point, second: Point, third: Point) -> bool:
  """
    Whether the path turns clockwise (strictly right) at `second`.\n
  """
  return orientation(first, second, third) < 0

Segment intersection

When do two segments $\overline{AB}$ and $\overline{CD}$ cross? The slope-and-solve approach drags in division and degenerate cases; orientation makes it a handful of sign comparisons. Straddling drives it: segment $\overline{AB}$ straddles the line through $C$ and $D$ when its endpoints fall on opposite sides of that line, that is, when $\mathrm{ccw}(C,D,A)$ and $\mathrm{ccw}(C,D,B)$ have opposite signs.

The two ccw evaluations are exact integers, so this test has no rounding error and never computes the intersection point itself.

The boundary, when some ccw is $0$ and three points are collinear, needs care, and is where naive implementations break. The case logic:

• All four products nonzero (the generic case): intersect iff both products are strictly negative, as above.
• Exactly one $\mathrm{ccw}$ is $0$, say $\mathrm{ccw}(A,B,C)=0$: then $C$ lies on the line $AB$; the segments touch iff $C$ lies on the segment$\overline{AB}$, i.e. $C$'s coordinates are within the bounding box of $A$ and $B$ (an on-segment check: $\min (a_x,b_x)\le c_x\le \max (a_x,b_x)$ and likewise for $y$).
• The segments are collinear (all four ccw vanish): they overlap iff their $1$-D projections onto the $x$-axis (or $y$, if vertical) overlap, again a bounding-box / interval-overlap test, no geometry beyond comparisons.

Run the generic case on integers. Let $\overline{AB}$ go from $A=(0,0)$ to $B=(4,4)$ and $\overline{CD}$ from $C=(0,4)$ to $D=(4,0)$ — the two diagonals of a square, which plainly cross at $(2,2)$. The four orientations are

Both products $d_1{d}_2=-256<0$ and $d_3{d}_4=-256<0$ are negative, so each segment straddles the other's line and the test reports a proper crossing — all without ever solving for the point $(2,2)$. Now shorten the second segment to $C=(0,4)$, $D=(1,3)$, a stub that stops well short of the diagonal $AB$. Both $C$ and $D$ lie on the same side of line $AB$ (the line $y=x$): $d_3=\mathrm{ccw}(A,B,C)=+16>0$ and $d_4=\mathrm{ccw}(A,B,D)=(4)(3)-(4)(1)=+8>0$, so $d_3{d}_4>0$ and the test reports no crossing, again with no coordinates computed.

from orientation import orientation
from vector import Point

def on_segment(start: Point, end: Point, query: Point) -> bool:
  """
    Whether `query`, already known to be collinear with the segment\n
    start -> end, actually lies within it — a bounding-box check.\n
  """
  within_x: bool = (
    min(start.x, end.x) <= query.x <= max(start.x, end.x)
  )
  within_y: bool = (
    min(start.y, end.y) <= query.y <= max(start.y, end.y)
  )
  return within_x and within_y

def segments_intersect(
  first_start: Point,
  first_end: Point,
  second_start: Point,
  second_end: Point,
) -> bool:
  """
    Whether segments first_start->first_end and second_start->second_end\n
    share at least one point — proper crossings plus collinear/touching.\n
  """
  # straddle signs: each endpoint against the other segment's line.
  first_vs_second_start: int = orientation(
    second_start, second_end, first_start
  )
  first_vs_second_end: int = orientation(
    second_start, second_end, first_end
  )
  second_vs_first_start: int = orientation(
    first_start, first_end, second_start
  )
  second_vs_first_end: int = orientation(
    first_start, first_end, second_end
  )

  # generic case: each segment straddles the other's supporting line.
  straddles_first: bool = first_vs_second_start * first_vs_second_end < 0
  straddles_second: bool = second_vs_first_start * second_vs_first_end < 0
  if straddles_first and straddles_second:
    return True

  # boundary cases: a vanishing ccw means a point sits on the other line,
  # so it touches iff it lies within that segment's bounding box.
  if first_vs_second_start == 0 and on_segment(
    second_start, second_end, first_start
  ):
    return True
  if first_vs_second_end == 0 and on_segment(
    second_start, second_end, first_end
  ):
    return True
  if second_vs_first_start == 0 and on_segment(
    first_start, first_end, second_start
  ):
    return True
  if second_vs_first_end == 0 and on_segment(
    first_start, first_end, second_end
  ):
    return True

  return False

Polygon area: the shoelace formula

Given a polygon as an ordered list of vertices $P_0,P_1,\dots ,P_{n-1}$ (indices modulo $n$), its area is the shoelace formula:

Each term $x_i{y}_{i+1}-x_{i+1}{y}_i$ equals the cross product $P_i\times P_{i+1}$, so the area is a sum of cross products, halved.

Drop the absolute value and the sign of the sum carries orientation: positive means the vertices are listed counterclockwise, negative means clockwise. This is the cheapest way to detect the winding direction of a polygon, and because the sum is an integer for integer vertices, the area comes out as an exact half-integer.² The full derivation — the trapezoid sum that proves the formula — is in Polygons & Proximity.

For example, take the square listed counterclockwise, $P_0=(0,0)$, $P_1=(2,0)$, $P_2=(2,2)$, $P_3=(0,2)$. The four cross terms $x_i{y}_{i+1}-x_{i+1}{y}_i$ are $0\cdot 0-2\cdot 0=0$, $2\cdot 2-2\cdot 0=4$, $2\cdot 2-0\cdot 2=4$, and $0\cdot 0-0\cdot 2=0$, summing to $8$. Halved, the area is $4$ — the $2\times 2$ square, exactly. Reverse the vertex order and the sum is $-8$; the sign flip is the only difference, and it reports the winding direction.

from fractions import Fraction
from typing import Sequence

from cross_product import cross
from vector import Point

def twice_signed_area(polygon: Sequence[Point]) -> int:
  """
    Twice the signed area, sum of P_i x P_{i+1} around the boundary.\n
    Kept doubled so it stays an exact integer; positive iff the vertices\n
    wind counterclockwise.\n
  """
  # sum each edge's cross product with its successor (wrapping at the end).
  count: int = len(polygon)
  return sum(
    cross(polygon[index], polygon[(index + 1) % count])
    for index in range(count)
  )

def signed_area(polygon: Sequence[Point]) -> Fraction:
  """
    The signed area of the polygon, an exact half-integer.\n
    Positive when the vertices are counterclockwise, negative when clockwise.\n
  """
  return Fraction(twice_signed_area(polygon), 2)

def polygon_area(polygon: Sequence[Point]) -> Fraction:
  """
    The unsigned area enclosed by the polygon, an exact half-integer.\n
  """
  return abs(signed_area(polygon))

def is_counterclockwise(polygon: Sequence[Point]) -> bool:
  """
    Whether the vertices are listed in counterclockwise winding order.\n
  """
  return twice_signed_area(polygon) > 0

Point in polygon

Is a query point $q$ inside a polygon? Two exact strategies, both built from the primitives above. Ray casting shoots a ray from $q$ in a fixed direction (say $+x$) and counts how many polygon edges it crosses: an odd count means $q$ is inside, even means outside — the Jordan-curve parity argument. Each ray-vs-edge crossing is decided with the same straddle/orientation tests, with careful tie-breaking when the ray grazes a vertex (count an edge only if exactly one endpoint is strictly above the ray). The winding number alternative sums the signed angles the polygon's edges subtend at $q$ (computed from cross- and dot-product signs, not actual angles); a total winding of $0$ means outside, $\pm 1$ means inside, and unlike parity it stays correct for self-intersecting polygons. For a convex polygon both can be sped up to $O(\log n)$: binary-search the vertex fan around $P_0$ to find the wedge containing $q$ using orientation tests, then one final ccw against the bounding edge decides inside vs. outside.³

For a convex polygon the fan of diagonals from $P_0$ cuts the interior into triangular wedges in angular order, so a single binary search on $\mathrm{ccw}(P_0,P_i,q)$ locates the wedge $⟨P_k,P_{k+1}⟩$ that could contain $q$; one last orientation test against the edge $\overline{P_k{P}_{k+1}}$ decides inside vs. outside.

from typing import Sequence

from orientation import orientation
from segment_intersect import on_segment
from vector import Point

def _on_boundary(polygon: Sequence[Point], query: Point) -> bool:
  """
    Whether `query` lies exactly on some edge of the polygon.\n
  """
  count: int = len(polygon)
  for index in range(count):
    # an edge contains the query iff it is collinear and within the segment.
    start: Point = polygon[index]
    end: Point = polygon[(index + 1) % count]
    if orientation(start, end, query) == 0 and on_segment(start, end, query):
      return True

  return False

def in_polygon_ray_casting(polygon: Sequence[Point], query: Point) -> bool:
  """
    Inside-test by +x ray-crossing parity; boundary counts as inside.\n
    An edge is counted only when exactly one endpoint is strictly above the\n
    ray, which breaks ties cleanly when the ray grazes a shared vertex.\n
  """
  if _on_boundary(polygon, query):
    return True

  # toggle parity once per edge the +x ray crosses to its right.
  inside: bool = False
  count: int = len(polygon)

  for index in range(count):
    start: Point = polygon[index]
    end: Point = polygon[(index + 1) % count]

    # half-open vertical span: exactly one endpoint strictly above the ray.
    straddles_ray: bool = (start.y > query.y) != (end.y > query.y)
    if not straddles_ray:
      continue

    # x of the edge at the ray's height, compared without division by
    # cross-multiplying; the edge lies to the right of the point iff the
    # query is on the correct side of the directed edge.
    lower, upper = (start, end) if start.y < end.y else (end, start)
    side: int = orientation(lower, upper, query)
    if side > 0:
      inside = not inside
  return inside

def winding_number(polygon: Sequence[Point], query: Point) -> int:
  """
    The number of times the polygon winds counterclockwise around `query`.\n
    Zero means outside; a nonzero value means inside even for a\n
    self-intersecting boundary. Computed from ccw signs, not real angles.\n
  """
  winding: int = 0
  count: int = len(polygon)

  for index in range(count):
    start: Point = polygon[index]
    end: Point = polygon[(index + 1) % count]
    side: int = orientation(start, end, query)

    # upward edge passing left wraps +1; downward edge passing left wraps -1.
    if start.y <= query.y and end.y > query.y and side > 0:
      winding += 1
    elif start.y > query.y and end.y <= query.y and side < 0:
      winding -= 1

  return winding

def in_polygon_winding(polygon: Sequence[Point], query: Point) -> bool:
  """
    Inside-test by nonzero winding number; boundary counts as inside.\n
  """
  if _on_boundary(polygon, query):
    return True
  return winding_number(polygon, query) != 0

def in_convex_polygon(polygon: Sequence[Point], query: Point) -> bool:
  """
    O(log n) inside-test for a convex polygon given counterclockwise.\n
    Binary-search the fan of diagonals from P_0 to find the wedge that\n
    could hold `query`, then one ccw on that wedge's far edge decides it.\n
    Boundary points count as inside.\n
  """
  count: int = len(polygon)
  if count < 3:
    # degenerate: fall back to the collinear on-segment check.
    if count == 2:
      return on_segment(polygon[0], polygon[1], query) and (
        orientation(polygon[0], polygon[1], query) == 0
      )
    return count == 1 and polygon[0] == query

  anchor: Point = polygon[0]

  # the query must be left of the first fan ray and right of the last,
  # otherwise it is outside the angular span the fan covers.
  first_side: int = orientation(anchor, polygon[1], query)
  last_side: int = orientation(anchor, polygon[count - 1], query)
  if first_side < 0 or last_side > 0:
    # still allow points lying exactly on the two bounding fan rays.
    if first_side == 0:
      return on_segment(anchor, polygon[1], query)
    if last_side == 0:
      return on_segment(anchor, polygon[count - 1], query)
    return False

  # binary-search for the wedge <P_low, P_high> containing the query.
  low: int = 1
  high: int = count - 1
  while high - low > 1:
    middle: int = (low + high) // 2
    if orientation(anchor, polygon[middle], query) >= 0:
      low = middle
    else:
      high = middle

  # inside iff the query is left of (or on) the wedge's far edge.
  return orientation(polygon[low], polygon[high], query) >= 0

Robust predicates and adaptive precision

The lesson's exact-integer stance works whenever inputs are integers and the predicates stay low-degree — orientation is degree $2$, the in-circle test used by Delaunay triangulation is degree $4$. But real geometric software must handle floating-point inputs (coordinates from sensors, CAD, GIS), and there the orientation determinant $(b_x-a_x)(c_y-a_y)-(b_y-a_y)(c_x-a_x)$ can suffer catastrophic cancellation: when three points are nearly collinear the two products are nearly equal, and subtracting them in double can return a value whose sign is wrong. A wrong sign is worse than a small numerical error: it can make a hull algorithm loop forever, or a triangulation report a non-planar mesh. Integer arithmetic avoids this failure entirely.

The standard fix is adaptive-precision exact predicates, developed by Jonathan Shewchuk (1997).⁴ The idea reconciles speed with correctness. First evaluate the determinant in fast floating point together with an error bound on the roundoff; if the computed value exceeds that bound in magnitude, its sign is certified correct and we return immediately — the common case, at nearly the cost of the naive test. Only when the value falls inside the error bound (the points are close to collinear, exactly when the sign is in doubt) does the code fall back to slower exact arithmetic, and even then it computes just enough extra digits to resolve the sign, not the full exact value. This staged strategy is what makes industrial-strength libraries reliable.

The alternative philosophy, followed by CGAL's exact geometric computation paradigm, is to make every predicate exact from the start using number types that carry as much precision as needed, accepting a constant-factor slowdown for a guarantee of topological consistency.⁵ The choice between epsilon comparisons and exact predicates is not pedantry: downstream combinatorics depends on the orientation sign as a discrete decision, so a single flipped sign corrupts the whole structure. Keeping inputs integer, as we do here, is the simplest way to get exactness for free.

Takeaways

• Treat points as vectors; the difference $B-A$ is the displacement that almost every primitive is built from, and integer inputs stay integer.
• The dot product $\vec{a}\cdot \vec{b}=a_x{b}_x+a_y{b}_y$ has the sign of $\cos \theta$: positive/zero/negative $\Rightarrow$ acute/right/obtuse — it handles projection, perpendicularity, and angle.
• The cross product $\vec{a}\times \vec{b}=a_x{b}_y-a_y{b}_x$ gives signed parallelogram area in its magnitude and orientation in its sign.
• The orientation test $\mathrm{ccw}(A,B,C)=\mathrm{sign}((B-A)\times (C-A))$ reports left/right/collinear in $O(1)$ exact integer arithmetic — the primitive that hull, intersection, and point-location all reduce to.
• Segments cross iff each straddles the other's line (opposite ccw signs on both); collinear and touching cases fall to on-segment bounding-box checks.
• The shoelace formula $\frac{1}{2}\left|\sum P_i\times P_{i+1}\right|$ is a sum of cross products giving area, and its sign reveals CCW vs. CW winding.
• Point-in-polygon is ray-casting parity or winding number; a convex polygon admits an $O(\log n)$ orientation-based binary search.