Tuesday, January 3, 2017

Implementing 4D rotation

In my last post, I explored various options for allowing a player to control their orientation in 4D space. Now, I want to look at how 4D rotations can be implemented, regardless of the control scheme.

For my 4D maze game, the only thing that need to be rotated is the virtual camera that provides the player's first-person point of view. Dealing only with that eliminates a lot of accidental complexity that would be involved in applying arbitrary rotations to models in the world, but the procedures described here can be extended for that purpose.

The first thing to decide is how to represent the player's current orientation. This could be done with hyperspherical coordinates, which work pretty much just like regular spherical coordinates but with some extra angles added. However, working with hyperspherical coordinates gets really confusing really fast, and they suffer from the same limitations as spherical coordinates regarding gimbal lock. So, we'll use a different representation: a set of four unit vectors which define a basis for an egocentric coordinate system, in terms of the game map coordinate system.

While the game map coordinates are given in terms of x, y, z, and w axes, the egocentric coordinate system we use for orienting the camera translates these into forward, right, up, and ana (or F, R, U, A) components. It's important to make sure these basis vectors for the egocentric coordinate system always remain mutually perpendicular (i.e., that they form an orthonormal basis); that means there is some redundancy in this system. It would be possible to use a leaner representation, such as 3 basis vectors where the 4th can be recovered computationally as needed, or a pair of right-isoclinic and left-isoclinic quaternions, but despite the redundancy the use of explicit basis vectors makes the algorithms for arbitrary rotations, and the math employed therein, much easier to follow, and thus less prone to implementation error and easier to debug.

The basic 3D rotations- pitch, roll, and yaw- correspond to movement in the forward-up (FU), up-right (UR) and forward-right (FR) planes, respectively. In order to execute, for example, a pitch up, we would want to rotate the forward and up vectors in the FU plane so that the new forward vector points a little more in the old up direction, and new up vector points a little more in the old negative forward direction. Thus, rotation in a plane can be expressed in terms of rotating each of the basis vectors (or any vector that falls in he same plane) towards or away from another vector in that plane. Vectors in the perpendicular plane are unaffected by rotations, so the remaining bases of our egocentric coordinate system can just be left alone, and there aren't any other vectors to worry about!

Therefore, the basic building-block we'll need to implement arbitrary rotations of our egocentric coordinates is a function that can rotate some vector in the plane defined by itself and a second, perpendicular vector. It turns out that this is much simpler than the general formula for rotating an arbitrary point parallel to an arbitrary plane. The code looks like this:

// Rotate a vector v in the vk plane, // by angle t, assuming v & k are orthogonal function vec_rot(v, k, t){ let cos = Math.cos(t); let sin = Math.sin(t); return { x: v.x*cos + k.x*sin, y: v.y*cos + k.y*sin, z: v.z*cos + k.z*sin, w: v.w*cos + k.w*sin }; }

We gain quite a lot of simplifications when we can assume that all of our vectors are normalized and orthogonal!

A complete rotation of the camera requires two calls to this function- one for each of the basis vectors in the plane of rotation. If we did just one, they would no longer be orthogonal, and that would cause all kinds of problems! The code for executing a complete pitch rotation thus looks like this:

Player.prototype.pitch = function(t){ let fprime = vec_rot(this.f, this.u, t); let uprime = vec_rot(this.u, this.f, -t); this.f = fprime; this.u = uprime; };

Note that the first basis vector moves towards the second, while the second basis vector moves away from the first; thus, the second half-rotation uses a negative angle. Also note that this implementation re-calculate the same sine and cosine values (up to a sign difference) twice, so a slight efficiency improvement can be gained by inlining vec_rot and eliminating the extra calculations, at the cost of some code duplication and a little loss of clarity. The implementations for roll and yaw are the same, simply replacing f and u with u and l or f and l, respectively.

Having defined our 3D rotations in terms of planes rather than axes, the generalization to 4D rotations is trivial: just add three more rotation functions that act in the ana-axis planes: FA, UA, and LA. Arbitrary rotations in non-basis planes can then be achieved simply by calling these basic functions in series. It's not quite as efficient as building a 4D rotation matrix to do everything at once, but it's good enough most of the time, and very easy to understand and hack on.

There is, however, one problem with the system as it stands- due to the accumulation of floating-point errors, after a large number of rotations the basis vectors for our egocentric coordinates can become denormalized and non-orthogonal, which can cause a variety of strange rendering errors. Thus, it is prudent to periodically re-normalize our coordinate system. This consists of two parts: normalizing each basis vector, and ensuring that the basis vectors are all orthogonal to each other.

Normalizing a single vector is fairly easy:

function normalize(v){ let {x,y,z,w} = v; let len = Math.sqrt(x*x+y*y+z*z+w*w); return {x:x/len, y:y/len, z:z/len, w:w/len}; }

And it's also not that hard to take a vector and get the next closest vector to it that is perpendicular to some other vector:

// subtract the projection of v onto k from v function orthogonalize(v,k){ let {x:vx,y:vy,z:vz,w:vw} = v; let {x:kx,y:ky,z:kz,w:kw} = k; let kk = kx*kx+ky*ky+kz*kz+kw*kw; let vk = vx*kx+vy*ky+vz*kz+vw*kw; let scale = vk/kk; return { x: vx - kx*scale, y: vy - ky*scale, z: vz - kz*scale, w: vw - kw*scale }; }

If we know that k is pre-normalized, this can be further simplified:

// subtract the projection of v onto k from v, // where k is assumed to be normalized function orthogonalize(v,k){ let {x:vx,y:vy,z:vz,w:vw} = v; let {x:kx,y:ky,z:kz,w:kw} = k; let vk = vx*kx+vy*ky+vz*kz+vw*kw; return { x: vx - kx*vk, y: vy - ky*vk, z: vz - kz*vk, w: vw - kw*vk }; }

When combining these, it is important to remember to orthogonalize first, then normalize. Otherwise, orthogonalization will undo normalization. Orthonormalizing any two vectors is thus made very simple:
let fprime = normalize(this.f); let rprime = normalize(orthogonalize(r, fprime));

But what about the remaining plane? In 3D, we could just generate a new one by calculating the cross-product, but a binary cross-product doesn't exist in higher dimensions, for the simple reason that there is no longer a unique (up to a sign) direction perpendicular to any two vectors. In the general case, we would need an algorithm for computing the null space of a matrix to figure out what the perpendicular vectors should be... but once again, we can take advantage of the fact that the vectors we need to be perpendicular to are themselves already orthonormalized to do some major simplifications.

Fixing the last two axes just consists of repeatedly calling orthogonalize for each of the previously-orthonormalized axes, and then normalizing the result.

// Ensure the egocentric basis is orthonormal Player.prototype.renormalize = function(){ let (f, r, u, a} = this; let fprime = normalize(f); let rprime = normalize(orthogonalize(r, fprime)); let uprime = normalize(orthogonalize(
orthogonalize(u,fprime),rprime)); let aprime = normalize(orthogonalize(
orthogonalize(a,fprime),rprime),uprime)); this.f = fprime; this.r = rprime; this.u = uprime; this.a = aprime; };

The fact that this works may not be completely obvious at first glance. After all, how do we know that orthogonalizing the basis vector for the U axis against rprime won't mess up the previous orthogonalization against fprime? Well, we know that rprime is guaranteed to be orthogonal to fprime, because we constructed it to be so. That means that the projection of u onto rprime will also be orthogonal to fprime. Adding or subtracting a component orthogonal to fprime clearly cannot introduce any components to uprime that are are parallel to fprime, and so we can conclude that orthogonalizing against rprime, or any vector which is also orthogonal to fprime, cannot cause uprime to become non-orthogonal to fprime. The same reasoning can be extended to prove the validity of the orthogonalization procedure for the A axis, and this procedure trivially generalizes to arbitrary higher dimensions as well.

As with the rotation itself, it is possible to gain some efficiency here at the cost of some code duplication, but since this operation need be performed only very rarely it has yet to become a bottleneck in my code.

Next: Building a 4D World