Rethinking geometry libraries
20150812
You may think that the following are obviously true, or you may find them controversial. Yet the fact that I’ve seen these issues in real code and that they caused much pain makes this, nevertheless, important.
So here we go.
Vectors

Points—Affine spaces are a kind of a mathematical abstraction that has no right to exist in code. In real life we first choose our coordinate system—origin and basis—then encode the points with vectors within the formed vector space. There is no other way around. Thus separating points from vectors in a computer is more of a burden. Let’s forget what they taught in the linear algebra class about the distinction of linear and affine spaces, and go with a single vector type.
Some libraries provide a point type without a complete set of vector operations, and no vector type. This leads to much verbosity when doing even simple math with those types. Instead, please, provide a full set of vector operations. Then call them vectors, because seeing a velocity vector or a rotation axis of type
Point
raises an eyebrow. 
operator ==
shall do exact comparison (total order or not). Comparing up to epsilon voids transitivity (which does break code). Furthermore, there isn’t any meaningful choice of epsilon because it depends on the algorithm used at the caller. There is no need for anApproximatelyEqual
function because it can be stated more concisely by norm(a − b) ≤ epsilon, making both the norm used and the threshold explicit.Funny that people at VM got it exactly backwards:
Defining
operator <
for your vectors just because you need amap
with them for keys is fine. It’s a bit strange though, as it has no geometrical sense.[1] You can use a custom comparator instead. A hashmap will serve you even better. They shall be PODs and be layout compatible with a simple array of size n. You have no right to necessitate me calling your factory function to create a dynamically allocated concrete point instance, derived from some abstract class, then setting each of its coordinates by a separate virtual setter call—all this just to convert between two coordinate systems (offender: BMG). Also there is no excuse to swap x and y (offender: KDU).
Binary operators shall be free functions.[2] This way they are symmetric wrt. overload resolution, and can be moved to a separate header. See How nonmember functions improve encapsulation.
Elementwise operations are surprisingly frequent. Rather than defining
operator *
as dot product (or cross product?), let*
and/
be elementwise in consistence with+
and
.
You might think that I like GLSL—this is true. However, I’d followed these conventions even before GLSL existed. And this is not a coincidence. GLSL was designed for writing practical highperformance computer graphics code rather than to satisfy somebody’s hallucinations of an OO pure world.
Boxes (axis aligned)
 Axisaligned boxes are mostly used for boundingbox calculations and space partitions. For both applications a representation as two corners is the correct one, because you cannot reliably tile the plane with floating point boxes represented as position plus extent. The exact bounding box of two boxes may be not representable in the later case. It’s also a matter that from calculating the boundingbox till the rendering of the rectangle on the screen, say, no part of it needs the extents.
 Boxes shall be half open—People smarter than me had already given an explanation. For geometry libraries in particular, it guarantees consistency across integer and floatingpoint boxes. It makes remapping of plane partitions easier.
 For that reason don’t use min and max for the bounds. Low and high, or a and b are other names to consider.
 Don’t use the top and bottom terms. Their interpretation depends on the direction of the yaxis. If in two dimensions there are only two popular choices, generalizing for higher dimensions is even harder. In three dimensions there is a plethora of conventions. Is y → +∞ the top, bottom, front or back face? Use indices for axes instead.
 Don’t enforce low ≤ high. Using [+∞,+∞,−∞,−∞] boxes is useful as the starting value in bounding box calculations, facilitating branchless code.
 A function that returns the bounding box of two boxes is not a union. There can be a true union function, along with difference, which returns the settheoretic union represented by multiple disjoint boxes. These are useful in some circumstances.
Offender: Qt.
Ordertheoretic way
Rather than starting from their geometric meaning as special subsets of ℝ^{n}, start with the ordertheoretic notion of the lattice (ℝ^{n}×ℝ^{n}, contains^{−1}) where:
 contains(x, y) iff ∀i x_{0,i} ≤ y_{0,i} and y_{1,i} ≤ x_{1,i}.
The corresponding meet and join operators are
 meet(x, y) = (max(x_{0}, y_{0}), min(x_{1}, y_{1}))
 join(x, y) = (min(x_{0}, y_{0}), max(x_{1}, y_{1})).[3]
The geometric interpretation of the box x is given by
 S(x) = { v ∈ ℝ^{n}  ∀i x_{0,i} ≤ v_{i} < x_{1,i} }
 empty(x) iff ∀i x_{1,i} ≤ x_{0,i}.
Defer the geometric interpretation until necessary. It is implied that:
 S(x) ∩ S(y) = S(meet(x, y))
 S(x) ∪ S(y) ⊆ S(join(x, y)) is the bounding box of x and y.
 empty(x) iff S(x) = ∅.
 contains(x, y) implies S(y) ⊆ S(x).[4]
Rotation matrices
This code is real. Do you see the problem? This is much like the Circle And Ellipse Problem. All those functions working on Matrix3d
can easily invalidate the assumption that RotationMatrix
represents a rotation, even remotely. Separating RotationMatrix
from Matrix3d
won’t save it from being ruined by floating point operations anyway.
For rotations, either use a matrix with a code prepared to digest any matrix fed in, or a quaternion.[5] Don’t try to solve this by hiding the representation. Matrices and quaternions have different numerical properties, and such properties will leak through the interfaces, much like operations over std::vector
and std::list
have different guarantees. If you just want to use RotationMatrix
for selfdocumentary purpose, then a simple old typedef
would do.
By the way, the use of OrthogonalMatrix
was virtually null.
Footnotes
 For dimensions ≥ 2.
 Liberate operators movement.
 Here min and max are elementwise.
 The opposite is true if x and y aren’t empty.
 Or a so(3) vector for infinitesimal rotations.