You have to convert the gdi+ image to a gdi+ bitmap and then to an OpenCV matrix. There is no easier way to do the first conversion than creating a bitmap and painting the image into it, as far as I know.

To perform the second conversion (Gdiplus::Bitmap -> cv::Mat), note that Mat constructors order their parameters rows then columns so that is Bitmap height then width. The actual memory layout is row major, however, so we can just copy the data, but there is no need to do the actual copy yourself. You can use one of the Mat constructors that will wrap existing data without copying and then force it to copy by calling the clone member function.

The other trouble here however is that gdi+ supports, in theory at least, loads of exotic pixel formats so to handle them all would be a chore. The following handles the basic case:

Gdiplus::Bitmap* GdiplusImageToBitmap(Gdiplus::Image* img, Gdiplus::Color bkgd = Gdiplus::Color::Transparent)
{
	int wd = img->GetWidth();
	int hgt = img->GetHeight();
	auto format = img->GetPixelFormat();
	Gdiplus::Bitmap* bmp = new Gdiplus::Bitmap(wd, hgt, format);

	if (bmp == nullptr)
		return nullptr; // this might happen if format is something exotic, not sure.

	auto g = std::unique_ptr<Gdiplus::Graphics>(Gdiplus::Graphics::FromImage(bmp));
	g->Clear(bkgd);
	g->DrawImage(img, 0, 0, wd, hgt);

	return bmp;
}

cv::Mat GdiPlusBitmapToOpenCvMat(Gdiplus::Bitmap* bmp)
{
	auto format = bmp->GetPixelFormat();
	if (format != PixelFormat24bppRGB)
		return cv::Mat();

	int wd = bmp->GetWidth();
	int hgt = bmp->GetHeight();
	Gdiplus::Rect rcLock(0, 0, wd, hgt);
	Gdiplus::BitmapData bmpData;

	if (!bmp->LockBits(&rcLock, Gdiplus::ImageLockModeRead, format, &bmpData) == Gdiplus::Ok)
		return cv::Mat();

	cv::Mat mat = cv::Mat(hgt, wd, CV_8UC3, static_cast<unsigned char*>(bmpData.Scan0), bmpData.Stride).clone();

	bmp->UnlockBits(&bmpData);
	return mat;
}

cv::Mat GdiplusImageToOpenCvMat(Gdiplus::Image* img)
{
	auto bmp = std::unique_ptr<Gdiplus::Bitmap>(GdiplusImageToBitmap(img));
	return (bmp != nullptr) ? GdiPlusBitmapToOpenCvMat(bmp.get()) : cv::Mat();
}

The short answer is that in the hexagonal case the best analog of Conway’s Game of Life — in my opinion as someone who has been a CA hobbyist for 30 years or so — is an original creation which I will describe for the first time in detail in this blog post. (Note: The links to cellular automata in this post go to a custom web-based player which may not run correctly on mobile devices.)

Regarding a hexagonal Game of Life, a key thing to understand is that Conway Life isn’t just the rule; it is the rule (Life 2333) running on a square grid with an 8-cell neighborhood i.e. the neighborhood of four squares directly adjacent to a square plus those diagonally adjacent. You can, of course, apply the same rule on a hexagonal grid using the natural six hexagon neighborhood but what you will get won’t look anything like Life. It will look like this: Life 2333 on a hexagon grid.

So if the above does not qualify as a hexagonal Game of Life then what would? Well, at the very least we need a glider. Carter Bays, a professor at the University of South Carolina, presented a Conway-like rule that admits a glider on the hexagonal lattice in 2005, Life 3,5/2 in his notation. However, by Bays’ own admission “this rule is not as rich as Conway’s Life” and indeed when we run Bays’ Hex Life on random input we do not see gliders: Life 3,5/2 on random input. The problem is that its glider is too big to occur randomly. Its glider is in a sense artificial. Part of the beauty of Conway Life is that gliders are frequently spontaneously generated. The other thing you’ll notice about Bays’ Life 3,5/2 is that it descends into still-lifes and oscillators too quickly. Conway Life is dynamic. It sprawls and grows, descends into bounded chaos, before finally decaying into still-lifes and oscillators.

To summarize, we want a hexagonal cellular automaton that

1. Has a glider that is frequently generated by random initial input.
2. Frequently exhibits bounded growth from random initial input.

It is my contention that there is no simple totalistic rule over two states and the hexagonal grid using the natural 6-cell neighborhood that exhibits both 1. and 2.

In order to achieve what we want, we need to drop one of the constraints. Using my cellular automata breeding application Lifelike, I explored dropping the constraint that the rule must be a simple totalistic rule over two states. What I have come up with is a cellular automaton that uses a simple totalistic rule over three states,  states 0 to 2. One way to view this move is to view the live cells in Conway Life as counters — beans, pennies, whatever — and imagine dropping the constraint that a cell can only contain one counter. Anyway, my rule is as follows:

• Take the sum S of the 6-cell neighborhood where death=0, yellow=1, and red=2.
  • If the cell is currently dead it comes to life as yellow if S is exactly 4.
  • If the cell is yellow it goes to red if S is 1 to 4 inclusive or is exactly 6.
  • If the cell is red it stays red if S is 1 or 2 and goes to yellow if S is 4.
  • Otherwise it is dead in the next generation.

The above, which I call “Joe Life”, is as rich as Conway life: click this link to see it run. It exhibits bounded growth with about the same burn rate as Conway Life and features two gliders, the little fish and the big fish, which are frequently generated spontaneously.

The little fish

The big fish

I concede Conway’s rule is more elegant than Joe Life’s rule but if one thinks about it, Conway Life’s neighborhood is less natural than the six hexagon neighborhood in that it is kind of weird to include the the diagonally adjacent squares. So in my opinion, elegance that Joe Life lacks in its state table it makes up for in its neighborhood.

A friend of mine, Jack Kutilek, wrote a web-based player for the CA rules format that I used in Lifelike — basically simple JSON files that describe a state table and the meta-information you need to execute the CA the state table encodes. It uses WebGL and a fragment shader and as such is very fast.

I made a little web page interface here, some cellular automata, which just embeds Jack’s fragment shader thing via an iframe. The naked player is here, running Joe Life by default. You can run it on different rules by providing it  URI-encoded Lifelike JSON as an HTTP query string.

A couple of weeks ago I got interested in this project but wanted full control of the code, wanted to know exactly what it is doing, wanted a bunch of features like the ability to import and export CA rules, and wanted to have the process not be seeded by cellular automata already featuring gliders (which the web app seems to be). To this end I have pushed a project to github that does something similar but, I hope, more transparently.

I’m calling it Lifelike Cellular Automata Breeder. It is a (C# WinForms) application in which given some settings a user can artificially select and breed cellular automata; i.e., it performs a genetic algorithm in which the user manually provides the fitness criteria interactively.

I decided I wanted to only allow a reproduction step in which I scramble together state tables in various ways, guessing that using “DNA” more complex than commensurate 2D tables of numbers wouldn’t work well for a genetic process in this case. I characterize CA rules as applying only relative to a given number of states and a given, what I call, “cell structure” and “neighborhood function”. Cell structure just means a lattice type and neighborhood –e.g. square, with four neighbors; square, with eight neighbors; hex, with six, etc. “Neighborhood function” is an arbitrary function that given the  states of the n cells in some cell’s neighborhood returns an integer from 0 to r where r is dependent on the neighborhood function and possibly number of states. For example, Conway Life uses the neighborhood function I refer to as “alive cell count”, and for an n-cell neighborhood, r equals n because the greatest number of alive cells that can surround a cell is just the size of the neighborhood. If the user has selected s total states, the state tables will be s by r.

Lifelike works as follows

1. The user selects a number of states, cellular structure, and neighborhood function and kicks off the genetic process.
2. Lifelike sets the current generation to nil, where by “generation” we just mean a set of cellular automata that have been tagged with fitness values.
3. While the user has not clicked the “go to the next generation” button,
   • If the current generation is nil, Lifelike randomly generates a cellular automata, CA, from scratch by making an s by r state table filled with random numbers from 0 to s. (The random states are generated via a discrete distribution controllable by the user). If the generation is not nil, Lifelike selects a reproduction method requiring k parents, selects k parents from the current generation such that this selection is weighted on the fitness of the automata, generates CA using the reproduction method and parents, and then possibly selects a random mutation function and mutates CA, selecting the mutation function via a discrete distribution controllable by the user and applying it with a “temperature” controlled by the user.
   • Lifelike presents CA in a window.
   • The user either skips CA in which case it no longer plays a role in the algorithm or applies a fitness value to it and adds it to the next generation.
4. When the user decides to go to next generation, the selections the user just made become the new parent generation and processing continues.

Here’s a video of me playing with an early version of the application.

Results, Musings, etc.

Briefly, Lifelike works.

You can produce interesting cellular automata with it and there is a weird feeling that is hard to describe when you first see a tiny glider wiggle across the screen; however, the way it works is somehow more mundane than I thought that it would be. I wonder if all genetic algorithms are like this. Most of what it produces is garbage. When you see something that isn’t garbage you can select for it. However Lifelike doesn’t do magic. It doesn’t magically find phenomena you want just because you have a fancy framework implemented to find such phenomena. For example, it is my belief at this point that there is no simple hex grid analog of Conway Life using alive cell count, a six cell neighborhood, and 2-states. I think you could probably prove the non-existence of gliders in this configuration but it would be a boring proof by exhaustion and running Lifelike on that configuration is boring as well. Just because you, the user, are a step in a “genetic algorithm” doesn’t somehow make it interesting.

The simple hex neighborhood negative result led me to ask the following question: What is the smallest change you can make to hex/6-cell/alive cell count/2 states to allow Conway Life-like behavior? If you google “Hex Life” you will find that it is well-known that interesting things can happen if you go to a 12-cell  Star of David shaped neighborhood, but this seems inelegant to me because the simple hex neighborhood is so nice. The question then is are gliders possible in the simple hex lattice and neighborhood if we add one state and modify “alive cell count” in a trivial way? The answer to this question turns out be yes. There are beautiful rules that live in the hex lattice with the natural neighborhood if we have 0 = death, 1 = Alive-A, 2 = Alive-B and instead of the simple alive cell count we use its natural analog when states can have values greater than one: sum of states.

Below is a such a rule set and is probably the best thing that has come out of my work with Lifelike as far as I am concerned (So I am naming it Joe Life, assuming that it is unknown in the literature).

The above has a nice quality that Conway Life also has that I call “burn”. This a qualitative thing that is really hard to define but it is what I look for now when I play with Lifelike: burn + gliders = a Life-like cellular automaton. “Burn” is the propensity of a cellular automata configuration to descend into segregated regions of chaos that churn for awhile before ultimately decaying into gliders, oscillators, and still lifes. Some CAs burn faster than others; the above has a nice slow burn. CAs that exhibit steady controlled burn turn out to be rare. Most CAs either die or devolve instantly into various flavors of unbounded chaos.

However there does turn out to be another quality that is not death or unbounded chaos that is sort of like the opposite of burn. See for example

(The above is hex 6-cell, four states, and using a neighborhood function I call “state-based binary”)  which I have been calling “Armada” and generally have been referring to these kind of CAs as being armada-like. Armada-like cellular automata quickly decay completely into only weakly interacting gliders. For example, one from the literature that I would characterize as armada-like is Brian’s Brain. Armada-like rules turn out to be more common than life-like rules. They’re impressive when you first start finding them but they are ultimately less interesting, to me at least. The best thing about armada-like rules is that they indicate that life-like rules are probably “nearby” in the space you are exploring in the genetic process.

Also they can breed weird hybrids that defy classification, such as the following which are all burn with large blob-like gliders and seem sometimes to live around the boundary between armada-like rules and life-like rules.

Magic Carpets (square, 4-cell/ 4 states / sum of states)

or Ink Blots (hex, 6-cell/ 3 states / “0-low-med-high”)

My other major result is that life-like rules exist in the square 4-cell neighborhood if we allow an extra state and use the simple sum of states as the neighborhood function, but they can be boring looking so instead here is an armada-like square 4-cell CA that is on the edge of being lifelike:

The above uses the neighborhood function I call “2-state count” which enumerates all possible combinations of c1, number of neighboring cells in state 1 and c2, number of neighboring cells in state 2 or above, in an n-cell neighborhood i.e. c1 + c2 ≤ n.

You can weave an icosahedron from 10 strips of paper — or anyway you can weave a construction that has icosahedral geometry; it is actually more of a snub icosahedron.

From about 11″ strips of printer paper folded lengthwise to have two layers with a little overlap to lock into rings, as pictured below, locks well and is rigid:

Construction follows the pattern implied by the following:

Below is a modular origami construction known as the 120-unit sonobe ball. It has the geometry of an icosidodecahedron but with the 12 pentagonal faces of the icosidodecahedron split into five of the little equilateral pyramids that sonobe units demand.

The above was made from 120 Post-it notes, three sets of forty in each color. I haven’t seen anyone post a picture of one of these 3-colored like this, i.e. in three colors such that no two units of the same color are adjacent to each other, so I thought it was of some interest…

A schematic of the coloring I found is below. I have no idea if it is unique:

I’m the kind of reader of Neal Stephenson who likes Snow Crash, Diamond Age, and Anathem but has never succeeded in plodding through Cryptonomicon (I tried twice!) and who has never even tried to read the Baroque Cycle, not being a masochist. I’m a reluctant fan. I wish he would  realize that just because the research he has been doing on the historical use of mercury in silver mining (say) is interesting to him, it isn’t necessarily interesting to everyone else.

His books can often be characterized this way, by the odd topics he shoehorns into them: Snow Crash contained a lot about the Sumerian language and various bits of esoterica from computer programming, Cryptonomicon had Enigma machines and the Theory of Computation, Anathem was ultimately about the Many Worlds interpretation of quantum mechanics mixed with Roger Penrose’s anti-computational theory of mind. The success of these books depends on whether or not the bits of trivia he is throwing at you begin to stick, whether they add up to anything — Does their interest become clear and does their inclusion attain a kind of inevitability?

Seveneves is concerned with orbital mechanics, the usefulness of swarming behavior in robotics, and the physics of tethers and chains — the latter of which Stephenson is determined to convince us is not only fascinating but shockingly under-studied. Somewhat surprisingly, this time around, it all comes together. The book works because ultimately it is about the nature of technology. The fun science facts are just the lens through which he examines that nature. The book is about the extent to which the human race needs technology to survive despite the fact that the technology itself is mostly mundane: ultimately the difference between a space tether and a loop of string is a few millenia plus the fact that former is a lot bigger.

The book is two stories in one. The first part is a defiantly retro “space cowboy” yarn, a disaster story in which a cast of characters has two years to set up a sustainable community in near earth orbit because the earth is about to become uninhabitable for five thousand years. It is a throwback to days when NASA was culturally relevant and space-walks and rocket launches were cool. In the Acknowledgements Stephenson mentions a staff member at Planetary Resources, an actual asteroid mining company he visited for research purposes, being pleasantly surprised that “someone was producing science fiction where the asteroid mining company was, for once, the good guys” and this pretty much sums up this portion of the book. It reminds me of hard science fiction from the 1970s, something like Arthur C. Clarke’s Rendezvous with Rama but with Neal Stephenson’s typical geek-chic type characters.

The second part deals with the descendants of the first part and the first part casts a shadow on this story in interesting ways. There is a kind of inevitability to it. For example, at one point it is important for a particular character to know Morse code and given the preceding six hundred pages it is endearingly believable that he would. It makes perfect sense that a Morse code knowing tradition would have survived because of plot details in the space cowboy portion of the book and it leads one to think about technology as cultural transmission, as useful memes. Much of Seveneves is like this. It is a book in which over the course of a eight hundred pages a melodrama set on the International Space Station leads naturally to a genetically engineered neanderthal wielding a whip made out of tiny robots, which is what science fiction is all about.

Below is gameplay video of a prototype of the puzzle game “Draak” that I am working on. Draak will be an iPad-only iOS game; the prototype below is written to Windows in C# on top of MonoGame.

The initial idea for this game was about the art. I thought it would be cool to make a game that looks like an animated version of one of M. C. Escher’s tessellations, particularly one that involves a similarity symmetry. So I lifted mechanics from a 1990s era game that I always felt was a good game trapped in the typical square lattice for which it was particularly ill-suited. The old game’s Euclidean lay-out wasted a lot of space.

A lot of the work in the above may not be immediately apparent (which is good I think). Specifically, there are actually two shapes of butterflies in the above not just one. They look like this:

and are arrayed as in a checkerboard with a 90 degree rotation applied to cells of opposite parity — the geometry might be clearer here in which I create these butterfly tiles unadorned with my tool EscherDraw (before I had beziers working in EscherDraw) Thus any time a butterfly moves to a cell with opposite parity I need to not only handle the scale and rotation tweening imposed by the spiral lattice, I also have to simulatenously do the 90 degree rotation imposed by the butterfly tiles and simultaneously morph the actual sprites above from one to the other.

I created the above sprites as vector art using tile outlines exported from Escher draw as SVG and wrote a Python script that morphs SVG as the publically available software for morphing vector graphics is surprisingly non-existent. I am going to post about my vector morphing code here soon, I just need to clean it up for use by people other than me.

Creating the draak logo in EscherDraw with bezier curves:

I just added full cubic bezier curve support to the vector tessellation creation tool I am developing (see here). I’ll include some video at a later date (Update: here) showing this feature in action but first I want to document and make publicly available some of the code involved: namely, a bezier class in C# that supports two functions that are not the easiest things in the world to implement. Both of these functions come up as soon as you try to do anything non-trivial with beziers and the literature online can be unhelpful.

My code is here: link.

The first function is splitting a cubic bezier into two curves at a given t parameter such that the concatenation of the two curves is congruent to the original. This turns out to be easy; however, I could find no actual code on offer, rather just descriptions of code. There is good coverage elsewhere on the theory so I won’t go into it deeply here, but, briefly, splitting a bezier is easy because how one does so follows naturally from the sort of recursive definition of the bezier that is implied by De Casteljau’s algorithm. Anyway here is my code:

...
 public Tuple<Bezier,Bezier> Split(double t) {
     PointD E = Interpolate(t, CtrlPoint1, CtrlPoint2);
     PointD F = Interpolate(t, CtrlPoint2, CtrlPoint3);
     PointD G = Interpolate(t, CtrlPoint3, CtrlPoint4);
     PointD H = Interpolate(t, E, F);
     PointD J = Interpolate(t, F, G);
     PointD K = Interpolate(t, H, J);

     return new Tuple<Bezier, Bezier>(
         new Bezier(CtrlPoint1, E, H, K),
         new Bezier(K, J, G, CtrlPoint4)
     );
 }

 static private PointD Interpolate(double t, PointD pt1, PointD pt2) {
     double x = (1.0 - t) * pt1.X + t * pt2.X;
     double y = (1.0 - t) * pt1.Y + t * pt2.Y;
     return new PointD(x, y);
 }
....

The other function is considerably harder: finding the nearest point on a bezier curve to a given point. If you scan the literature (i.e. perform the google search) you will find that what most people end up using is some code from an article in the book Graphics Gems I, “Solving the Nearest-Point-On-Curve Problem” by Philip J. Schneider. I, however, have a couple of problems with this code (1) I don’t really understand it, (2) it is long and I need to port it to C#, and (3) someone on the internet is claiming that it is incorrect but given (1) I can’t really evaluate the claim; see the comment that begins “There seem to be errors in the ControlPolygonFlatEnough function […]” in the C code here.

But more relevantly I just don’t get this code on a fundamental level. The nearest point on a bezier problem is difficult to solve mathematically but it is easy to formulate mathematically. I think that the Graphics Gems code obscures the straightforward aspect of how this problem can be formulated. Further, I don’t have a problem using code that I don’t understand; however, if I am going to port an opaque block of code from one language to another I’d like to keep that portion of code to a minimum and I can do this easily by structuring my code around the formulation of the problem that I understand.

The formulation I am talking about is as follows, if B(t) is a particular cubic bezier function, B’(t) is its derivative, and P is an arbitrary point then [B(t) – P] ⋅ B’(t) = 0 when t is the parameter of the closest point on B to P. [B(t) – P] is a vector that points from P towards some point on the curve and B’(t) is a vector that is tangent to the curve at this point, if the distance is going to be a minimum then these two vectors will be perpendicular and thus their dot product will be zero. Since B(t) is cubic, its derivative will be quadratic and thus when we take the scalar product of [B(t) – P] and B’(t) we will end up with a quintic equation of one variable. This leads naturally to a nearest point function that looks like the following on the top level:

 public Tuple<PointD, double> FindNearestPoint(PointD pt)
 {
     var polyQuintic = GetNearestBezierPtQuintic(
         _points[0].X, _points[0].Y,
         _points[1].X, _points[1].Y,
         _points[2].X, _points[2].Y,
         _points[3].X, _points[3].Y,
         pt.X, pt.Y
     );
     List<Complex> roots = FindAllRoots(polyQuintic);

     // Filter out roots with nonzero imaginary parts and roots
     // with real parts that are not between 0 and 1.
     List<double> candidates = roots.FindAll(
         root => root.Real > 0 && root.Real <= 1.0 && Math.Abs(root.Imaginary) < ROOT_EPS
     ).Select(
         root => root.Real
     ).ToList();

     // add t=0 and t=1 ... the edge cases.
     candidates.Add(0.0);
     candidates.Add(1.0);

     // find the candidate that yields the closest point on the bezier to the given point.
     double t = double.NaN;
     PointD output = new PointD(double.NaN,double.NaN);
     double minDistance = double.MaxValue;
     foreach (double candidate in candidates)
     {
         var ptAtCandidate = GetPoint(candidate);
         double distance = DistSqu(ptAtCandidate, pt);
         if (distance < minDistance)
         {
             minDistance = distance;
             t = candidate;
             output = ptAtCandidate;
         }
     }

     return new Tuple<PointD, double>(output, t);
 }

which reduces the problem to implementing GetNearestBezierPtQuintic() and a numeric root finding function for polynomials, given that quintic equations cannot be solved via a closed-form formula like the quadratic equation.

GetNearestBezierPtQuintic() –which returns the coefficients of [B(t) – P] ⋅ B’(t) when fully expanded given P and the control points of B — turns out to be a little two much algebra to work out comfortably by hand, so I ran the open source symbolic mathematics application Sage on this Sage script yielding the following when translated into C#: (For the first time on the internet! … as far as I can tell)

static private List<Complex> GetNearestBezierPtQuintic(double x_0, double y_0, double x_1, double y_1,
            double x_2, double y_2, double x_3, double y_3, double x, double y)
{
    double t5 = 3 * x_0 * x_0 - 18 * x_0 * x_1 + 27 * x_1 * x_1 + 18 * x_0 * x_2 - 54 * x_1 * x_2 + 27 * x_2 * x_2 -
        6 * x_0 * x_3 + 18 * x_1 * x_3 - 18 * x_2 * x_3 + 3 * x_3 * x_3 + 3 * y_0 * y_0 - 18 * y_0 * y_1 +
        27 * y_1 * y_1 + 18 * y_0 * y_2 - 54 * y_1 * y_2 + 27 * y_2 * y_2 - 6 * y_0 * y_3 + 18 * y_1 * y_3 -
        18 * y_2 * y_3 + 3 * y_3 * y_3;
    double t4 = -15 * x_0 * x_0 + 75 * x_0 * x_1 - 90 * x_1 * x_1 - 60 * x_0 * x_2 + 135 * x_1 * x_2 -
        45 * x_2 * x_2 + 15 * x_0 * x_3 - 30 * x_1 * x_3 + 15 * x_2 * x_3 - 15 * y_0 * y_0 + 75 * y_0 * y_1 -
        90 * y_1 * y_1 - 60 * y_0 * y_2 + 135 * y_1 * y_2 - 45 * y_2 * y_2 + 15 * y_0 * y_3 - 30 * y_1 * y_3 +
        15 * y_2 * y_3;
    double t3 = 30 * x_0 * x_0 - 120 * x_0 * x_1 + 108 * x_1 * x_1 + 72 * x_0 * x_2 - 108 * x_1 * x_2 +
        18 * x_2 * x_2 - 12 * x_0 * x_3 + 12 * x_1 * x_3 + 30 * y_0 * y_0 - 120 * y_0 * y_1 +
        108 * y_1 * y_1 + 72 * y_0 * y_2 - 108 * y_1 * y_2 + 18 * y_2 * y_2 - 12 * y_0 * y_3 + 12 * y_1 * y_3;
    double t2 = 3 * x * x_0 - 30 * x_0 * x_0 - 9 * x * x_1 + 90 * x_0 * x_1 - 54 * x_1 * x_1 + 9 * x * x_2 -
        36 * x_0 * x_2 + 27 * x_1 * x_2 - 3 * x * x_3 + 3 * x_0 * x_3 + 3 * y * y_0 - 30 * y_0 * y_0 - 9 * y * y_1 +
        90 * y_0 * y_1 - 54 * y_1 * y_1 + 9 * y * y_2 - 36 * y_0 * y_2 + 27 * y_1 * y_2 - 3 * y * y_3 + 3 * y_0 * y_3;
    double t1 = -6 * x * x_0 + 15 * x_0 * x_0 + 12 * x * x_1 - 30 * x_0 * x_1 + 9 * x_1 * x_1 - 6 * x * x_2 +
        6 * x_0 * x_2 - 6 * y * y_0 + 15 * y_0 * y_0 + 12 * y * y_1 - 30 * y_0 * y_1 + 9 * y_1 * y_1 -
        6 * y * y_2 + 6 * y_0 * y_2;
    double t0 = 3 * x * x_0 - 3 * x_0 * x_0 - 3 * x * x_1 + 3 * x_0 * x_1 + 3 * y * y_0 - 3 * y_0 * y_0 -
        3 * y * y_1 + 3 * y_0 * y_1;

    return new List<Complex> { (Complex)t0/t5, (Complex)t1/t5, (Complex)t2/t5, (Complex)t3/t5, (Complex)t4/t5, (Complex)1.0 };
}

and decided to use Laguerre’s Method to solve the equation.

I chose Laguerre’s Method because it is optimized for polynomials, is less flaky than Newton-Raphson, is relatively concise, and is easy to port from C++ given .Net’s System.Numerics’s implementation of complex numbers. I ported the implementation of Laguerre’s Method found in Data Structures and Algorithms in C++ by Goodrich et. al. (here). This solver works well enough for my purposes; however, I think the optimal thing to do to get the optimal nearest-point-on-a-bezier algorithm would be to implement one of a class of newish algorithms that are specifically targeted at solving the quintic, e.g. “Solving the Quintic by Iteration” by Doyle and McMullen. No implementations of these algorithms seem to be publicly available, however, and it would be a research project that I don’t have time for right now to translate that paper from academese into code.