This is an interesting problem if one insists on a solution constructable using purely geometric operations.
For simplicity in describing the solution, we will choose a coordinate system such that
the two foci $A, B$ are $(1,0)$ and $(0,1)$ respectively. We will also assume the circle
$\mathscr{C}$ we are dealing with is centered at $C = (\alpha,\beta)$ in the first
quadrant and the radius of $\mathscr{C}$ is $r$.
It is known that the family of ellipses and hyperbolas having $A$ and $B$ as foci are given by:
$$\frac{x^2}{1+\lambda} + \frac{y^2}{\lambda} = 1\quad\text{ and is }
\begin{cases}
\text{ an ellipse, }&\text{ if } \lambda > 0\\
\text{ a hyperbola, } &\text{ if } 0 > \lambda > -1
\end{cases}\tag{*0}$$
For a point $p = (x,y)$ on this ellipse/hyperbola, the normal is in the direction $(\frac{x}{1+\lambda}, \frac{y}{\lambda})$. i.e. the normal line passing through $(x,y)$ is given by:
$$\left\{\;\left( x ( 1 + \frac{t}{1+\lambda} ),\,y ( 1 + \frac{t}{\lambda} ) \right) : t \in \mathbb{R}\;\right\}$$
If this ellipse/hyperbola is tangent to any circle centered at $C$ at point $p$, the above normal line with pass through $C = (\alpha,\beta)$ and one can show that $(x,y)$ lies on
a cubic curve $\mathscr{P}$:
$$(\alpha y - \beta x )\left(x(x-\alpha) + y(y-\beta)\right) - (x-\alpha)(y-\beta) = 0\tag{*1}$$
This is the red curve depicted in above picture.
For $r$ not too large, $\mathscr{C}$ and $\mathscr{P}$ in general will intersect at 4 points. Two of them are "tangent" points to two ellipses and the other two are tangent points to two hyperbolas. All these four ellipses/hyperbolas are having $A, B$ as their foci. For these four points, one can simplify $(*1)$ and shows that they lie on a conic
$\mathscr{Q}$:
$$(\alpha y - \beta x)(\alpha x + \beta y - \gamma) - (x-\alpha)(y-\beta) = 0\tag{*2}$$
where $\gamma = \alpha^2 + \beta^2 - r^2$. This is the cyan curve containing
$C, F, G, J, K, L$ in above picture.
Some of these points are relatively easy to construct.
- $C$ - of course we know where $C$ is. In fact, one can show that $\mathscr{Q}$ is tangent
to the line $OC$ at point $C$. Counting multiplicity, this gives "two" points to determine $\mathscr{Q}$.
- If one draws a circle using $OC$ as diameter, it will intersect $\mathscr{C}$ at two points $D$ and $E$. The line $DE$ in turn intersects the horizontal line containing $C$ at $F$ and the vertical line containing $C$ at $G$.
A conic is determined by any five points on it. If for whatever means we get an extra
point on $\mathscr{Q}$, then we can construct $\mathscr{Q}$ and its intersections with $\mathscr{C}$, e.g. the point $L$ in the above picture. These are the tangent points we want.
One can show that the conic $\mathscr{Q}$ intersects the $x$-axis at $J = (\mu,0)$
and $K = (1/\mu,0)$ where $\mu$ is a root of the quadratic equation: $$\alpha x^2 - (\gamma +1) x + \alpha = 0\tag{*3}$$
If one doesn't mind a little bit of algebra, we are done. Otherwise, we can convert this
formula to a geometric construction of $J$ and $K$ in 3 steps:
- Construct an auxiliary hyperbola $xy + 1 = 0$ (the olive green hyperbola in above picture). The foci of this hyperbola are $(\sqrt{2},-\sqrt{2})$ and $(-\sqrt{2},\sqrt{2} )$ and it contains the point $(1,-1)$. All of them can be constructed using compass and ruler start from $O$ and $A$.
- Construct an auxiliary horizontal line $\beta y + 1 = 0$ (the blue horizontal line containing $I$ in above picture). One first constructs a line segment joining $Y = (0,\beta)$ to $B = (-1,0)$. One then constructs another line segment parallel to it
start at $Y' = (0,1)$. This new segment will intersect the $x$-axis at $B' = (-\frac{1}{\beta},0)$. Another quarter circle will bring us to the point $I = (0,-\frac{1}{\beta})$ and the horizontal line at $I$ is the one we need.
- Let $H$ be the intersection of the line $DE$ with the blue line in previous step.
Project it onto $x$-axis as $H_1$. Now construct a ray start from $H_1$ at $-135^{\circ}$ with respect to the $x$-axis. This ray will intersect the "olive green" hyperbola at $H_2$ and $H_3$. If one projects $H_2$ and $H_3$ back to the $x$-axis, we get $K$ and $J$ respectively.
Update
Since Stefan asks, here are some details how to derive the various equations.
If the normal passes through $C = (\alpha,\beta)$, we have a $t$ such that:
$$x (1 + \frac{t}{1+\lambda}) = \alpha \quad\text{ and }\quad y (1+ \frac{t}{\lambda}) = \beta$$
Eliminating $t$ give us:
$$\begin{align}(1+\lambda)(\frac{\alpha}{x}-1) = \lambda (\frac{\beta}{y} - 1 )
\iff & (1+\lambda)(x-\alpha)y - \lambda(y-\beta)\alpha = 0\\
\iff & (x-\alpha)y - \lambda( \alpha y - \beta x ) = 0
\end{align}$$
This gives us
$$
\lambda = \frac{(x - \alpha)y}{\alpha y - \beta x}
\quad\text{ and }\quad
1 + \lambda = \frac{(y - \beta)x}{\alpha y - \beta x}
$$
Substitute this back into $(*0)$, we get $(*1)$:
$$\begin{align}
& (\alpha y - \beta x)\left(\frac{x}{y-\beta} + \frac{y}{x-\alpha}\right) = 1\\
\iff & (\alpha y - \beta x)\left(x(x-\alpha)+y(y-\beta)\right) - (x-\alpha)(y-\beta) = 0
\end{align}$$
If $(x,y)$ also lies on $\mathscr{C}$, the $2^{nd}$ factor in above expression can be simplified as:
$$\begin{align}
x(x-\alpha)+y(y-\beta) = & (x-\alpha)^2 + (y-\beta)^2 + \alpha (x-\alpha)+\beta(y-\beta)\\
= & \alpha x + \beta y - (\alpha^2 + \beta^2 - r^2)\\
= & \alpha x + \beta y - \gamma
\end{align}$$
Substituting this into $(*1)$ gives us the equation of $\mathscr{Q}$ in $(*2)$:
$$(\alpha y - \beta x)(\alpha x + \beta y - \gamma) - (x-\alpha)(y-\beta) = 0$$
Now $(*3)$ is really the equation of intersecting $\mathscr{Q}$ with the $x$-axis.
We just set $y = 0$ in $(*2)$, eliminate the common factor $\beta$ and we are done.
I assume that the variables x and y in the program are
arrays of numbers, because without that assumption it makes no sense
to fit an ellipse to them.
Then for each pair of values $(x_i,y_i)$
you compute $\theta_i = \mathrm{atan2}(y_i,x_i).$
This appears to be a fatal flaw, because the angle from
the origin to an $(x,y)$ point after the circle has been distorted
into an ellipse and translated away from the origin
has very little to do with the parameter $\theta$ that you need
in order to make the rest of your formulas work.
If you knew (or could reliably guess)
the angle $\theta$ corresponding to a point on the circle
before all the transformations that create the ellipse were applied,
then I think you could use calculations like the ones you showed
in order to convert the $(x,y)$ coordinates of the point on the ellipse
to which your original point was mapped--but on the other hand
if you know $\theta$ then you can plot the point on the circle without
even looking at the coordinates on the ellipse.
I suggest trying a simple example, such as the ellipse generated
by these parameterized equations:
\begin{align}
x &= 2\cos(\theta) + 5, \\
y &= \sin(\theta) + 2. \\
\end{align}
Plot these points for a few dozen equally-spaced of $\theta$
ranging from $0$ to $2\pi,$
confirm that you get an ellipse,
then give these sets of coordinates to your function
and see if it gives you back a circle.
I will be surprised if the result makes any sense at all.
To do the algorithm correctly, as noted in a comment you just have to
use the information from fitEllipse to create a matrix that
scales the plane by a factor of $1/a$ parallel to the first axis of the
ellipse and by a factor of $1/b$ parallel to the second axis of the ellipse.
Translate all the points of the ellipse so that they are centered around
the origin, then use the matrix to scale them onto a circle.
Best Answer
If point $O$ is to be mapped to the centre of the circle, let $EF$ be the diameter of the ellipse through $O$. Draw tangents $AB$, $CD$ through $F$ and $E$ and draw $GH$ through $O$ parallel to both tangents. Draw then tangents $BC$, $DA$ through $G$ and $H$. There is a homography mapping the ellipse to the circle, point $O$ to the center $O'$ of the circle and trapezoid $ABCD$ to square $A'B'C'D'$.
To find the transformed of a point $P$ inside the ellipse, you can exploit the invariance of cross-ratios in a homography. For instance, draw line $OP$ intersecting two opposite sides of the trapezoid at $R$ and $Q$. From $(A,B;F,Q)=(A',B';F',Q')$ you can find $Q'$ and then from $(R,Q;O,P)=(R',Q';O',P')$ you can find $P'$.
As this is a linear transformation (in homogeneous coordinates) you can also find a transformation matrix from four couples of corresponding points (e.g. points $EFGH$ and $E'F'G'H'$).