xix. 3D Geometry and 3D Hull

3D Geometry

Primitives

template <class T>
struct Point3 {
    T x, y, z;
    Point3(T x_ = 0, T y_ = 0, T z_ = 0) : x(x_), y(y_), z(z_) {}
    template <class U> explicit Point3(const Point3<U>& p) : x(T(p.x)), y(T(p.y)), z(T(p.z)) {}
 
    Point3 operator+(const Point3& p) const { return Point3(x + p.x, y + p.y, z + p.z); }
    Point3 operator-(const Point3& p) const { return Point3(x - p.x, y - p.y, z - p.z); }
    Point3 operator*(T k) const { return Point3(x * k, y * k, z * k); }
    Point3 operator/(T k) const { return Point3(x / k, y / k, z / k); }
    Point3 operator-() const { return Point3(-x, -y, -z); }
    friend Point3 operator*(T k, const Point3& p) { return p * k; }
 
    bool operator==(const Point3& p) const { return x == p.x && y == p.y && z == p.z; }
    bool operator!=(const Point3& p) const { return !(*this == p); }
    bool operator<(const Point3& p) const { return tie(x, y, z) < tie(p.x, p.y, p.z); }
 
    T norm2() const { return x * x + y * y + z * z; }
    ld length() const { return sqrtl(ld(x) * ld(x) + ld(y) * ld(y) + ld(z) * ld(z)); }   // no overflow for i64
    Point3<ld> unit() const { return Point3<ld>(*this) / length(); }
    bool isZero() const { return x == 0 && y == 0 && z == 0; }   // exact for integer T
 
    friend istream& operator>>(istream& in, Point3& p) { return in >> p.x >> p.y >> p.z; }
    friend ostream& operator<<(ostream& out, const Point3& p) { return out << p.x << ' ' << p.y << ' ' << p.z; }
};
 
template <class T> T dot(const Point3<T>& a, const Point3<T>& b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
 
template <class T>
Point3<T> cross(const Point3<T>& a, const Point3<T>& b) {
    return Point3<T>(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
}
template <class T>   // (a - o) x (b - o): normal of triangle oab, length = twice its area
Point3<T> cross(const Point3<T>& o, const Point3<T>& a, const Point3<T>& b) { return cross(a - o, b - o); }

The dot product is unchanged from 2D. The cross product becomes a vector: perpendicular to both inputs, with length equal to the parallelogram area, direction by the right-hand rule.

The scalar triple product

$[a,b,c]=a\cdot (b\times c)=\left|\begin{array}{ccc}{a}_x& a_y& a_z\\ b_x& b_y& b_z\\ c_x& c_y& c_z\end{array}\right|$

This is the 3D orientation test:

Value	Meaning
$>0$	$d$ is above the plane $abc$ (right-hand rule)
$=0$	the four points are coplanar
$<0$	below

And $|[a,b,c]|$ is 6x the tetrahedron volume, or the volume of the parallelepiped spanned by $a,b,c$.

// six times the signed volume of tetrahedron abcd: > 0 iff d is on the side cross(a, b, c) points to.
template <class T>
T volume6(const Point3<T>& a, const Point3<T>& b, const Point3<T>& c, const Point3<T>& d) {
    return dot(cross(a, b, c), d - a);
}
// +1: d above plane abc (abc counter-clockwise seen from d), -1: below, 0: coplanar.
template <class T>
int orient(const Point3<T>& a, const Point3<T>& b, const Point3<T>& c, const Point3<T>& d) {
    T v = volume6(a, b, c, d);
    return (v > T(0)) - (v < T(0));
}

Proof: the triple product is signed volume, and orient's sign

$a\cdot (b\times c)$: $b\times c$ has magnitude equal to the area of the parallelogram spanned by $b,c$ and direction normal to that parallelogram (right-hand rule). Dotting with $a$ gives $|a|\cos \theta \cdot |b\times c|$ where $\theta$ is the angle between $a$ and the normal — exactly (base area) $\times$ (height of $a$ above the $b,c$-plane, signed by which side), i.e. the signed volume of the parallelepiped spanned by $a,b,c$. A tetrahedron is $\frac{1}{6}$ of the parallelepiped spanned by its three edge vectors from one vertex (the parallelepiped decomposes into 6 congruent-volume tetrahedra by symmetry — a standard solid-geometry fact), giving the volume formula; volume6 computes $[b-a,c-a,d-a]$, six times the tetrahedron $abcd$‘s volume, positive exactly when $d$ is on the side that $(b-a)\times (c-a)$ points to (the same sign convention as the parallelepiped’s orientation). $■$

Overflow warning: this is a degree-3 expression. With coordinates up to $10^6$ it reaches $10^{19}$ — already past long long. Use __int128 unless the coordinates are small.

Planes

A plane is $\{p:n\cdot p=d\}$ with normal $n$. From three points: $n=(b-a)\times (c-a)$, $d=n\cdot a$.

// the points p with dot(n, p) == d; n is a normal vector (not normalized).
template <class T>
struct Plane {
    Point3<T> n;
    T d;
    Plane(const Point3<T>& n_, T d_) : n(n_), d(d_) {}
    Plane(const Point3<T>& a, const Point3<T>& b, const Point3<T>& c) : n(cross(a, b, c)), d(dot(n, a)) {}  // not collinear
 
    T side(const Point3<T>& p) const { return dot(n, p) - d; }           // > 0 on the side n points to, 0 on the plane
    ld dist(const Point3<T>& p) const { return side(p) / n.length(); }    // signed
    Point3<ld> project(const Point3<T>& p) const {
        return Point3<ld>(p) - Point3<ld>(n) * (ld(side(p)) / (n.length() * n.length()));
    }
};
 
// line through p and q (p != q) against a plane: {1, intersection point}, {0, p} if parallel,
// {-1, p} if the line lies in the plane.
template <class T>
pair<int, Point3<ld>> linePlaneIntersection(const Plane<T>& pl, const Point3<T>& p, const Point3<T>& q) {
    T sp = pl.side(p), sq = pl.side(q);
    if (sp == sq) return {sp == T(0) ? -1 : 0, Point3<ld>(p)};
    return {1, Point3<ld>(p) + Point3<ld>(q - p) * (ld(sp) / ld(sp - sq))};
}

Task	Formula
Signed distance from $p$	$\frac{n\cdot p-d}{|n|}$
Which side is $p$ on?	sign of $n\cdot p-d$ (exact in integers)
Projection of $p$ onto the plane	$p-\frac{n\cdot p-d}{|n{|}^2}n$
Reflection across the plane	$p-2\frac{n\cdot p-d}{|n{|}^2}n$
Plane-plane intersection	a line with direction $n_1\times n_2$
Angle between planes	the angle between the normals

Proof: Plane's formulas

$n\cdot p-d=n\cdot p-n\cdot a=n\cdot (p-a)$, the component of $p-a$ along the unit normal $\hat{n}$ scaled by $|n|$ — exactly $|n|$ times the signed distance from $p$ to the plane (by definition of the plane as the level set $n\cdot p=d$, moving a unit step along $\hat{n}$ changes $n\cdot p$ by $|n|$). Dividing by $|n|$ gives the signed distance directly. Projection subtracts exactly that signed-distance multiple of $\hat{n}=n/|n|$ from $p$, landing on the plane (verified by substituting into $n\cdot (\cdot )=d$); reflection subtracts twice that, landing on the mirror point. linePlaneIntersection: parametrize the line $p+t(q-p)$; it hits the plane when $\mathrm{side}(p)+t(\mathrm{side}(q)-\mathrm{side}(p))=0$ (since side is affine in the point), giving $t=\mathrm{side}(p)/(\mathrm{side}(p)-\mathrm{side}(q))$ — matching the code, with the degenerate case ($\mathrm{side}(p)=\mathrm{side}(q)$) meaning the line is parallel to the plane, either disjoint from it (side(p) != 0) or lying entirely within it. $■$

Distances

// distance from p to the line through a and b (a != b).
template <class T> ld lineDist(const Point3<T>& p, const Point3<T>& a, const Point3<T>& b) {
    return cross(Point3<ld>(b - a), Point3<ld>(p - a)).length() / (b - a).length();
}
// distance from p to the segment ab (a == b allowed).
template <class T> ld segmentDist(const Point3<T>& p, const Point3<T>& a, const Point3<T>& b) {
    Point3<ld> d(b - a), ap(p - a);
    ld len2 = dot(d, d);
    ld t = len2 == 0 ? 0 : clamp(dot(ap, d) / len2, ld(0), ld(1));
    return (ap - d * t).length();
}
// distance between the segments ab and cd (degenerate segments allowed).
template <class T>
ld segmentDist(const Point3<T>& a, const Point3<T>& b, const Point3<T>& c, const Point3<T>& d) {
    ld res = min({segmentDist(a, c, d), segmentDist(b, c, d), segmentDist(c, a, b), segmentDist(d, a, b)});
    if (cross(b - a, d - c).isZero()) return res;    // parallel: the minimum is attained at an endpoint
    Point3<ld> u(b - a), v(d - c), w(a - c), n = cross(u, v);
    ld uu = dot(u, u), uv = dot(u, v), vv = dot(v, v), uw = dot(u, w), vw = dot(v, w), den = dot(n, n);
    ld s = (uv * vw - vv * uw) / den, t = (uu * vw - uv * uw) / den;   // closest points a + u s and c + v t
    if (0 <= s && s <= 1 && 0 <= t && t <= 1) res = min(res, fabsl(dot(w, n)) / n.length());
    return res;
}

Pair	Formula
Point to line	$\frac{|(p-a)\times d|}{|d|}$
Point to plane	$\frac{|n\cdot p-d|}{|n|}$
Line to line (skew)	$\frac{|(a_2-a_1)\cdot (d_1\times d_2)|}{|d_1\times d_2|}$
Point to segment / triangle	project and clamp
Point to sphere	$||p-c|-r|$

Proof: the distance formulas, and segment-segment via the skew-line reduction

Point to line. $|(p-a)\times d|$ is (parallelogram area spanned by $p-a$ and $d$) $=|d|\cdot (\text{height of }p\text{ above the line})$; dividing by $|d|$ leaves the height, i.e. the perpendicular distance — the direct 3D analogue of the 2D cross-product distance formula.
Segment to point. Projecting and clamping $t\in [0,1]$ finds the closest point on the segment exactly as in 2D (proved in Primitives); the 3D version is identical since the projection formula only ever uses dot products, unchanged in any dimension.
Skew line-to-line. $d_1\times d_2$ is a common perpendicular direction to both lines (perpendicular to each by the cross product’s definition). Projecting the connecting vector $a_2-a_1$ onto this common perpendicular gives the (constant) separation between the two lines along their shared normal — this is the unique distance that doesn’t depend on where along either line you measure from, since moving along either line changes the connecting vector only by a multiple of $d_1$ or $d_2$, both of which are perpendicular to $d_1\times d_2$ and so don’t affect that projection.
Segment to segment. The true minimum distance between two segments is either attained at an endpoint-to-opposite-segment distance (covered by the four segmentDist(point, segment) calls) or, when the segments are not parallel, at the unique pair of points realizing the skew-line distance — but only if that pair’s parameters $s,t$ both land within $[0,1]$ (otherwise the true closest points on the finite segments are clamped to an endpoint, which the first four terms already cover). Taking the minimum of all these candidates is exactly right: if the unconstrained skew-line closest points are both within their segments, that is the true minimum (no endpoint pair can do better, since it’s a global minimum over the infinite lines); otherwise the constrained optimum on the finite segments occurs at a boundary of the feasible $(s,t)$ box, which is exactly captured by one of the four endpoint-to-segment distances. Parallel segments (cross(b-a,d-c) is the zero vector) have no unique skew-line distance, so the minimum is necessarily attained at an endpoint, which the four-term min already covers without needing the skew computation at all. $■$

Tested: lineDist and segmentDist (point-to-segment) against brute-force fine sampling along the line/segment; segmentDist (segment-to-segment) against a fine 2D grid sample over $(s,t)\in [0,1{]}^2$ taking the minimum distance, on random segment pairs including parallel and skew configurations.

Volumes

Solid	Formula
Tetrahedron $abcd$	$\frac{1}{6}|[b-a,c-a,d-a]|$
Parallelepiped	$|[a,b,c]|$
Convex polyhedron	$\frac{1}{3}{\sum }_{\text{faces}}(\text{face centroid}\cdot \hat{n})A_{\text{face}}$ (divergence theorem)
Polyhedron from triangles	$\frac{1}{6}\sum volume6(a_0,a_i,b_i,c_i)$ over outward-oriented triangles, fixed apex $a_0$
Sphere	$\frac{4}{3}\pi r^3$
Intersection of two spheres	a lens; two spherical caps

The triangle-sum volume formula is the 3D shoelace: sum the signed tetrahedron volumes from a fixed apex to each face, and everything outside cancels — exactly the 2D polygon-area proof, one dimension up.

3D convex hull

$O(n\log n)$ expected with randomised incremental construction plus a conflict graph; $O(n^2)$ with the simple incremental algorithm below. Output is $O(n)$ faces (Euler: $V-E+F=2$).

Substantially harder than 2D — the degeneracies (coplanar points, collinear edges) are much nastier, and the orientation predicate overflows easily (see above).

Payoff: the 3D hull of points lifted onto the paraboloid $(x,y)\to (x,y,x^2+y^2)$ projects to the 2D Delaunay triangulation.

Rotations

Representation	Note
Rotation matrix ($3\times 3$ orthogonal)	composes by multiplication; 9 numbers
Euler angles	intuitive but suffers gimbal lock
Axis-angle (Rodrigues)	$p^{\prime }=p\cos \theta +(k\times p)\sin \theta +k(k\cdot p)(1-\cos \theta )$
Quaternions	4 numbers, no gimbal lock, composes by quaternion multiplication, interpolates smoothly

Quaternions are the standard in graphics and robotics. In competitive programming, rotations by fixed axis-aligned amounts ($90°$ about an axis) are much more common, and those are simple coordinate permutations with sign flips — there are 24 of them for a cube.

Proof: Rodrigues' rotation formula

Decompose $p$ into components parallel and perpendicular to the rotation axis $k$ (unit vector): $p=p_{\parallel }+p_{\perp }$ with $p_{\parallel }=k(k\cdot p)$ and $p_{\perp }=p-p_{\parallel }$. Rotation about $k$ fixes $p_{\parallel }$ and rotates $p_{\perp }$ within the plane perpendicular to $k$ by angle $\theta$: that plane has an orthonormal basis $\{p_{\perp },k\times p_{\perp }\}$ (since $k\times p_{\perp }$ is perpendicular to both $k$ and $p_{\perp }$, with the same length as $p_{\perp }$ because $k\perp p_{\perp }$), so the rotated perpendicular component is $p_{\perp }\cos \theta +(k\times p_{\perp })\sin \theta$. Noting $k\times p_{\perp }=k\times p$ (since $k\times p_{\parallel }=k\times k(k\cdot p)=0$), the full rotated vector is $p_{\parallel }+p_{\perp }\cos \theta +(k\times p)\sin \theta =p\cos \theta +(k\times p)\sin \theta +k(k\cdot p)(1-\cos \theta )$ (substituting $p_{\perp }=p-p_{\parallel }$ and collecting terms), matching the formula. $■$

Spherical geometry

Task	Formula
Great-circle distance	$r\arccos (\hat{u}\cdot \hat{v})$, or the numerically better haversine formula
Spherical triangle area	$r^2(A+B+C-\pi )$ — the spherical excess
Convert lat/long to Cartesian	$(\cos \varphi \cos \lambda ,\cos \varphi \sin \lambda ,\sin \varphi )$

Use the haversine form for great-circle distance: acos of a dot product loses precision badly for nearby points (the same ill-conditioning as acos near $\pm 1$, proved for the 2D angle formula in Primitives).

Proof: the spherical excess formula

This is Girard’s theorem: on a unit sphere, a spherical triangle’s three great-circle arcs, each extended to a full great circle, divide the sphere into 8 regions; summing the areas of the triangle and its “antipodal-reflection” partners via the angles $A,B,C$ (each great circle contributes a lune of area $2\times \text{angle}$ at the two points where two of the circles cross) and using that the total sphere area is $4\pi$, algebraic bookkeeping (each lune of angle $2A$ covers the triangle plus one adjacent region, and summing all three lunes covers the triangle three times plus each adjacent region once, while the lunes’ total area is $2(A+B+C)\cdot 2=4(A+B+C)$… the clean statement) gives triangle area $=(A+B+C-\pi )r^2$ on a sphere of radius $r$ (unit sphere scaled by $r^2$ for area). The full combinatorial derivation is Girard’s classical lune-counting argument; the result is standard and is what’s tested below.

Tested numerically: the spherical excess formula checked against subdividing a spherical triangle into many small flat triangles and summing (converging to the true spherical area) for several triangles on the unit sphere.

3D Hull

Convex hull of points in 3D as outward-oriented triangles (incremental, exact), with surface area and volume.

vector<Point3<i64>> p(n);
vector<array<int, 3>> faces = convexHull3D(p);   // {a, b, c}: cross(p[a], p[b], p[c]) points outwards
ld area = hullArea(p, faces);
i64 vol6 = hullVolume6(p, faces);                 // six times the volume, exact

Complexity: $O(nF)$ time, $F=$ number of hull faces ($O(n^2)$ worst case), $O(n)$ memory.

Handles duplicate, collinear and coplanar points exactly: a point is added only if it lies strictly outside some face, so points inside the hull or on its surface are skipped. Consequently a flat facet with $k>3$ corners is returned as coplanar triangles, but no face has zero area and the faces form a closed surface (merge coplanar neighbours if the problem counts facets). All-coplanar input returns {} (handle flat input with the 2D hull instead). Exact with T=i64 while $|\text{coord}|\le 10^5$ (orient multiplies three coordinates); beyond that, use Point3<i128>. $F\le 2n-4$.

// triangular faces {a, b, c} (indices into p) of the convex hull, CCW seen from outside, so
// cross(p[a], p[b], p[c]) points outwards. returns {} if all points are coplanar.
template <class T>
vector<array<int, 3>> convexHull3D(const vector<Point3<T>>& p) {
    int n = p.size(), b = -1, c = -1, d = -1;
    for (int i = 1; i < n && b < 0; i++) if (p[i] != p[0]) b = i;
    for (int i = 1; i < n && b >= 0 && c < 0; i++) if (!cross(p[0], p[b], p[i]).isZero()) c = i;
    for (int i = 1; i < n && c >= 0 && d < 0; i++) if (orient(p[0], p[b], p[c], p[i]) != 0) d = i;
    if (d < 0) return {};
    if (orient(p[0], p[b], p[c], p[d]) > 0) swap(b, c);
    vector<array<int, 3>> faces = {{0, b, c}, {0, d, b}, {b, d, c}, {c, d, 0}}, kept;
    vector<i64> edges;   // directed edges u -> v of the faces visible from the new point, as u * n + v
    for (int i = 0; i < n; i++) {
        kept.clear();
        edges.clear();
        for (auto& f : faces) {
            if (orient(p[f[0]], p[f[1]], p[f[2]], p[i]) > 0) {
                for (int k = 0; k < 3; k++) edges.push_back(i64(f[k]) * n + f[(k + 1) % 3]);
            } else {
                kept.push_back(f);
            }
        }
        if (edges.empty()) continue;   // p[i] is inside the hull or on its surface
        sort(edges.begin(), edges.end());
        for (i64 e : edges) {
            int u = e / n, v = e % n;
            if (!binary_search(edges.begin(), edges.end(), i64(v) * n + u)) kept.push_back({u, v, i});   // horizon edge
        }
        swap(faces, kept);
    }
    return faces;
}
 
template <class T> ld hullArea(const vector<Point3<T>>& p, const vector<array<int, 3>>& faces) {
    ld res = 0;
    for (auto& f : faces) res += cross(p[f[0]], p[f[1]], p[f[2]]).length() / 2;
    return res;
}
// six times the volume enclosed by outward-oriented faces (exact for integer T).
template <class T> T hullVolume6(const vector<Point3<T>>& p, const vector<array<int, 3>>& faces) {
    T res = 0;
    for (auto& f : faces) res += volume6(p[faces[0][0]], p[f[0]], p[f[1]], p[f[2]]);
    return res;
}

Proof: convexHull3D is correct

Initialization. $p[0],p[b]$ are the first two distinct points; $p[c]$ is the first point not collinear with them (making $\mathrm{cross}(p_0,p_b,p_c)$ nonzero, a valid triangle normal); $p[d]$ is the first point not coplanar with that triangle (orient != 0). The four points $0,b,c,d$ then form a valid tetrahedron, oriented so all four initial faces point outward (the swap(b,c) corrects the orientation if the initial triangle $0,b,c$ was wound so that $d$ ended up on its positive side, which would make the face list inconsistently oriented).
The incremental step (invariant). Before processing point $i$, faces is the correct hull of $p[\mathrm{0..}i)$ (points already processed), each face’s vertex order outward-oriented. For the new point $p[i]$: a face is “visible” from $p[i]$ exactly when $p[i]$ is strictly on its outer side (orient(...) > 0) — such a face can no longer be on the hull once $p[i]$ is included, since $p[i]$ now lies outside it. Discarding all visible faces and connecting $p[i]$ to the horizon — the boundary between visible and non-visible faces — restores a valid convex hull containing $p[i]$. The horizon is found by the standard “directed edge cancellation” trick: every edge of a visible face is recorded in one direction; an edge shared between two visible faces appears in both directions (once per face, oppositely oriented, since adjacent faces traverse a shared edge in opposite directions by the outward-orientation convention) and cancels; only horizon edges (shared between one visible and one non-visible face, or on the silhouette) survive uncancelled, and each becomes a new triangle with apex $p[i]$.
Degenerate points are skipped correctly. If no face is visible (edges.empty()), $p[i]$ is inside the hull or exactly on its surface — adding it would change nothing, so it’s correctly skipped, preserving exactness for duplicate and surface points.
hullArea/hullVolume6. Each face’s triangle area is half the cross product’s length (the standard area formula); summing gives total surface area. For the volume, summing signed tetrahedra from a fixed apex $p[faces[0][0]]$ to every face is the 3D shoelace: for a convex (or any closed, consistently outward-oriented) polyhedron, this telescopes to exactly $6\times$ the enclosed volume, by the identical origin-independence argument as the 2D polygon-area proof — regions “behind” the apex relative to some faces contribute negative volume that exactly cancels the double-counted part, leaving the true volume regardless of which vertex is chosen as apex. $■$

Tested: convexHull3D checked on a cube, regular tetrahedron, and octahedron (correct vertex/face counts via Euler’s formula $V-E+F=2$, correct volume and surface area against the known closed forms); against a brute-force 3D hull check (a candidate triangle is a hull face iff all other points are weakly on one side of its plane) on small random point sets; and on inputs with duplicate points, coplanar subsets, and all-coplanar input (correctly returns {}).