I’ve been thinking lately about what makes good engineering. Actually, let’s be honest. As engineers, that’s not what we focus on. We get hung up on the failures—which always seem to outnumber the successes. Of course once you get out into the real world of programming for a livelihood, success is measured in terms of project completion and client satisfaction, but as easy as it is to put your head down and focus on the next ticket, that doesn’t absolve us of having to master our craft, to learn new things and better understand and utilize the old.
But I won’t try to say anything intelligent on good engineering. Instead, I’ll try to communicate the regrets that surfaced last weekend when I plugged in my old laptop and sifted through a decade-old folder called projects.
How I started programming
The first thing I went looking for was one of my first major programming efforts, the result of the reason I started programming. A pretty long time ago, I came across a paper called Visual Simulation of Smoke by some pretty impressive people at Stanford. Research like this is what makes those Pixar movies possible. That you can simulate such complex dynamics from a set of simple principles fascinated me. This was about halfway through high school, so I decided to learn vector calculus and linear algebra just to make this happen.
It was a long, slow process. Timestamps indicate that I finally got something working five years later in about 2005 (Junior year of college). Frankly, I still feel pretty good about the result!
+The finished product, five years in the making.
Down the Rabbit Hole
In hindsight, this shouldn’t really amount to more than a week or two of work, but learning is the process of making hindsight obvious. One way or another though, I had something special: knowledge. And experience. At this point, I had three options:
Learn from what I’d done and improve
Learn from what others have done and improve
Go down the rabbit hole and see how complicated I can make things
+Given my experience at the time, it’s hard to say I chose the wrong option, but let’s be clear: I chose (3), and I didn’t just slowly wade into it. I sought out complexity. I enhanced the complexity. No, seriously. I took the basic method and made it as absolutely opaque as I could. I wanted to enter it into the International Obfuscated C Code Contest:
This has basically all the features of the video above, except it’s dramatically less useful. I even RLE-encoded the University of Michigan ‘M’ on line 22 and a small part of line 21. It has a X-Windows interface. It doesn’t run correctly unless line 1 starts with a tab.
So am I proud of this code? Yeah, a little. (It’s on github by the way.) But that’s vastly outweighed by the regret I have over not having used this period of my life more effectively. It’s likely I may never again have the combination of time, motivation, and inspiration that it takes to dive so deep into this particular confluence of math, science, and programming.
Still more
Actually, I went a bit farther. The next stop on this journey through my time capsule was into the pejoratively named field of Colorful Fluid Dynamics (CFD) (not to be confused with the more useful field of Computational Fluid Dynamics). Colorful Fluid Dynamics is what you get when you leave behind rigor and accuracy and just try to make pretty pictures.
In 2006, I got my Incomplete LU-Preconditioned Bi-Conjugate Gradient Stabilized (ILU-Preconditioned BiCG-Stab) solver working and used it to solve Laplace’s equation (think: stretched membrane) on a two-dimensional finite volume grid. I’m pretty certain it remains the most complicated thing I’ve ever programmed successfully, but by the time I finally got it working, I was so exhausted and burned out on my search for complexity that the rest of my career in Aerospace Engineering was really just a slow coming-to-terms with the fact that I no longer wanted anything to do with it.
I found the complexity I was looking for, but it didn’t have the meaning I’d hoped for.
So did I learn anything?
I’m not trying to bash all my old projects. I learned a lot. I made some things I’m proud of. But for such a large amount of effort, there are so many ways I could’ve put this effort and energy to better use. So has anything changed? I’d definitely say so—at least I like to think I’m finally failing less—but for someone who codes full-time now, it’s hard to look at my github page and convince myself that I’m actually doing a better job of integrating meaningfully into a larger community of developers. Of course to be frank, that’s not the goal I’m after anyway. Tools are fun to use, learn, and create, but they’re not inherently meaningful. They are means to an end. So my real goal?
goal:
communicate passion and knowledge about things I find meaningful
At Standard Cyborg, we use 3D scans from a number of different sources to design sockets for prosthetic devices (btw we’re always hiring if this stuff interests you). Even in the best of circumstances in which we’re able to keep track of the physical “up” direction from acquisition through to the design phase, “physical up” isn’t really a useful direction since the limb could have been in any orientation during scanning (or any orientation relative to the scanner). We’d like to automatically orient incoming scans vertically so that they’re easy to work with.
+ Not an actual scan of a knee, but the below the knee (BK) scans users process tend to look very similar and tend to enter the system in no particular orientation.
+
Working with 3D models isn’t exactly new territory. Before throwing math at this, we should stop to consider whether a nice arcball camera (or—shudder—x/y/z rotation handles!) would allow users to orient scans as they see fit, removing orientation as any concern of ours. We know so much about this problem though! In a broad sense, we know what the scans look like and how users will be modifying them, and we know that even the most grizzled power users find extra degrees of rotational freedom cumbersome and frustrating when they’re not required for the task at hand. At the very worst, we find automatic alignment a great preprocessing step that helps users, doesn’t hurt what’s already arbitrary, and in many cases nails it right away.
The question remains then what alignment could possibly mean. There are an infinite number of valid meanings and corresponding solutions. The solution I describe here only addresses a particular meaning that happens to solve our little micro-problem quite well. It’s not new or novel—an alternate title for this article was “In which I discover the ellipsoid!”—and I’m only taking the trouble to describe it because I was so delighted to pick a heuristic out of the sky, find cause to break out some math, and actually end up with a function which runs robustly in a couple milliseconds.
Choosing an axis
Most of the scans for which people use our software (and limbs in general, really) are basically cylindrical tubes so that the longest axis makes a decent first cut for orienting the scan.
+ Orienting a scan by the longest axis isn’t entirely without merit, but it’s not very robust.
+
This orientation seems friendlier than what we started with, but it doesn’t take long to spot some problems. For one, we haven’t said anything about how to actually compute the longest axis (Principal Component Analysis (PCA) feels relevant?). More importantly though, if the scan had roughly equal proportions, the longest axis would be entirely arbitrary even for a cylinder with plainly obvious orientation.
length=0.3 width=0.7 depth=0.5
+ We’d like it if our algorithm always picked the axis of the cylinder, but the longest axis doesn’t accomplish that.
+
Failure of this basic sanity check suggests the overall orientation of these scans isn’t so much defined by the position of the surface as by its orientation. Without agonizing over why, I decided a better option would be to select an alignment axis as perpendicular as possible to the surface normals. Hazarding a guess at stating that mathematically, I’d call that the axis which minimizes the sum of the squares of the dot products of the alignment axis with the surface normal vectors.
+ As an educated guess, we instead aim to find an alignment axis as perpendicular as possible to the surface normal vectors.
+
As perpendicular as possible
The statement above is a mouthful which requires a bit of unpacking. If we’re going to tackle this as a minimization problem, we at least know we’ll need to roll up the ideas above into an objective function.
Let’s start with the dot products. Recall the dot product between vectors a and b is equal to a⋅b=abcosθ, where a and b are the magnitudes of the two vectors, respectively, and θ is the angle between them. All we really need to know here is that if two vectors are perpendicular, their dot product is zero.
We can talk about a single surface normal vector, but somehow we need to aggregate information across all faces. Let’s call a surface normal vector of the ith mesh face ni and a candidate axis of alignment ξ (the Greek letter “xi”, pronounced ”ksee″, which I’m selecting because it’s fun to write, isn’t likely to get confused with anything, and is fun to call “tornado” instead). My supposition is that if we dot the two, square the result, sum over the faces and call it f(ξ), i.e. f(ξ)=i=1∑nfaces(ξ⋅ni)2, then the best alignment is the one which minimizes f(ξ).
We might have a reasonable objective function here, but to see why it feels like it should work, consider a cylinder. The axis of the cylinder is always perpendicular to the surface normal vectors. Assuming for simplicity that the vectors are all normalized, then the magnitudes drop out and f(ξ)=∑cos2(90∘)=∑0=0. The axis of a cylinder minimizes f(ξ) even when it’s not the longest axis, thus fixing the failed sanity check above. (If you want to be fancy, I think you could say we’re solving the same principal axis problem but in the tangent space instead, though I don’t think that interpretation is likely to help most people.)
(Why the square? On a strictly mathematical basis, the dot product may be either positive or negative which would cause the minimization to diverge to −∞. The square keeps f(ξ) non-negative so that we can meaningfully minimize it.)
+ The surface normal vectors of a cylinder are everywhere perpendicular to its axis so that f(ξ)=0. A cylinder passes the test!
A bit more precisely, if the faces comprising the mesh aren’t uniformly distributed, the sum will be biased toward clusters of vertices and their associated normals. Instead of a sum over normal vectors ni, what we really want is an area-weighted sum. In fact what we really want is just an integral over the surface (call it S) with respect to the differential area vector (call it dA). We define dA≡ndA as parallel to the surface normal but with magnitude equal to the area of a differential surface element. The continuous limit of f(ξ) is then f(ξ)=∫S(ξ⋅dA)2.
While we’re being precise, we assumed implicity that the axis of alignment was a nonzero vector, but let’s now make that explicit in order to avoid the trivial solution ξ=0=(0,0,0) which always minimizes f(ξ)=0. Constraining ξ to be a unit vector will do just fine.
Fully stating our problem, we want to find the argument ξ which minimizes f(ξ) subject to the constraint that ξ is a unit vector: 0argminξ∈R3∫S(ξ⋅dA)20subjectto∣ξ∣=1.
For piecewise constant faces with surface normal dAi (magnitude equal to the face’s area, recall), we can recast this as a discrete summation and arrive at our final problem statement, 0argminξ∈R3i=1∑nfaces(ξxdAi,x+ξydAi,y+ξzdAi,z)20subjecttoξx2+ξy2+ξz2=1.
As for the areas, Eric Arnebäck has a nice article about Computing the Area of a Convex Polygon. It covers triangles. And for you geometry sorcerers and sorceresses, the answer is yes. We’re fitting an ellipsoid now. The rest of the article is me realizing I’m looking for an ellipsoid.
Computing it
The problem above is a constrained optimization problem. Those can be a bit challenging to solve since you often only want to explore the solution space in directions which keep the constraints satisfied. It took me a while to recall, but if I learned one thing about constrained optimization in engineering (sadly I didn’t learn much more), I learned that the method of Lagrange multipliers exists to transform constrained optimization problems into unconstrained problems. The method works like this. Instead solving the problem minimizef(ξ)subjecttog(ξ)=0, we solve the problem minimizeL(ξ,λ)=f(ξ)−λ⋅g(ξ) where λ is an auxiliary parameter (the “Lagrange multiplier”) that drives objective function toward satisfying the constraint. With just a bit of handwaving, we can demonstrate that setting the partial derivatives of L equal to zero yields 0=∂λ∂L(ξ,λ)=−g(ξ)0=g(ξ) which confirms the constraint is satisfied, and 0=∇L(ξ,λ)=∇f(ξ)−λ∇g(ξ)∇f(ξ)=λ∇g(ξ)=0, with the final leap of faith equality to zero taken since g(ξ)=0 is a stationary point. This step then enforces the original objective function, though I haven’t adequately justified it here. Wikipedia actually has a pretty good explanation which I’d be foolish to try to outdo.
It only takes the tiniest modification to state our problem in the canonical form of a Lagrange-multiplier-ready problem, 0argminξ∈R3i=1∑nfaces(ξxdAi,x+ξydAi,y+ξzdAi,z)20subjecttog(ξ)=ξx2+ξy2+ξz2−1=0. Applying the method, we arrive at the unconstrained problem 00argminξ∈R3λ∈Ri=1∑nfaces(ξxdAi,x+ξydAi,y+ξzdAi,z)2−λ(ξx2+ξy2+ξz2−1).
Taking the partial derivatives with respect to ξx, ξy, and ξz as well as λ and equating to zero isn’t particularly tedious. The result is a system of four simultaneous equations, ⎩⎪⎪⎨⎪⎪⎧0=ξx∑dAi,xdAi,x+ξy∑dAi,ydAi,x+ξz∑dAi,zdAi,x−λξx0=ξx∑dAi,xdAi,y+ξy∑dAi,ydAi,y+ξz∑dAi,zdAi,y−λξy0=ξx∑dAi,xdAi,z+ξy∑dAi,ydAi,z+ξz∑dAi,zdAi,z−λξz0=ξx2+ξy2+ξz2−1
It suddenly feels hopeless, especially since the fourth equation is a bit nonlinear in ξ. Let’s cut down on the visual noise by defining Axy=i=1∑nfacesdAi,xdAi,y as well as the analogous definitions for all pairwise combinations of axes. With these definitions, the above equation looks a bit more manageable, yielding ⎩⎪⎪⎨⎪⎪⎧0=ξxAxx+ξyAyx+ξzAzx−λξx0=ξxAxy+ξyAyy+ξzAzy−λξy0=ξxAxz+ξyAyz+ξzAzz−λξz0=ξx2+ξy2+ξz2−1.
Neglecting the last equation for a moment, we can state the first three as a matrix multiplication, ⎣⎡AxxAxyAxzAyxAyyAyzAzxAzyAzz⎦⎤⎣⎡ξxξyξz⎦⎤=λ⎣⎡ξxξyξz⎦⎤.
This is just the standard form of an eigenvalue problem, Aξ=λξ, and what’s more, its eigenvectors ξ are normalized by convention, which implicitly satisfies the constraint ∣ξ∣=1. Eigenvalues are simple and easy to compute, even in JavaScript. We’ve solved it! Upon solving, we get three eigenvalues and corresponding unit eigenvectors which are identically the model axes and associated inverse strengths along the respective eigenvectors.
As a final bonus, recall—or discover today!—that the eigenvalues of a symmetric positive-definite matrix are real and orthogonal, i.e. mutually perpendicular. And there are three of them. So we don’t just get unit vectors out of this, we get a three dimensional rotation matrix which can be applied directly to the model.
Does it work?
Wonderfully! Robustly! Efficiently! The only nontrivial numerical part is the eigenvalue computation, but it’s only a small 3x3 matrix you can farm out to any old numerial library.
The main caveat is that eigenvalues are only unique up to a sign so that we need to check for reflections and apply some slightly ad-hoc heuristics to disambiguate the sign. In particular, I’m just using the total summed area vector to see if we can put the open end in a consistent direction. There’s room for improvement.
You can see the final result below. Note that the two remaining axes also align the knee!
Update: Eric Arnebäck asked about noise. I’ve added a noise slider below and have removed a square root in the scaling so that the magnitudes are a bit more separated. The noise is not IID noise so take it with a grain of salt, but it hopefully gives some indication of the approach’s ability to reject noise.
curvature=0.5 radius=1 length=1 noise=0.1
+ Adjust the sliders and observe the effect it has upon the computed alignment.
+
Conclusions
At the end of the day, I rather suspect I’ve rederived a pretty standard technique for talking about the shape of a surface. I hope you’ll forgive me if my satisfaction isn’t diminished though since opportunities to legitimately break out Lagrange multipliers are so rare! And as part of my day job no less.
There’s room for improvement in the final disambiguation of signs, but frankly once we’ve solved the main problem of figuring out a rough alignment, the subsequent algorithms have a significantly easier time making sense of the scan.
This post uses idyll and regl. They’re great projects! You should check them out! You can find the article source here and an implementation of the algorithm here.
Questions? Comments? Corrections? Drop me a line @rickyreusser!
