Rethinking geometry libraries

2015-09-10 Permalink

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 an ApproximatelyEqual 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:

  const double EPS_DOUBLE = 1e-10; // !??
  bool operator ==(const Vec &u, const Vec &v); // Compares up to EPS_DOUBLE
  bool ExactlyEqual(const Vec &u, const Vec &v);

• Defining operator < for your vectors just because you need a map with them for keys is fine. It’s a bit strange though, as it makes no geometrical sense.[1] You can use a custom comparator instead. A hash-map 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 non-member 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 -.

Boxes (axis aligned)

• Axis-aligned boxes are mostly used for bounding-box 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 bounding-box 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 floating-point 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 y-axis. 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 set-theoretic union represented by multiple disjoint boxes. These are useful in some circumstances.

Offender: Qt.

Order-theoretic way

Rather than starting from their geometric meaning as special subsets of ℝⁿ, start with the order-theoretic notion of the lattice (ℝⁿ×ℝⁿ, contains⁻¹) 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₀, y₀), min(x₁, y₁))
• join(x, y) = (min(x₀, y₀), max(x₁, y₁)).[3]

The geometric interpretation of the box x is given by

• S(x) = { v ∈ ℝⁿ | ∀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

class Matrix3d {
public:
	double mat[3][3];
	/* ... */
};
class OrthogonalMatrix : public Matrix3d { /* ... */ };
class RotationMatrix : public OrthogonalMatrix { /* ... */ };

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 self-documentary purpose, then a simple old typedef would do.

By the way, the use of OrthogonalMatrix was virtually null.

Footnotes

1. For dimensions ≥ 2.
2. Liberate operators movement.
3. Here min and max are elementwise.
4. The opposite is true if x and y aren’t empty.
5. Or a so(3) vector for infinitesimal rotations.