Rotations with Arcball

So after trying to figure out how to rotate a 3D world by clicking on a 2D object, I think I found a decent implementation. There are bits and pieces all over the Web about Arcball, but there wasn't one in WebGL, so I decided to make one using glMatrix (http://github.com/toji/gl-matrix). I learned lots about 3D programming, the pitfalls and the possibilities. I also learned numerical stability is a bitch. Since I couldn't find a tutorial on the Web that goes from A to B with Arcball, I'm going to write one.

First off, the way to picture the Arcball is a sphere centered in the middle of your viewing space. You're looking down at it from the top, so you only see the outer surface. The way to make this happen is:

1. Normalize the screen coordinates so that the X and Y coordinates are in [-1,1] (what I'm going to call view coordinates) and that the X and Y directions are to the right and up, respectively. (Usually in windows, the direction with increasing Y might be flipped).
2. We want a sphere with radius 1, so that means that sqrt(x*x + y*y + z*z) = 1. Thus, z = sqrt(1 - x*x - y*y). Make sure this vector is normalized so the click point maps to the sphere's surface.

Now we have a map from a 2D screen click to a 3D unit vector in view coordinates. Let's get on to rotation.

1. When we click for the first time, let's call that vector U = (u1, u2, u3).
2. If we want to change the rotation on mouse drag, then every time we move the mouse, we need to capture a drag vector V = (v1, v2, v3). 
3. To get the rotation axis W for these vectors, we can define W = U x V, the cross product.
4. To get the rotation angle a, we can define a = acos(U . V), the dot product.
1. Since the dot product might give us a value in the domain of Math.acos, we can use the min(1.0, U . V) function to keep it in the bounds.

Note: W and a are in view coordinates still. In this post, I'm treating them as if they were projected into modelview coordinates. To do this "the right way", we would need to project them into the model by a method like gluUnproject or multiplying by the inverse of the combined projection and modelview matrices. Since I'm not that far yet with my implementation, I'm just going to assume for the rest of the post that my W vector is in modelview coordinates and that it passes through the modelview origin.

Now we get into the heart of this post. In previous posts, I detailed rotation in 3D in 2 ways: with matrices and with quaternions. Now, I said that matrix multiplication isn't as numerically stable as quaternions. I learned that one the hard way. When doing the 7 matrix multiplies, I got some weird values and I couldn't hunt the problem down. Plus, it wasn't as easy as quaternions. So, I used the quaternion method I detailed below to do all my rotation and translation. But, there is one trick I did that allowed me to get a good, proper rotation matrix.

The trick is that I simply rotated the X, Y and Z unit vectors around W. If the quaternion rotation is accurate (which it is), then the vectors should stay orthogonal to each other, so no weird scaling and skewing. I constructed a rotation matrix from these vectors as follows:

R = [a[0] a[1] a[2] 0; b[0] b[1] b[2] 0; c[0] c[1] c[2] 0; 0 0 0 1]

This seems to give me what I need. I know the 4th column is related to the "translation" of the origin, but I'm not there yet. Hope this helps anyone.

Cheers!

Rotations in 3D with Quaternions

This is part 2 of the 3D rotations post. After looking at rotations with matrices, let's now consider rotation with quaternions.

The real problem with matrix representations of rotations are the number of calculations in multiplication and the fact that they're not very numerically stable. As well, a spherical linear interpolation for smooth rotations (SLERP) isn't easily done. Quaternions are good for these reasons. Plus, it saves on the amount of information that needs to be stored.

Quaternions have "real" and "imaginary" parts. Although quaternions are not vectors, a lot of the same operations and terminology are shared between them. For example, a quaternion can be represented as:

Q = a + bi + cj + dk,

where a is the "real" (scalar) part and b, c, d are the "imaginary" coefficients. Quaternions can easily represent axis-angle rotations. The angle (in radians) would be the scalar and the vector (axis) would be the "imaginary" part. This is precisely why quaternions are so useful.

On other websites, you can see how quaternion math works, so I would post it here, but the steps to working with them will be shown. One other thing is to note that in naive matrix multiplication takes 128 FLOPs (2n^3) for 2 4x4 matrices. Multiplying 2 quaternions takes 28 FLOPs, which is much faster (although matrix multiplication algorithms like Strassen's can decrease the matrix multiplication time significantly).

So...on to the technique. The steps are:

1. Create a "translation" quaternion by putting 0 in the scalar part and the "translation" vector for the rotation line back to the origin in the imaginary part. If the rotation line vector v is defined from a point v0 to v1, then the quaternion would be T = 0 - v0.
2. Calculate the conjugate of T, which is just TConj = v0.
3. Create a quaternion for the point you want to rotate: P = 0 + pxi + pyj + pzk.
4. Create a quaternion for the rotation axis v and angle t and calculate its inverse: R = cos(t/2) + vx*sin(t/2)i + vy*sin(t/2)j + vz*sin(t/2)k.
5. Translate the point by adding the translation quaternion: A = P + T
6. Rotate the point by multiplying the new point by the rotation quaternion and its inverse: B = R*A*RInv
7. Translate the point back by adding the conjugate of the translation: P' = B + TConj

So the entire equation then becomes P' = (R * (P + T) * RInv) + TConj, where P' = S + x'i + y'j + z'k and x', y', z' are the coordinates of the rotated point. I think the scalar part might map to the point weight in homogeneous coordinates, but I'm not sure.

Anyway, there it is. I hope to learn about SLERP in the meantime and possibly post something up later about it. Cheers!

Rotations in 3D with Matrices

3D Rotation has always been a problem for me because I never really sat down and worked out the math to the point where I can be confident that it would work. Every page I read it on would confuse me and talk about rotations in general, gimbal lock, quaternions, Euler angles, etc. They focused more on the math concepts and less about the implementation. This post is all about just getting it done.

Quick summary: 3D programming uses matrices to turn the objects in the world. In order to rotate objects in OpenGL, for example, we need to manipulate the modelview matrix. Although I'm not all that knowledgeable about OpenGL, websites I have read have said that altering the projection matrix will mess with the lighting and things. I will show you a way to manipulate the matrix with matrices first, then with quaternions in a later post.

Usually, rotation about any line is accomplished by multiplying 7 4x4 matrices together to get the right matrix. To be clear, this matrix will only rotate points. If you are going to rotate vectors or lines, you have to rotate the definition points. The 7 matrices needed to rotate any vector about any line are:

1. A translation matrix to move the rotation line so that it passes through the origin.
2. A rotation matrix to rotate the line (axis) about the X-axis to the XZ-plane.
3. A rotation matrix to rotate the line (axis) about the Y-axis to the Z-axis.
4. A rotation matrix to rotate the desired angle about the Z-axis.
5. The inverse of the matrix in step (3).
6. The inverse of the matrix in step (2).
7. The inverse of the matrix in step (1).

Step (1) is only needed if the line you want to rotate about doesn't already pass through the origin. You can see that the rotation just moves the line to the Z-axis and does the rotation about that. The thing about these matrices is that we can avoid using trigonometric functions of the angles by repeatedly setting up the following equation, solving for the cosine and sine terms, and then substituting them back into the matrix:

[1, 0,       0,        0;  * [x;  = [xp;

 0, cos(t), -sin(t), 0;     y;      yp;

 0, -sin(t), cos(t), 0;     z;      zp;

 0, 0,        0,       1]    1]       1]

The translation matrix is the equivalent of adding subtracting a vector from the origin to one of the points on the line. The "inverse" is the addition of this vector back to the point. This code assumes the rotation line starts at v0 and ends at v1, and the line itself is the vector v. The translation matrix and its "inverse" is as follows (in Matlab code):

T = [1,0,0,-v0(1);

        0,1,0,-v0(2);

        0,0,1,-v0(3);

        0,0,0,1];

Tinv =[1,0,0,v0(1);

            0,1,0,v0(2);

            0,0,1,v0(3);

            0,0,0,1];

The code for the rotation matrices are:

Rx = [1, 0,                           0,                          0;

          0, v(3) / sqrt(v(2)^2+v(3)^2), -v(2) / sqrt(v(2)^2+v(3)^2), 0;

          0, v(2) / sqrt(v(2)^2+v(3)^2),  v(3) / sqrt(v(2)^2+v(3)^2), 0;

          0, 0,                           0,                          1];

Ry = [sqrt(v(2)^2+v(3)^2), 0, -v(1),               0;

          0,                   1, 0,                   0;

          v(1),                0, sqrt(v(2)^2+v(3)^2), 0;

          0,                   0, 0,                   1];

Rz = [cos(t), -sin(t), 0, 0;

          sin(t),  cos(t), 0, 0;

          0,       0,      1, 0;

          0,       0,      0, 1];

Note: The Rx matrix depends on the term sqrt(v(2)^2+v(3)^2). This term vanishes when the line is coincident to the X-axis. If this is the case, then the Rx matrix will need to be the identity matrix.

Rotating the point consists of multiply the matrices to the left side of the point. Since rotation matrices are orthogonal matrices, the inverse of the matrix is the transpose. Transforming a point P to P' would be:

P' = Tinv * RxInv * RyInv * Rz * Ry * Rx * T * P

That's it! Happy rotating!

2D Convex Hull: The Graham Scan

Although the Jarvis march is a fairly good algorithm, the Graham scan is somewhat better; it only scans through the entire set of points once, although it runs through several of those points more than once. The algorithm steps are as follows:

1. Find the point with the lowest y-coordinate. This is the extreme point we start with and call it P (like in the example).. This takes O(n) time.
2. Calculate the polar angles between the x-axis and the line through P and each other point. 
3. Sort the point set by ascending angles and include P last. This operation takes O(n log n) time and any sorting algorithm can be used.
4. Start at P and check the next point if it makes a "left turn", or if the cross product is positive. The cross product should be made up of 3 points, like (A-P) x (B-A). If it is, include the point in the convex hull (B). If not, remove the last point (A) in the convex hull and add the current point.
5. If a "right turn" is found, the points in the convex hull should be worked backwards until a "left turn" is found. This is an important step. The picture doesn't show it because the point set doesn't require it.
6. Continue the march until the current point is P.

This algorithm is slightly better since it runs in O(n log n) time. This isn't always better than the Jarvis march, which runs in O(nh) time. The Graham scan can be slower if the number of points on the hull h < log(n), however it isn't output-dependent.

This is a bit more complicated than the Jarvis march, but still somewhat simple. The advantage is its linearithmic time complexity. You can be fairly sure that the algorithm's time to completion won't blow up too large. The disadvantage is the polar angle comparisons. Although you can use the definition of the dot product and use just the cosine of the angle (since it's monotonically decreasing), the difference between the angles can be very small and error can be significant. Another disadvantage is that it can't be parallelized very easily.

There is another convex hull algorithm that can be parallelized better: the divide-and-conquer algorithm.  I'll cover that in another post.

The Convex Hull and Jarvis March

The convex hull of a set of points is important because it's related to the Delaunay Triangulation. According to Wikipedia, the convex hull of points in R^3 can be projected onto a plane as a DT in R^2 and a convex hull in R^4 can be projected as a DT in R^3. With this in mind, I've tried to study convex hull algorithms. There are a lot of algorithms that compute the convex hull and they all have their advantages and disadvantages. I will attempt to explain what I've discovered about some of their algorithms.

The Jarvis March (Gift wrapping algorithm)

This is probably the simplest algorithm to get the convex hull. We start with a set of points, like in the picture. The algorithm steps as follows:

1. Find an extreme point that is guaranteed to be on the convex hull. Many times the point can be the minimum or maximum value in a particular direction. For example, we can choose the point with the minimum value on the x-axis (not shown in the picture). 
2. For the first point, calculate the polar angle between the x-axis and the line connecting the extreme point and every other point in the set and the point that minimizes this angle is the next point on the convex hull.
3. For every other point, calculate the polar angle between the line connecting the last two points added to the convex hull and the line connecting the last point added to the hull and every other point in the set. Again, the point that minimizes this angle (without it being zero) is the next point in the convex hull.
4. Step 3 is repeated until it encounters the first extreme point, when the algorithm ends.

It's really a march around the set to find the convex hull. R.A. Jarvis came up with the algorithm in 1973. The advantage to the Jarvis march is that it's very simple. Simple to understand and simple to construct. However, the major disadvantage is related to its time complexity. Time complexity is a way of comparing the time an algorithm will take to the size of the input. It's expressed in big O notation, a way to describe a function's behavior in terms of other functions. The Jarvis march is output-sensitive, meaning that the time to completion depends on output as well as input. In this case, the Jarvis march runs in O(nh) time, where n is the total number of points and h is the number of points on the convex hull. You can see this easily because for every point on the hull, it checks all n points.

We often need to check on the worst-case scenario. For this specific algorithm in 2D, let's say all of the points are arranged in a perfect circle. That would mean that all of the points would be on the convex hull, or n = h. That would mean that the worst-case would be T(n) = O(n^2). Also this isn't necessarily bad, there would certainly be a way to improve on the algorithm. Others have thought this too and created better, more efficient algorithms. I'll cover that in another post.

Delaunay Triangulations

Apparently, Delaunay triangulation is harder than I thought. It seems like such a simple idea: draw triangles between points in a set so that the circle that runs through the triangle's points doesn't have any other points in it. Although there are apparently lots of solutions, but the simpler ones that tutorials exist for aren't very efficient and no one's talking about the efficient ones.

A 2D Delaunay Triangulation in MATLAB

Why Delaunay? Why can't I just be happy with drawing regular triangles? It's because the DT focuses on maximizing the internal angles. In layman's terms, it avoids long, skinny triangles which turns out to be bad for both rendering and finite element analysis. This will help when I eventually have to split surfaces up into triangles for OpenGL rendering. Of course, there might be a quick easy way to tessellate a Bezier surface without resorting to triangulation methods. If there are, please let me know ASAP!

Anyway, I'm trying to first start out with a 2D triangulation. It seemed like the simplest way to go. I implemented a Bowyer-Watson algorithm because it was easily explained and it was incremental. It is also supposedly one of the fastest DT algorithms out there. Here's the basic gist of how it works:

1. If you have a set of points, start out by drawing a "supertriangle" around all the points. The DT should just consist of the "supervertices" right now, which don't really belong there but we'll remove them later. Since the triangulation consists of only 3 points right now, you can say it's Delaunay.
2. Include each point from the point set into the DT. Check to see which triangles (by checking their circumcircles) are made non-Delaunay by the point.
3. Remove the triangles so you have a void in the DT.
4. Take all the points on the edges of that void and make triangles that include the new point as a vertex.
5. Repeat this until you've added all the points from the original set.
6. Now remove all triangles that have any "supervertices".

If you've followed this algorithm, then supposedly you should have a Delaunay triangulation of the point set. I made code that did this, but it wasn't very efficient. After hours of trial and error, it finally produced what it was supposed to! Well, kinda. I did get a Delaunay triangulation, but removing the "supervertices" from the set caused the DT to not represent the true convex hull of the points. It removes the triangles that contain points that, if connected, would build a convex hull. I guess the question is, "If you know what the problem is, why not fix it?" Good question. First of all, I don't know how to just limit my searches to triangles around the outer edge. I read about a quad-edge data structure that Guibas and Stolfi came up with to map the DT topologically (meaning each triangle would know who its neighbors, either triangles or points, were) and I could use that, but I didn't quite understand how to pinpoint only the DT edges and how to connect them or check to see if they were already connected. Second, I read that the "supertriangle" method doesn't always produce correct results, which worries me because I don't know any other method for producing a starting DT. This algorithm would be great if I already had a rough DT, but I haven't seen a good method for doing that either.

I'm sort of stuck on Delaunay triangulations, but I read about a technique that talked about a 2D DT being a projection of a 3D convex hull of a paraboloid onto a plane. I know I used some big words, but let me explain. A paraboloid is a 3D representation of a parabola. Imagine taking a normal parabola (like y= x^2) and revolving it about the axis that passes through the saddle and the focus (like x=0) and you get a paraboloid. For a set of 2D points in (x,y) format, we can find the "middle" of the point set (center of mass) and set that as (xm, ym). We then set the Z value of the 2D points equal to sqrt((x-xm)^2+(y-ym)^2) and that way we get the paraboloid. Then, we just compute the lower convex hull of the paraboloid and then we can project it onto a 2D plane. This will create a DT of the 2D point set.

I guess my next task is to learn how to get a convex hull in 3D. However, I think I'll start slow and learn how to do it in 2D. It seems there are a lot of techniques and I'll cover some of them in my next post. Happy sciencing!

Visions of Flight to Virtual Worlds

Since last summer, a lot of things have happened:

• I've hung up my wings. 
• I help develop some software that got me published in a mechanical engineering journal.
• I completed a project for Boeing in school and it might land me a job there.
• I finally finished my bachelor's degree.
• I got accepted to graduate school at BYU.

About giving up flying, I tried to fly one last time after my last published failures but I figured out I have a bad servo on one of the wings and with all my bad experiences so far, I'm done with flight for a while. I'm content to solve problems here on the ground. Since I hung my airplane from my light fixture last summer, I took a class on computational geometry just before I graduated. It focused mainly around curves and surfaces used in computer-aided engineering and the mathematics behind them. Dr. Thomas Sederberg was my teacher and he seems to have been involved quite heavily in computer-aided geometric design (CAGD). He really sparked a huge interest I've had in applied mathematics. I've always wanted to do more with mathematics, but I was never convinced that a math degree would get me a decent-paying job. All the math majors I know became teachers, including my cousin. CAGD, however, seems like it combines the best of math with computer science to come up with a product perfectly suited for engineering. In other words, I found my Holy Grail.

I didn't learn about OpenGL or any 3D programming, but I did learn about Bezier and NURBS curves and surfaces. I wish I learned more about the inner workings of CAD and animation systems so I could construct a visualization module, which is my next project. This one I might even finish. I'm going through NeHe's OpenGL tutorial and I've learned a lot. I'm still unsure about some things, but I'm planning to take a class in graduate school on computer graphics. Apparently, Utah has been a hotspot for computer graphics so there would likely be some great resources at BYU. But until I get a better handle on how 3D programming works, I've decided to learn about surface tessellation. I've tried starting out simple and learning about Delaunay triangulation. I'm going to post more about my efforts to learn about how this works. I've been searching the Web using every technique I know in order learn more about it, but there are no really good tutorials. Maybe my trail of breadcrumbs will inspire someone to write some.

Anyway, I hope I stay more on task with CAGD than I did with flying.