Tag Archives: Math

Unit Testing Ambiguity

This article was originally published on September 25, 2010 at Pulsar Engine.

It has been updated to employ the conventions we use in Sauce instead of the ones I had previously used in Pulsar — aside from that, the content remains unchanged. This article still reflects my views and code.

Background

This post is related to the previous post on testing a return value that consists of a set of items with no predefined order (see Unit Testing Unordered Lists)[1]. However, it differs in that this time the problem is a little more mathematical in nature: how do we properly test a routine that performs an operation that may have multiple correct solutions?

For those of you who have written core math routines and solvers before, you know exactly what I’m referring to, but even if you haven’t I still encourage you to continue reading. I’ve chosen a relatively simple, common, yet important example to work through — there’s a good chance you’ll make use of this knowledge somewhere.

Unit Testing Ambiguity

There have been times when I’ve found myself trying to impose unrealistic expectations on my routines. Of course, when I test for them, I end up with failed tests and my first thought is more often than not: the code being tested must be wrong. However, the reality is that I have actually written a bad test.

A good example of this happened to me when I was working on my ComputeOrthoVectors() routine[2]. The purpose of the function is: given a single Vector3, compute two additional Vector3s which form an orthonormal basis with the input Vector3. So, I wrote the following test:

TEST(CanComputeOrthoBasisVectors)
{
   Vector3 vecA;
   Vector3 vecB;
   ComputeOrthoBasisVectors(vecA, vecB, Vector3::UNIT_X);

   CHECK_ARRAY_CLOSE(Vector3::UNIT_Y, vecA, 3, EPSILON);
   CHECK_ARRAY_CLOSE(Vector3::UNIT_Z, vecB, 3, EPSILON);
}

The primary issue here is that the operation itself is soaked in ambiguity. More specifically, the two Vector3s that make up the return value can be produced by one of many valid solutions. To see this, let’s look at the math required for the implementation.

Computing an Orthonormal Basis

First, let’s be clear on the terminology: a set of three Vector3s that are each perpendicular to the other two is called an orthogonal basis. If we normalize the vectors in the set, we can call it an orthonormal basis. In other words: A 3D orthonormal basis is a set of three vectors that are each perpendicular to the other two and each of unit length.

The need for creating an orthonormal basis from a single vector is a fairly common operation (probably the most common usage is in constructing a camera matrix given a “look-at” vector).

Given two non-coincident[3] vectors, we can use the cross product to find a third vector that is perpendicular to both the input vectors. This is often why we say that two non-coincident vectors span a plane, because it is from these two vectors that we can compute the normal to that plane.

The thing to note here is that there are an infinite number of valid non-coincident vectors that span a plane. You can imagine grabbing the normal vector as if it were a rod connected to two other rods (the input vectors) and spinning it; any orientation you spin it to is a valid configuration that would produce the same normal vector. I have created an animation demonstrating this below:

Animation of two vectors orthogonal to input vector. The two vectors remain in the plane defined by the input vector (the plane normal).

In essence, this is why the results of the operation are valid, yet, for the lack of a better word, ambiguous. In the case of our routine, we are supplying the normal and computing two other vectors that span the plane. The fact that the two returned vectors are orthogonal to the input vector does not change the fact that there are an infinite number of configurations.

Testing the Geometry, Not the Values

Returning to the original testing dilemma, since there are an infinite number of possible solutions for our returned pair of Vector3s, if we write tests that check the values of those vectors, then our tests will be bound to the implementation. The result is a ticking time-bomb that may explode in our face later on (maybe during an optimization pass) because, although we may change the pair of vectors returned for a given input, the three vectors still form an orthonormal basis — the correct and desired result — in spite of the fact that we now have failing tests.

My solution was to check the following geometric properties:

  1. All three vectors are perpendicular to each other (in other words, they form an orthogonal basis).
  2. The two returned vectors are of unit length (input vector is assumed to be normalized).

The following are a few tests (straight from my codebase) that employ this approach:

TEST(CanComputeOrthoBasisVectors_UnitX)
{
   Vector3 vecA;
   Vector3 vecB;
   ComputeOrthoBasisVectors(vecA, vecB, Vector3::UNIT_X);

   CHECK_CLOSE(0.0f, vecA.Dot(vecB), EPSILON);
   CHECK_CLOSE(0.0f, vecA.Dot(Vector3::UNIT_X), EPSILON);
   CHECK_CLOSE(0.0f, vecB.Dot(Vector3::UNIT_X), EPSILON);
   CHECK_CLOSE(1.0f, vecA.Mag(), EPSILON);
   CHECK_CLOSE(1.0f, vecB.Mag(), EPSILON);
}


TEST(CanComputeOrthoBasisVectors_UnitY)
{
   Vector3 vecA;
   Vector3 vecB;
   ComputeOrthoBasisVectors(vecA, vecB, Vector3::UNIT_Y);

   CHECK_CLOSE(0.0f, vecA.Dot(vecB), EPSILON);
   CHECK_CLOSE(0.0f, vecA.Dot(Vector3::UNIT_Y), EPSILON);
   CHECK_CLOSE(0.0f, vecB.Dot(Vector3::UNIT_Y), EPSILON);
   CHECK_CLOSE(1.0f, vecA.Mag(), EPSILON);
   CHECK_CLOSE(1.0f, vecB.Mag(), EPSILON);
}


TEST(CanComputeOrthoBasisVectors_UnitZ)
{
   Vector3 vecA;
   Vector3 vecB;
   ComputeOrthoBasisVectors(vecA, vecB, Vector3::UNIT_Z);

   CHECK_CLOSE(0.0f, vecA.Dot(vecB), EPSILON);
   CHECK_CLOSE(0.0f, vecA.Dot(Vector3::UNIT_Z), EPSILON);
   CHECK_CLOSE(0.0f, vecB.Dot(Vector3::UNIT_Z), EPSILON);
   CHECK_CLOSE(1.0f, vecA.Mag(), EPSILON);
   CHECK_CLOSE(1.0f, vecB.Mag(), EPSILON);
}


TEST(CanComputeOrthoBasisVectors_RefPaperExample)
{
   Vector3 vecA;
   Vector3 vecB;
   const Vector3 unitAxis(-0.285714f, 0.857143f, 0.428571f);
   ComputeOrthoBasisVectors(vecA, vecB, unitAxis);

   CHECK_CLOSE(0.0f, vecA.Dot(vecB), EPSILON);
   CHECK_CLOSE(0.0f, vecA.Dot(unitAxis), EPSILON);
   CHECK_CLOSE(0.0f, vecB.Dot(unitAxis), EPSILON);
   CHECK_CLOSE(1.0f, vecA.Mag(), EPSILON);
   CHECK_CLOSE(1.0f, vecB.Mag(), EPSILON);
}

By testing these properties, as opposed to testing for the resulting vector values directly (as in the original test shown above), it doesn’t matter how the internals of the ComputeOrthoBasisVectors() produces the two returned Vector3s. As long as the input vector and the returned vectors all form an orthonormal basis, our tests will pass.

Final Thoughts

My hope is that the example presented in this article demonstrates one of the pitfalls of having tests that depend on the internal implementation of the routine being tested. As I have stated before, although it is important to write tests for the functionality, it can be difficult to recognize when a test is bound to the implementation.

A good place to start looking for this sort of scenario is in tests that explicitly check return values. Certainly, explicit checks is actually what you want in most cases, but for some operations it is not.

Footnotes

  1. Originally, these two topics were going to be addressed in a single article, but I decided against it in hopes that keeping them separate would allow for more clarity in each.
  2. In fact, the trouble I had in testing the ComputeOrthoVectors() routine is what inspired me to post both this article and the last.
  3. Two vectors are said to be coincident if they have the same direction when you discard their magnitudes. In other words, two vectors are coincident if you normalize both of them and the results are the same.

References

  • Möller, Tomas and John F. Hughes. Building an Orthonormal Basis from a Unit Vector. [ pdf ]

Covariance Matrix

This article was originally published on December 6, 2009 at Pulsar Engine.

It has been updated to employ the conventions we use in Sauce instead of the ones I had previously used in Pulsar — aside from that, the content remains unchanged. This article still reflects my views and code.

This weekend I was going through a section in Real-Time Collision Detection on computing the covariance matrix from a set of points (from hereon “point cloud”). Of course, while working through the implementation it’s always good to have an example or two to help solidify the concept, so I found myself working through a few samples for my unit tests and decided to post them here.

Example 1

Let’s start with a simple, somewhat uninteresting example where we have a single point on each of the cardinal axes.

+x:  1.0,  0.0,  0.0
+y:  0.0,  1.0,  0.0
+z:  0.0,  0.0,  1.0
-x: -1.0,  0.0,  0.0
-y:  0.0, -1.0,  0.0
-z:  0.0,  0.0, -1.0

Covariance Matrix:

0.333333,  0.0,  0.0
0.0,  0.333333,  0.0
0.0,  0.0,  0.333333

This is what we would expect because the spread is even along all the axes.

Example 2

Now, let’s look at a rotated version of the point cloud used in Example 1. Rotate the points 45 deg about the y-axis. I’m going to use 0.5 instead of sqrt(2)/2 to illustrate what happens when there is a dominant axis.

+x:  1.0,  0.0,  0.0   ->    0.5,  0.0,  0.5
+y:  0.0,  1.0,  0.0   ->    0.0,  1.0,  0.0  (stays the same)
+z:  0.0,  0.0,  1.0   ->   -0.5,  0.0,  0.5
-x: -1.0,  0.0,  0.0   ->   -0.5,  0.0, -0.5
-y:  0.0, -1.0,  0.0   ->    0.0, -1.0,  0.0  (stays the same)
-z:  0.0,  0.0, -1.0   ->    0.5,  0.0, -0.5

Covariance Matrix:

0.166667,  0.0,  0.0
0.0,  0.333333,  0.0
0.0,  0.0,  0.166667

This result should make sense because the x- and z-axis are no longer as spread as wide as the two points on the y-axis. As a result, the y component of the diagonal is larger than the x and z counterparts.

Also, it should be noted that even though the points in the xz-plane were not on the cardinal axes, the result is still a diagonal matrix. This is because the point cloud in this example is symmetric, such that each point cloud point in the xz-plane has a corresponding point (x,y) -> (-x,-y).

Example 3

In this example, we use the same point cloud as in Example 2, but translate the points all by 1.1, -0.4, 0.7.

 1.6, -0.4,  1.2
 1.1,  0.6,  0.7
 0.6, -0.4,  1.2
 0.6, -0.4,  0.2
 1.1, -1.4,  0.7
 1.6, -0.4,  0.2

Covariance Matrix:

0.166667,  0.0,  0.0
0.0,  0.333333,  0.0
0.0,  0.0,  0.166667

This example is a test to confirm a covariance matrix property: the covariance matrix remains the same if the point cloud is translated.

Example 4

Okay, so now let’s just take an arbitrary point cloud of eight points. The example here has no built-in symmetry, nor is it centered at the origin.

 1.2,  1.2,  1.2
-0.8, -0.8, -0.8
 0.7,  0.7,  0.5
 0.3,  0.4, -0.7
-0.2,  1.1,  0.5
 1.3, -0.8,  0.9
-0.1, -0.1, -0.3
 0.4, -0.5, -0.7

Centroid: 0.35, 0.15, 0.075

Covariance Matrix:

0.447500,  0.102500,  0.353750
0.102500,  0.582500,  0.283750
0.353750,  0.283750,  0.551875

Implementation Notes

The covariance matrix is symmetric. Therefore, only the upper triangular entries (including the diagonal) must be computed.

I included the centroid in the final example since it is subtracted off the points in the point cloud before computing the covariance matrix entries. We do this because we want our result to reflect the spread of the points in the local space of the point cloud.