Somebody with more experience, please confirm on hackernews (https://news.ycombinator.com/item?id=

]]>Somebody with more experience, please confirm on hackernews (https://news.ycombinator.com/item?id=22612139)

Now that Covid-19 is out, we need to figure out how to detect all infected people so that uninfected people can return to work/school/life. Further, we need continually test everyone who is assumed uninfected weekly to prevent future outbreaks. This test cannot be an antibody test, as an infected person may be contagious before the immune system has created the antibodies. The only current way to detect the virus in an asymptomatic individual is RT-PCR of the viral RNA.

Let's take Denmark, where I now live. That's 5 million people and therefore tests per week. If the current cost of a viral RNA test is $200, that's $1 billion per week or $52 billion per year. These tests are usually done on 96 well plates with 2 control wells. That leaves 94 tests per plate (which can be thought of as batch size for tests).

However, the assumption is that currently most people are not infected. This is especially true if our goal is to weekly test the people returning to daily life who are assumed to be uninfected. Therefore, most of the wells are "wasted" on uninfected individuals.

Here's how to test 4731 people per 96 well plate: take two samples from each individual, and mix the samples in a binary encoding scheme that assign 2 bits of a 94-bit number to each person.

To illustrate, consider first this how to test 6 people using 4 sample wells (forget the controls for now). If each sample well is either red for positive or blue for negative, there are 6 ways to color the wells such that exactly 2 sample wells are positive.

The stategy is to assign each person to one of the above colorings, and then in each well mix the samples from testees such that if one (and only one) person is positive, the resulting coloring will match the assigned person. To do this, place each person's two samples in the wells as shown below.

The encoding is equivelent to enumerating the 6 binary numbers between 0 and 15 with two set bits, and then assigning one person to each of those numbers. With 4 wells, we can test 4 choose 2 people = 6. With 94 wells, we can test 94 choose 2 people = 4731.

The PCR step (including preparation) is the bottle-neck in the testing process. For the 94 well method to work, each well would need to combine samples from 93 people, such that 8742 samples are combined to make a single 96 well plate (4731 people times 2 samples each = 8732 samples). The good news is that the test samples can be combined into the 94 groups (one for each well) right at the beginning of the process (literally, each group of 93 swabs could be put together into a large beaker and mixed). All other downstream operations then only need to be done once per well.

We need a mechanical device which can automatically sort the samples into well groups. Correctly combining 4731 samples is not something that can be done by hand. However, compared to the PCR machines and magical reagants, this sorting machine should be relatively simple to build. The reward is a 46.5 times increase in testing throughput.

Yes, there is also the possibility that more than one of the 4731 people is infected, and in this case you would get 3 or more positive wells. Still, this method can be used to determine a subset of samples which need individual retesting. However, if most samples are assumed to be uninfected, this is not a huge issue.

]]>Things finally slowed down as I ran out of administrative tasks to chase down. I powered up the robot again and got it talking on the wifi. The video issues I was experiencing seem to

]]>Things finally slowed down as I ran out of administrative tasks to chase down. I powered up the robot again and got it talking on the wifi. The video issues I was experiencing seem to have disappeared. Previously I had a suspected issue with dropped USB video packets, which resulted in image shearing. An electrical engineer I talked to was suspicious of the power input adjacent to the USB connectors causing interference, and my best guess now is that my new bedroom power socket has less noise than the previous power supply.

I also developed some maths for a robust calibration model for the robot's cameras. Most models assume uniform radial optical distortion, which is a reasonable for camera lenses. My optical setup includes an air-water interface separated by acrylic, and measurements from an underwater sequence I previously collected suggest that a simple radial model is inadequate. Week 3 was otherwise quite, except for pulling my quad playing football.

On Wednesday of week 4 a pitch training session was held for the September investor summit. I'm extremely grateful to Attila Hadnagy who has been doing my graphic design work. The Fishi pitch received very positive feedback, but I also received some suggestions which I think reflect the local mentality. In Silicon Valley if you don't talk big numbers no one is interested. Here they want you to get to the point, skip the big numbers and jump straight to ship inspections.

The smaller thinking also introduces what I'm calling the Danish discount. Since the area started focusing on robotics, there has been one $250 million exit ([Universal Robotics](https://www.universal-robots.com/)), and another ([MiR](https://www.mobile-industrial-robots.com/en/)) in the works. There simply haven't been enough outlier exits to support Silicon Valley style valuations for seed stage companies.

I still don't have the complete picture, but it seems at least some companies are giving up 30% to raise their first $100k. That hit encourages bootstrapping, which in turn encourages a bit of fake it before you make it. I talked to one company of three, with only a single technical person. They have one or two customers for whom they are basically piloting drones manually. The selling point of the company is of course that one day the drones will be autonomous, but I simply don't see how they can raise the funds to hire the people needed to make that a reality. They are operating in a crowded space where several SV companies have raised well over $10 million a piece, and the European market will only protect them for so long.

I was full aware that raising money in Europe would be no where near as easy as SV. To be critical of myself, I should have put a number on the Danish discount before I came. All things being equal, my guess is that an early stage Danish company can expect a valuation around 20% of what one would expect in SV.

I'm not exactly comfortable accepting that, but the good news for a technical founder is that bootstrapping really is an option. To help make up the difference, there are government grants available. Other funding sources, e.g. from the Danish Maritime Fund, provide high risk business loans to scale new technologies aligned with industry. Explosive early stage growth, though, seems tricky.

The other silver lining to being here is the proximity to [Svendborg Maritime Academy](http://www.simac.dk/en/about-simac/the-simac-organisation/) which effectively operates a teaching and research port, similar in function to a university hospital. I've been told that all major Danish shipping companies have a relationship with the academy, and the port could serve as a testing ground for the MVP. Verifying this claim is high on my priority list.

The investor summit is September 10th. The ideal would be to meet one or two strategic Danish angel investors. If the people and terms are right, I'll consider investment. My current mindset is to build the MVP myself, and in the mean time find the right angel investors for a later seed round.

While the summit will be the best chance I've had so far to meet the local players, it is also proving a distraction. For these last two weeks, I've been coding on and off on my camera calibration method. Last Friday I finally sat down in a cafe, killed my notifications and pretended like I had no other commitments. It was somewhat reassuring to have such a productive day.

In other news, I applied for the YC startup school and was accepted into the Advisor Track. Unofficial numbers seem to say 25% of the applicants were accepted, which is a nice little confidence booster.

Next week I have my screening interview for the local, city funded incubator. I'm also beginning the process of incorporating Fishi Robotics.

]]>This week has been mostly administrative. I finally have a lease, which was prerequisite for banking, health coverage, and obtaining idenity cards. I also had my first contact with the RoboHub and received some positive signals.

On Tuesday I went down to Cortex Park where Odense Robotics is based and

]]>This week has been mostly administrative. I finally have a lease, which was prerequisite for banking, health coverage, and obtaining idenity cards. I also had my first contact with the RoboHub and received some positive signals.

On Tuesday I went down to Cortex Park where Odense Robotics is based and received a thorough orientation. The region of 250,000 people has some pretty impressive statistics: 120 robotics companies, 500 supporting companies, 8 figure exits. I still need to apply to the incubators, but the next application window is September and I received some pretty positive hints. First, I learned that there are actually three separate incubators, one for robotics, one for drones, and one for maritime. Since I could technically fit into any of those, I was told the intake quotas could be juggled in my favour. Secondly, I've been slipped into the Odense Investor Summit which will put me in front of local and European investors.

We also discussed available grants. As an open trade block, the Eurozone has rules preventing governments from subsidising local industry. However, there is a de minimis rule allowing companies to receive up to €200,000 in state support. Not coincidentally, there is an InnoBooster fund with simplified rules for grants up to this ammount. The basic rule is that they will fund 33% of a project, but the accounting rules let you bill an unpaid founder at (conveniently) €200,000 per year as payment in kind. For a two cofounder startup, this is basically free money. In my case, I will need to get a little creative, likely by billing myself out over 18 months and raising €100,000 to match the rest of the grant. The incubator also provides advisors and free office space in a package worth around €80,000.

Back to the administrative side, on Monday I signed a 6 month sub-lease after a little drama. It happened that both rooms in the appartment were changing over at the same time, one permanently and one temperarily during a doctural exchange program. Danish rental law provides very strong protections for long term tenants, but there is no provision for wear and tear and any "deposit" is essentialy forfeit. For short term tenants, rather than losing the deposit it is common to take over an existing lease, but there are horror stories where the new tenant was exposed to the liabilites of the old. A new tenant was found for the sub-leased room and I was to take over the other room on the actual lease, but the sub-leasee bailed. Both the exchange student and I were under time pressure, so we went ahead with a sub-lease on his room. In hindsight this has worked out well, as I have a fully furnished room and wont lose my deposit.

The City of Odense has a special International Citizens Service center, which like the rest of Denmark seems to be on holiday for the summer. Luckly the usual Citizens Service was able to process me, and I'm now set up with health insurance and NemID. NemID is a national digital identity system which is used for online interactions with the government and banks, with a 2FA system using either codes sent via snail-mail or a phone app. While not yet mandatory, it is also possible to link a NemID to a public certificate for cryptographically signing digital documents. For banking, SydBank was recommended to me, but even when applying in person the paperwork is sent via DigitalPost, which is again linked to NemID. I won't be able to complete the process until my first codes arrive in the mail. Getting a phone number turned out to be tricky, seemingly because it is a conduit for credit card fraud. Various carriers asked for a European passport or 3 months of Danish pay stubs. Telia finally gave me a phone number, but it is also the most expensive option.

This week I'm hoping to crack open the robot and start writing code again. There are also preparations for the Investor Summit and the RoboHub is helping me with another iteration for the pitch deck.

]]>After a week in Odense I'm feeling pretty optimistic about my move. I have found a place to live, been put in touch with the right people, toured the local airport, and above all found everyone to be extremely helpful.

Odense has four universities/colleges serving a city

]]>After a week in Odense I'm feeling pretty optimistic about my move. I have found a place to live, been put in touch with the right people, toured the local airport, and above all found everyone to be extremely helpful.

Odense has four universities/colleges serving a city of 170,000 people, and therefore students make up a large percentage of the population. Online room listings move quickly, but a few days before arriving I heard back from a German Ph.D. studying bioinformatics who will become my new roommate. However, I've since learned from locals that the facebook group Lejligheder til salg og leje i Odense og omegn is where everyone posts their rooms, so skip the listings sites.

For now I'm staying at an airbnb with a cat for company while my host is in Sweden. Two of her friends have been dropping by to look after the cat, and they have both been helpful. The apartment is about as central as you can get and backs onto a large public square. Unfortunately, since Wednesday there has been an ongoing music event with a stage right outside of my window. The music seems to wrap up at about midnight, so I'm currently missing a few hours of sleep.

A social hack I've borrowed from San Francisco is to go find the local Aussie Rules Footy team, though remarkably in a squad of about 20 I'm one of 2.5 Aussies. Somehow the locals have even erected propper goal posts! Training is conducted in English for the sake of the few of us who wouldn't understand otherwise, but after practice beers are in dansk. Another foreign student here told me he won't go out unless the foreigners outnumber the Danes, otherwise they won't speak English. It is a bit swim or sink, but at least there's an invitation here to use some Danish while being laughed at. Most other attempts would be met with a reply in perfect English.

Once I had decided to move to Denmark, everyone I talked to pushed me towards Odense. Since at least 2013, the country and city have been investing heavily to turn Odense into a robotics hub. University Syd Denmark offers degrees in Mechatronics and Robot Systems, with a focus on drones. It has acquired a large hangar with offices at the local airport for both teaching purposes and use by affiliated companies. There is a one-of-its-kind UAV testing zone where drones can be flown autonomously outside line of sight view. It covers over 1000 square kilometers, extends out over the ocean, and is used by everyone from startups to Boeing and the military. There is also a city funded incubator, Odense Robohub for robotics startups that I will visit next week and hopefully become a part of.

A week earlier I started sending out emails to the listed contact addresses of nice to know people, including some at the Robohub. However the emails went unanswered, and on day one I truly had the sensation of "what have I done", as I found myself in a new city without a clear idea of what needed to be done first. In a bout of minor desperation, I shotgunned the founders of the other companies in the hub and applied to another program meant for non-robotics startups. One company, QuadSAT got back to me and the next day a founder gave me an airport tour. My application also received a reply, which referred me back to the Robohub but this time with an introduction to the correct person. I suppose this is one advantage of operating in a smaller community where everyone knows everyone.

QuadSAT has an enviably solid business case, drones that mimic satellites for calibrating antennae. Satellite connectivity is now so important to ships that contracts are suspended if connectivity is lost. Systems are currently tested by taking a ship out of port and motoring around to see what happens. It almost seems too obvious, but having a drone fly the pattern instead is orders of magnitude cheaper and convenient.

A final tidbit which I thought was just way too cool. Before leaving San Francisco I found that the holes on one of my 3D printed parts were incorrectly placed, but I didn't have time to print a corrected piece. As part of my shotgunning, I found FabLab which I thought was a local maker space. It's actually a college facility with a fairly impressive collection of laser cutters, CNC routers, milling machines and a small army of 3D printers. Because it's funded by public money, it sits in a catch 22 where it should be open to the public, but isn't allowed to compete with local businesses. Therefore it has one open day a week where anyone can just walk in and use the machines. They also didn't seem too concerned about who was paying for the filament so long as your part was within reason. I saw another part being printed that was much larger than mine. So yay, a public prototyping workshop! Go Denmark!

Next week I'll hopefully have a lease agreement which will let me start on administrative matters. I'll also be going down to talk with the Robohub to start my application there. As for the weekend, the footy team has invited me on a trip out to Malmö for a game against the Swedes. The bus trip back is portended to be a bit of a booze cruise.

]]>I want to share a plot which took me by surprise. In computer vision, a useful tool for estimating motion is image alignment. Commonly, one may wish to estimate the affine warp matrix which minimizes a measure between pixel values of two images. Below are two aligned images from the

]]>I want to share a plot which took me by surprise. In computer vision, a useful tool for estimating motion is image alignment. Commonly, one may wish to estimate the affine warp matrix which minimizes a measure between pixel values of two images. Below are two aligned images from the KITTI data set, which are demonstrative of forward motion.

An alternative, which I believe was pioneered by Semi-Direct Visual Odometry, is to align two images by estimating the six parameter transform between the cameras. Usually these are the parameters of the $(\mathfrak{se}(3)$), but could also be three Euler angles plus an $(\langle x, y, z\rangle$) Euclidean transform. Because this method essentially reprojects the first image into the second, the depth field also needs to be known a priori.

The mathematics involved is surprisingly simple. For a feature based system, several points of known depth will be chosen from the first image $(\mathrm{I}_1$). Around each point, a small rectangular patch $(\mathrm{P}$) of pixels will be extracted and then compared to the second image $(\mathrm{I}_2$) at an offset of $(\langle x, y\rangle$). The goal is for each patch to find an offset which minimizes the below sum.

$$

\begin{equation}

\text{obj}(x, y) = \sum_{i,j}^{\text{patch}} (\mathrm{P}[i, j] - \mathrm{I}_2[i+x, j+y])^2

\end{equation}

$$

Since the objective function is the square of residuals, Gauss-Newton provides a formulation for an iterative update scheme

$$

\begin{equation}

\mathbf{x}^{i+1} = \mathbf{x}^{i} - (\mathbf{J}^T \mathbf{J})^{-1}\mathbf{J}^T\mathbf{r}

\end{equation}

$$

where $(\mathbf{r}$) are the stacked pixel residuals and $(\mathbf{J} = \partial/\partial_{x, y}\mathbf{ r}$). The derivative of a residual,

$$

\begin{equation}

\mathrm{P}[i, j] - \mathrm{I}_2[i+x, j+y]

\end{equation}

$$

simplifies to the derivative of the only part that depends on $(x, y$), namely $(\partial/\partial_{x, y}\mathrm{ I}_2$). This can be approximated by applying an edge operator, e.g Sobel, to $(\mathrm{I}_2$) in each direction.

$$

\begin{align}

\mathbf{r} =& \begin{pmatrix}

\text{flatten}(\mathrm{P}[i, j] - \mathrm{I}_2[i+x, j+y])

\end{pmatrix}^T \\

\mathbf{J} =& \begin{pmatrix}

-\text{flatten}(\partial/\partial_x \mathrm{ I}_2) \\

-\text{flatten}(\partial/\partial_y \mathrm{ I}_2)

\end{pmatrix}^T

\end{align}

$$

For the sake of completeness, it is worth mentioning that in practice $(\partial/\partial_{x, y}\mathrm{ P}$) is used in place of $(\partial/\partial_{x, y}\mathrm{ I}_2$). Although the maths is technically incorrect, since $(\mathrm{P}$) does not depend on $(x, y$), it turns out that locally during each iteration the role of the patch and the image can be reversed. The trick is to then do the opposite of whatever the Gauss-Newton update step would suggest. The advantages are that $(\partial/\partial_{x, y}\mathrm{ P}$) need only be computed once, and that using a fixed Jacobian seems to make the algorithm more stable.

At any rate, the previous formulation is for a single patch. All patches can be combined into a single iteration as below.

$$

\begin{equation}

\begin{bmatrix}

\mathbf{J}_1 && && && \\

\vdots && && && \\

&& \mathbf{J}_2 && && \\

&& \vdots && && \\

&& && \ddots && \\

&& && \ddots && \\

&& && && \mathbf{J}_n \\

&& && && \vdots

\end{bmatrix}\begin{bmatrix}

\Delta x_1 \\

\Delta y_1 \\

\Delta x_2 \\

\Delta y_2 \\

\vdots \\

\Delta x_n \\

\Delta y_n

\end{bmatrix} = \begin{bmatrix}

\mathbf{r}_1 \\

\vdots \\

\mathbf{r}_2 \\

\vdots \\

\vdots \\

\vdots \\

\mathbf{r}_n \\

\vdots

\end{bmatrix}

\end{equation}

$$

To convert from this form, which individually aligns each patch, to one which globally aligns all patches with respect to a transformation $(\mathbf{T}$), we apply the chain rule.

$$

\begin{equation}

\mathbf{J}_\mathbf{T} = \frac{\partial\mathrm{I}_2}{\partial\mathbf{T}} = \frac{\partial\mathrm{I}_2}{\partial \langle x,y \rangle}\frac{\partial \langle x,y \rangle}{\partial\mathbf{T}}

\end{equation}

$$

Computing this is an exercise in projective geometry, but since our two images are separated predominantly by forward motion, let us further simplify things by instead only computing the image Jacobian with respect to this motion, which we will call $(\mathbf{J}_z$). We now have a function of a single variable, $(\text{ obj}(z)$), which can be plotted.

There is one last implementation detail, which is that the alignment algorithm is performed at multiple image scales. The smoother, lower curves show $(\text{obj}(z)$) plotted at reduced image sizes, where noise is spread across adjacent pixels. The least squares update, as computed using Gauss-Newton, is shown at the bottom.

Interestingly, despite multiple local minima, each estimated gradient crosses zero exactly once. This concludes the portion of the blog where I claim to know exactly what's going on.

I believe the key is that $(\mathbf{J}_z$) has been built using only estimates of image gradients, which have all been averaged together to compute $(\partial/\partial_z\text{ obj}(z)$). This has produced a gradient function which has been smoothed in a non-trivial way. The take-away is that this method of estimation seems substantially more robust than numerical differentiation.

It is important to be aware of this phenomena because of how many root finding algorithms work. Some, for example Levenberg-Marquardt, will reject updates which worsen the objective, effectively becoming stuck in a local minimum despite a strong gradient away from the current value of $(z$).

Other algorithms will opt to numerically differentiate with respect to a subset of variables being optimized, even if a method of computing the Jacobian is provided by the user. The rational is that numeric differentiation can be cheaper than computing the full Jacobian, however this post has shown that in some cases, the Jacobian may provide a better estimate of the gradient.

I was unable to get any off-the-shelf least squares solvers to find the correct minimum. I've currently settled for vanilla Gauss-Newton with a termination criterion for when the update direction reverses.

]]>Analytically finding the smallest distance between a point and an ellipse boils down to solving a quartic equation. The easiest way to see this is by considering that a circle can intersect an ellipse at most four times. In this view, the problem is to find the circle that intersects

]]>Analytically finding the smallest distance between a point and an ellipse boils down to solving a quartic equation. The easiest way to see this is by considering that a circle can intersect an ellipse at most four times. In this view, the problem is to find the circle that intersects the ellipse once. If formulated as a quartic, two of the roots will be shared at the intersection, and the other two roots will be complex.

Although a closed solution to the quartic does exist, it is unsurprising that iterative methods would be easier to implement and less prone to numeric instabilities. Generic root finders will work, but a robust implementation would be advised. This StackOverflow post shows what can go wrong when using a naive implementation of Newton's Method. Specialised iterative methods also exist, such as described by Maisonobe.

The paper's premise is that an oversized circle will intersect the ellipse at least two times. The iteration to minimise the distance between point `a`

on the ellipse and `p`

is as follows:

- Choose a point
`a`

on the ellipse - Compute the distance
`r`

between`a`

and`p`

- On the circle centred at
`p`

with radius`r`

, find the other intersection`b`

- Set
`a'`

to the midpoint of`a`

and`b`

on the ellipse - Return to step 2

Visually `a'`

is obviously a very good innovation over `a`

, however computing `b`

and `a'`

is non trivial. In brief, the quartic still has 4 roots, but because the initial guess `a`

is a solution, the quartic can be reduced to a cubic. The cubic is solved to find `b`

. An additional quadratic is then solved to find `a'`

. As noted, this complexity buys a very good estimate for `a'`

, and therefore the paper's claim that this algorithm is robust and converges quickly is very believable. What follows is an alternative, much simpler computation for `a'`

.

Consider the previous algorithm when the ellipse is a circle. `a'`

then trivially lies on the line between `c`

and `p`

. It can also be seen that `a'`

is optimal. The idea will become to locally approximate the ellipse as a circle, but first we must touch on curvature.

As a smooth curve, every point on the ellipse has a perpendicular normal. Additionally, every point has a radius of curvature. Without going into specifics, a tighter curve will result in a smaller radius of curvature, while a shallow curve will result in a larger radius of curvature. The centre of curvature of point `e`

on the ellipse is the centre of the circle with radius equal to the radius of curvature and perpendicular to the ellipse at `e`

. In other words, it is the centre of the circle that locally approximates the ellipse at `e`

. The centre of curvatures can be plotted to form the evolute of the ellipse, which Wolfram and Wikipedia are better suited to explain.

It's now time to define the curves mathematically. Although everything can derived in Cartesian coordinates, a parametric definition avoids infinities about the vertices. For the ellipse use the standard definition.

```
x = a * cos(t)
y = b * sin(t)
```

The ellipse evolute also has a parametric form.

```
ev_x = (a*a-b*b)/a * cos^3(t)
ev_y = (b*b-a*a)/b * sin^3(t)
```

This is extremely useful for approximating the ellipse as a circle, since we know the centre of curvature for any `t`

. Recall that when the ellipse is a circle, `a'`

lies on the line between the circle centre `c`

and `p`

. Now modify this claim slightly, to state that `a'`

lies on the line between the centre of curvature `c`

and `p`

. This is not as exact as the original `a'`

from the paper since we've approximated the ellipse as a circle, but the beauty is that the closer we get to the optimal `a'`

, the better the approximation becomes. Therefore I would expect similar convergence and robustness properties.

Intersecting a line with an ellipse is a much simpler problem, however it still requires solving a quadratic. Since we're already estimating, an exact solution isn't required. From the previous figure, define two vectors `r`

and `q`

. Very loosely speaking, `r`

is the radius of curvature and `q`

intersects the ellipse where we would like the radius of curvature to have been. Since we are still using our circle approximation, we can compute the arc length between `a`

and `a'`

as the angle between `r`

and `q`

times the radius of curvature.

```
(r x q)
sin(Δc/|r|) ≈ -------
|r||q|
```

Additionally, since `Δc`

is small, we could further approximate by dropping the sine.

Now we would like to know how much to vary `t`

by to achieve the same arc length delta on the ellipse. With a bit of calculus:

```
dx/dt = -a * sin(t)
dy/dt = b * cos(t)
dc/dt = sqrt((dx/dt)^2 + (dy/dt)^2)
Δc/Δt ≈ sqrt(a^2 * sin^2(t) + b^2 * cos^2(t))
Δt ≈ Δc / sqrt(a^2 * sin^2(t) + b^2 * cos^2(t))
Δt ≈ Δc / sqrt(a^2 - a^2 * cos^2(t) + b^2 - b^2 * sin^2(t))
Δt ≈ Δc / sqrt(a^2 + b^2 - x^2 - y^2)
a'x = a * cos(t + Δt)
a'y = b * sin(t + Δt)
```

This proves to be a good enough approximation, but it does suffer at the vertices (pointy ends) of the ellipse. Although in my experiments convergence was not affected, in the code listing I choose to confine `t`

to the first quadrant and then correct the signs of the final result. For clarity, some obvious optimisations have been left unimplemented.

```
def solve(semi_major, semi_minor, p):
px = abs(p[0])
py = abs(p[1])
t = math.pi / 4
a = semi_major
b = semi_minor
for x in range(0, 3):
x = a * math.cos(t)
y = b * math.sin(t)
ex = (a*a - b*b) * math.cos(t)**3 / a
ey = (b*b - a*a) * math.sin(t)**3 / b
rx = x - ex
ry = y - ey
qx = px - ex
qy = py - ey
r = math.hypot(ry, rx)
q = math.hypot(qy, qx)
delta_c = r * math.asin((rx*qy - ry*qx)/(r*q))
delta_t = delta_c / math.sqrt(a*a + b*b - x*x - y*y)
t += delta_t
t = min(math.pi/2, max(0, t))
return (math.copysign(x, p[0]), math.copysign(y, p[1]))
```

One remaining question is why `t`

is initialised to a constant instead of `atan2(py, px)`

. This is only a good guess if `p`

is outside of the ellipse. On the inside it is optimally bad. If `p`

is inside the ellipse initialise to `atan2(px, py)`

. It turns out that a close to optimally bad guess will require additional iterations to move away from the initial guess in a random fashion. An optimally bad guess will get stuck. I've shown that for rather detailed reasons, the algorithm converges as long as you don't start near the optimally bad guess. You will be fine using `atan2`

if you correctly case for being inside or outside of the ellipse.

The plots show good convergence after three iterations over a wide range of eccentricities. However, at extreme eccentricities the algorithm will eventually fail. At some cutoff the ellipse should be considered a line.

The code for the plots is available on GitHub.

]]>