原创 Random thoughts (Part 2)

[Continued from Random thoughts (Part 1)]

A 3×3 matrix class?

It's this character flaw that keeps me from being able to show you a completed 3×3 matrix class. Until recently, I've been stuck in one of those while-loops, and this time the ex-wife isn't here to assert the interrupt. I wasn't able to exit the loop until I began writing this column. I'd like to tell you what I have, and see what you think.

Clearly, a 3×3 matrix is not just any old matrix, even any old matrix with one dimension of three. It's a square matrix, and square matrices have special properties. For one thing, they have well-defined diagonals. A square matrix might, in fact, be diagonal. That is, its off-diagonal elements may all be zero. That's important, because a diagonal matrix has a trivial inverse:

A square matrix might also be symmetric. If so, it can easily be transformed to a diagonal matrix. When I'm doing arithmetic with 3×3 matrices, they're almost always associated with coordinate rotations. For example, to convert a vector from a rotating coordinate system to a fixed one, I might write:

The matrix T is even more specialized than most. It's a rotation matrix, also called an orthonormal matrix. It has two very important attributes: Its determinant is 1, and its inverse is the same as its transpose:

This is very important, because a transpose of a 3×3 matrix is trivial to generate, whereas an inverse takes a lot longer. Also, because the determinant is known, we don't have to worry about the matrix being singular.

Finally, there exist more than one way to represent a rotation. The matrix is one choice, but not the only one, or even the best one. Other choices include Euler angles and quaternions. Euler angles are useful because people tend to be able to visualize the rotations better. But they're terrible choices for computations. Quaternions are the best and most compact, but good luck trying to visualize them.

Side Comment: You know you've arrived as a guru when you realize you can visualize a quaternion. Since it's a four-dimensional vector, that takes a little practice.

Decisions, decisions

Now that I've given you the background, I can tell you why I was stuck in the while-loop. I know that I want to use 3×3 matrices mostly to effect coordinate rotations. But is that always going to be the case? Maybe not. And if it is, may I include constructors to convert a set of 3-vectors to a matrix? Or a set of Euler angles, or quaternions? If I demand that the matrix must be orthonormal, I'd better not try to create one from vectors.

For that matter, if I'm dealing exclusively with rotation matrices, maybe I need a function to make sure they stay orthonormal.

Here's why. When we're writing dynamic simulations involving rotations, we have to numerically integrate the elements that make up the state vector. If one of those elements is a rotation matrix (or equivalent structure), numerical roundoff errors can cause the matrix to drift away from its nominally orthonormal state. Those of us who write such simulations often include functions to re-normalize the matrix. It's not an easy thing to do, and doing it in an optimal way is even harder.

So a re-normalizing function is going to be almost a necessity for rotation matrices, but it's going to look really, really strange in a general-purpose matrix package—even one specialized to 3×3 matrices.

Of course, you know the conceptual solution to problems such as this: inheritance. Write a C++ base class for 3×3 matrices. Then define a derived class that specializes the class even further, to rotation matrices. Most of the operations will be the same, but a few functions, such as re-normalization, can be added, and the inverse function can be changed to invoke the simpler transpose function.

There's only one problem with this approach. A rotation in 3-space is a very unique thing, and the function that it performs is vastly different from the more generic matrix product. Heck, the member data for a "rotation matrix" need not even include a matrix. Rotations can also be described by a quaternion or a set of Euler angles.

And because there are three possible ways to describe the rotation, we're also going to need conversion functions to transform one set of descriptors to another.

FYI, a couple of years ago I wrote a library of Matlab functions capable of converting between angle, matrix, and quaternion representations. In doing so, I even surprised myself, coming up with killer algorithms that resulted in code that was tight, accurate, and bulletproof. I'll be sharing the algorithms with you soon.

Most people who work every day with rotation-related problems choose the quaternion, because it's more efficient to do so (no angles involved, and fewer elements to store). One has to ask: what's the point having a rotation matrix be a derived class of the 3×3 matrix, if it doesn't even hold a matrix anywhere in its structure?

Off and on, I've wrestled with questions like this many times. I can't say that I've ever come up with the definitive solution. More often, I find myself in need of that priority interrupt.

But this time, I've finally got the structure I need. The trick is, there should be no class called RotationMatrix. Instead, there's a class called Rotation. When applied to a vector or another matrix (or Rotation), the operation may look like a matrix product, but the appearance is only superficial. Internally, a Rotation should be a separate type of object, worthy of its own class. Heck, it may not—and probably won't—even have a matrix as part of its member data. And in the best tradition of object-oriented programming, the way the rotation class is implemented should be transparent to the user. If I decide to change the internal representation from matrix to quaternion, that decision shouldn't matter one iota to users of the class.

When you think about it, the fact that I can define a rotation operation that looks like a matrix product is no different than defining products of integers or real numbers. I can declare:

int i, j, k;

double x, y, z;

Vector X, Y;

Matrix A;

Rotation R;

and write:

k = i*j;

z = x*y;

Y = A*X;

Y = R*X;

Yes, the statements look the same and mean much the same. That's the elegance of operator overloading, and it makes life incredibly easy for the programmer. But we have no particular reason for a rotation object to inherit the '*' operator from a matrix, any more than we might have it do so from class int. Separating the class Rotation into a full-blown separate class is beginning to make a lot of sense.

Wrapping up

I'm glad we had this little talk. So often in life, I've found that as I try to explain a knotty problem to someone else, I realize that I suddenly know the solution. That seems to be what has happened here. Thanks for listening.