Fit Anything

Perhaps every programmer has a collection of hammers that they use to hit the nails they encounter all the time. Some hammers are more powerful, but it is always good to have more of them handy. I almost lost one of my hammers, having not used it for a while (and the code being lost due to the unfortunate demise of my last laptop). Fortunately, I recently picked it back up (along with a even more powerful version of it), when my colleague Yufeng presented the following question:

"Given any polygon, how to determine how similar it is to a triangle?" (He is trying to find all the squares in NYC that are in fact triangular)

Along with the problem, he also proposed a simple solution: to run polygon simplification (more specifically, Douglas-Peuker, or DP) on the convex hull iteratively until there're only 3 vertices left: compare the area of the resultant triangle with that of the initial polygon -- the closer they are, the more the initial shape resembles a triangle. The algorithm worked reasonably well on many cases, but he was wondering if there's a more sophisticated one.

Initially we're thinking about whether the ways of computing an oriented bounding box (OBB) could be easily generalized to "bounding triangles" too. Meanwhile, I started to vaguely recall that I might have solved similar problems in the past, and the keyword is "Least Squares". However, what I didn't manage to recall was how exactly to formulate the problem in this case. So I started waving my hands and saying something like "make this matrix with the inputs, and solve for this other matrix of parameters such that when you multiplied the two you get the matrix of the outputs..."

With my curiosity unsatiated, I went home and started texting another friend of mine (ChatGPT) to jog my memory. According to them, what I really need instead is "Non-Linear Least Squares", or NLLS. While we were discussing the implementation details, I gradually started to remember how I have been solving this type of problem without requiring the "non-linear" part, i.e. linear least squares. I ended up implementing both, and it turned out that either of them has its own merits and advantages over the other. Let me start by describing the old algorithm that I previously came up with.

Given an input polygon, resample it to have exactly N vertices (say, 100). Make a "template" polygon (the shape you want to fit to, and in this case, a triangle), also resample it to 100 vertices. Now, solve for the best transformation matrix that, when multiplied to the 100 points of the template, gives 100 points that are closest to the 100 points of the input polygon. That's it: the transformed template is the optimal fitting, and the sum-of-squares residual is the similarity score for the original question. Let's explain some of the concepts.

A transformation matrix is a bunch of numbers packed in a way that it encodes a transformation, such as rotation, translation, scale, distortion... or any combination of some such. When you "apply" the transformation matrix by multiplying it to a point, it rotates, translates, scales or otherwise "transform" the point. In our case of triangles, a 2D "affine" transformation matrix, which consists of 6 parameters, is sufficient to encode anything we can possibly do to the shape.

To formulate the least squares problem, which can be expressed as solving for "x" in "A * x = v", we "unroll" the affine matrix application into two rows of "A" and "v". e.g. [X[i] 0 Y[i] 0 0 1; 0 X[i] 0 Y[i] 0 1; ...] * [a; b; c; d; e; f] = [X'[i]; Y'[i]; ...], where (X[i],Y[i]) are from the template points, (X'[i],Y'[i]) are from the input points, and [a c e; b d f] is the affine matrix parameters we're looking for.

While this can be solved in one fell swoop of least squares algorithm, there're a few catches. Firstly, "correspondence" must be known: since this formula explicitly maps point 1 of the first shape to point 1 of the second shape, and point 2 of the first shape to point 2 of the second shape, and so on, this requires that the vertices of the input polygon be ordered (which it is in our case), and start at the same "phase" as the template polygon, to find the global optimal fitting. In other words, if the template is a triangle and the first point one of the three corners, the first point of the input polygon must also be on a corner (that is if the input shape already resembles a triangle, and if not, be on a certain point, which will end up best being matched to a corner).

While there're method like Iterative Closest Point (ICP) and Hungarian Algorithm that can deal with correspondence, in our constrained case the solution can be much simpler. To generate the optimal matching for arbitrary polygons with arbitrary "phase" in relation to the triangle template, we simply rotate the points of the input polygon N/3 times and try our algorithm at each possible configuration (and pick the best result). N/3 is the derived from the symmetry of the triangle; for squares you'll need N/4, for ellipses you'll need 1, for unknown shapes you'll need N-1.

The other problem is that if you're trying to match, say, a very skinny triangle to an equilateral triangle, the resampling will give one side of the skinny triangle way too few points, and the other two sides way too many points, whereas the equilateral triangle will have same number of points on each side. Therefore, some points will have to match with points on a neighboring side -- it will not make a great fit. This also applies to skinny shapes in general: some portions of it will get not enough resampling, whereas others get too much.

To eliminate this problem, we simply squash the skinny guy so that it is as fat as it is tall. We first find the oriented bounding box (OBB) of the shape with principle component analysis (PCA), and scale it in the direction of the eigen so that the aspect ratio is 1:1. Then we do the resampling, then we do the fitting, and finally we un-squash.

The main benefit of this least squares method is that it's simple, fast and robust. You might need to try a finite number of starting points, but the whole thing is deterministic and the result is guaranteed global optimum (within resolution).

Now let's look at the non-linear alternative: NLLS.

The core idea of NLLS solving algorithms roughly goes as follows: We're trying to find a set of optimal parameters such that we minimize some "cost" that is affected by these parameters. (In our particular case, the parameters is a representation of the triangle, and the cost is the dissimilarity between the triangle and our input shape). To do so we need to have an initial "guess" of the parameters, and a "cost" function to tell us how bad any guess is. At each iteration we check how off our current guess is by evaluating the cost function, obtaining the residuals "r". Then we check how the residuals will change if we "nudge" our guess ever so slightly in every possible direction and store it in the Jacobian "J", a matrix of first order partial derivatives. We solve for "delta" in "J * delta = -r" using linear least squares, to move our guess by delta. When delta approaches zero, it means that we're likely near a local optimum. This essentially describes the Gauss-Newton algorithm for NLLS; more sophisticated algorithms such as the Levenburg-Marquardt (which I end up implementing) adds some improvements (trust region, damping etc.) to prevent the algorithm from going crazy too often, but mostly follows the same idea.

To compute the Jacobian, some would use automatic (or manual) differentiation, but since we don't have that many parameters, we use a numerical technique called finite differences, which literally nudges the guess by a tiny epsilon and checks what happens. It's like standing on some terrain trying to descend but not knowing where to go, and trying to take a small step in every direction and seeing if that helps lower our altitude.

Initially I parameterized the triangle by its three vertices: 6 parameters of 3 (x,y) pairs. The cost is computed as the sum of the distances from each point to the closest point on the triangle (min of point-segment distances). While this worked well, I began thinking about how one could generalize it to shapes with more vertices than triangles (I imagined that there'll be nothing preventing the fitted shape from becoming a self-intersecting, entangled mess). It is then that I remembered my previous method of applying a transformation matrix to a template, and realized that the same parameterization can be adapted for NLLS too. It turned out that the 3x2 affine transformation parameterization is a bit "harder to crack" than the vertex parameterization and made the solution explode more often (understandably so). That's why I upgraded my Gauss-Newton algorithm to Levenburg-Marquardt, and added a bound and a regularization term (that favors non-crazy solutions to crazy solutions).

The remaining problem is that the NLLS algorithms find a local, not the global optimum. That is to say, the initial guess matters a lot. I mean, really a lot. If you just give it a default starting configuration without putting an effort into guessing, it will just get stuck in the middle of nowhere and be like "yeah I tried". If you know your shape beforehand, e.g. triangles in our case, it is not difficult to come up with a heuristic for the initial guess (e.g. bounding boxes, convex hulls, ...). But I also wanted to generalize it without making too many assumptions. Therefore, I simply make it try a bunch of starting configurations and pick the best one. In the end, I first normalize the translation and scale, and try eight non-collapsing versions of the upper-left 2x2 matrix with each element being one of (-1,0,1) (basically axis aligned rotations and flips). This end up working quite well. If we know the shape is symmetrical in some way, we can further reduce the number of starting configuration needed.

Let's take a step back and look up what we essentially have created: An algorithm to fit any shape to any shape. No matter what your shapes are, as long as it can be achieved through an affine transformation, we can fit it. But let's take even further step back, and marvel at how awesome NLLS is. With linear LS, you have to twiddle your problem to squeeze it into this special matrix layout, and many problems simply won't fit. NLLS doesn't care, you just give me some numbers, and tell me how bad they are, and I'll just make your numbers better! It literally can solve almost anything that needs solving!
At the same time, the "requiring good initial guess" part is a bit disappointing, for some problems at least. It's like saying, "if you give me 1 billion dollars, I'll give you back 1.25 billion dollars.". Sure it's great that it can make me a quarter billion dollars, but where do I get the 1 billion dollar to start, and, if I already have 1 billion in my pockets, why would I even bother. Of course, the numbers and proportions in the analogy would vary for different problems, and it'll probably be worth it most of the times -- I'm yet to find out in my future experiments.

Another interesting difference between my old resample+LS method and the NLLS method is that the latter doesn't require you to finish a shape to start giving good fits. If you draw, say an arc that is a quarter of a circle, the NLLS will happily finish the whole circle for you, whereas resample+LS will try its best to fit the arc to a circle -- the result is of course still reasonable (it does what you asked for), but not as "smart" as NLLS (figures out your heart's desires).

But hey! my method isn't completely useless. For many situations where NLLS would be overkill my simpler method would suffice; it also relies less on initial guesswork. Previously I've made quite a few interesting projects with it: I remember one time when another collegue showed me this "circle drawing game" that was going viral, where players compete to draw the most perfect circle. I was frustrated that it forces a fixed origin for you before you draw, and no matter how perfect your circle is by itself, if the center does not perfectly align with the fixed origin, you'll receive massive panelty. I was like "cool idea, but this is not the correct way to code that!" and went on to prove my point by writing a least squares solver for circles using this method. Unfortunately my version was not nearly as popular as the original game, perhaps due to the fact that I completely forgot to publish it. Just as I was (re)developing the algorithms for the triangle-fitting problem, my advisor prof. Zach Lieberman asked me to port some of his pedagogical examples to p5.js, and guess what? one of them is a perfect circle drawing game! It dynamically computes a centroid and checks standard deviation of distances -- so already massively better than the fixed-origin implementation while still being relatively simple. Zach asked me to showcase my more sophisticated algorithms during the class in case anyone is interested in fitting ellipses and other shapes in addition to circles.

When I tried ranting to my math-savvy friends about how incredible Levenburg-Marquardt is and how I only just learned about it, they were like "yeah that's a classic, that's what we've always been using". Oh well, better late than never! Now if you'll excuse me, I have a lot of newfound nails to hammer with my newfound hammer!