Circle tangent to rotated ellipse and horizontal line

Question

I would like to find the position for the center of a circle $(x_0, y_0)$ that is tangent to both an ellipse and a horizontal line. The ellipse is positioned at $(0,0)$ and is defined by major axis $a$ and minor axis $b$ and can be rotated by $\alpha$. The line $h$ is horizontal and can be moved on the $y$-axis.

The Python snippet in this answer works well but doesn't allow for the ellipse to be rotated.

Edit: I'm a technical artist in the game industry working on a procedural tool that involves calculating the center of a circle to create an almost accurate 3D model. I need to understand the math behind it in order to generate code to solve it. My math knowledge ranges from high school level to math applicable to basic 3D graphics.

Answer 1

From

$$ \left(\frac xa\right)^2+\left(\frac yb\right)^2-1=0 $$

changing of coordinates by rotating $-\theta$ we have

$$ \frac{(X \cos (\theta )+Y \sin (\theta ))^2}{a^2}+\frac{(Y \cos (\theta )-X \sin (\theta ))^2}{b^2}-1=0 $$

or

$$ \left(\frac{\cos ^2(\theta )}{a^2}+\frac{\sin ^2(\theta )}{b^2}\right)X^2+\left(\frac{1}{a^2}-\frac{1}{b^2}\right) \sin (2 \theta )XY+\left(\frac{\sin ^2(\theta )}{a^2}+\frac{\cos ^2(\theta )}{b^2}\right)Y^2-1=0 $$

and the normal vector $\vec n$ is given by

$$ \vec n = \left\{ \begin{array}{c} \frac{ \cos ^2(\theta ) (Y \tan (\theta )+X)}{a^2}+\frac{ \sin (\theta ) (X \sin (\theta )-Y \cos (\theta ))}{b^2} \\ \frac{ \sin (\theta ) (X \cos (\theta )+Y \sin (\theta ))}{a^2}+\frac{ \cos (\theta ) (Y \cos (\theta )-X \sin (\theta ))}{b^2} \\ \end{array} \right. $$

now changing to polar coordinates $$ E(t)=\cases{ X = \rho(t)\cos t\\ Y = \rho(t)\sin t } $$

we have

$$ \cases{ \rho ^2(t) \left(\frac{\cos ^2(t-\theta )}{a^2}+\frac{\sin ^2(t-\theta )}{b^2}\right)-1=0\\ \vec n(t) = \left\{ \begin{array}{c} \frac{\rho(t) \left(\left(b^2-a^2\right) \cos (t-2 \theta )+\left(a^2+b^2\right) \cos (t)\right)}{a^2 b^2} \\ \frac{\rho(t) \left(\left(a^2+b^2\right) \sin (t)+(a-b) (a+b) \sin (t-2 \theta )\right)}{a^2 b^2} \\ \end{array} \right. } $$

now taking

$$ \rho(t) = \frac{ab\sqrt{2}}{\sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}} $$

and substituting into

$$ \mathcal{E}_r(t)=(\rho(t)\cos t,\rho(t)\sin t) + r\frac{\vec n}{\|\vec n\|} $$

where $\mathcal{E}_r(t)$ is $E(t)$ expanded by $r$, we get at

$$ \mathcal{E}_r(t)=\frac{ab\sqrt{2}}{\sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}}\left(\cos t, \sin t\right)+r\frac{\vec n}{\|\vec n\|} $$

where

$$ \vec n = \left\{ \begin{array}{c} \frac{\sqrt{2} \left(\left(b^2-a^2\right) \cos (t-2 \theta )+\left(a^2+b^2\right) \cos (t)\right)}{a b \sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}} \\ \frac{\sqrt{2} \left(\left(a^2+b^2\right) \sin (t)+(a-b) (a+b) \sin (t-2 \theta )\right)}{a b \sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}} \\ \end{array} \right. $$

now, if the horizontal line is given by $Y=h$, the tangent circle center is located at the solutions for

$$ (0,1)\cdot \mathcal{E}_r(t) = h + r $$

where $r$ is the circle radius.

NOTE

Assuming the parametric values $a=1,b=3,h=1,r=0.5,\theta = -\frac{\pi}{6}$ we have

$$ (0,1)\cdot \mathcal{E}_r(t) - (h + r) = \frac{3 \sin (t)}{\sqrt{2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5}}+\frac{3 \sin (t)+2 \sqrt{3} \cos (t)}{3 \sqrt{2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5} \sqrt{\frac{4 \left(3 \sin (t)+2 \sqrt{3} \cos (t)\right)^2}{9 \left(2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5\right)}+\frac{4 \left(2 \sqrt{3} \sin (t)+7 \cos (t)\right)^2}{9 \left(2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5\right)}}}-\frac{3}{2}=0 $$

to solve for $t$. Follows a plot shoving in blue $E(t)$ in dashed light brown $\mathcal{E}_{\frac 12}(t)$ also in dashed red $Y=1+\frac 12$. The solution for $t$ can be obtained by binary search.

Follows a python script to find the solution angles and circle centers.

import numpy as np
import matplotlib.pyplot as plt


n     = 150
x0    = np.linspace(0,np.pi,n)
err   = 0.00001
a     = 1
b     = 3
h     = 1
r     = 1/2
theta =-np.pi/6

def expanded_ellipse(t):
    st   = np.sin(t)
    ct   = np.cos(t)
    stt  = np.sin(t+theta)
    ctt  = np.cos(t+theta)
    c2tt = np.cos(2*(t+theta))
    st2t = np.sin(t+2*theta)
    ct2t = np.cos(t+2*theta)
    den1 = np.sqrt((b*ctt)**2 + (a*stt)**2)
    den2 = np.sqrt(2)*np.sqrt(a**4 + b**4 + (b**4 - a**4)*c2tt)
    numx = ((a**2 + b**2)*ct + (b**2 - a**2)*ct2t)*r
    numy = ((a**2 + b**2)*st + (a**2 - b**2)*st2t)*r
    return [a*b*ct/den1 + numx/den2, a*b*st/den1 + numy/den2]


def residual(t):
    st   = np.sin(t)
    stt  = np.sin(t+theta)
    ctt  = np.cos(t+theta)
    c2tt = np.cos(2*(t+theta))
    st2t = np.sin(t+2*theta)
    num  = r*((a**2 + b**2)*st + (a**2 - b**2)*st2t)
    den  = np.sqrt(2)*np.sqrt(a**4 + b**4 + (b**4 - a**4)*c2tt)
    return a*b*st/np.sqrt((b*ctt)**2+(a*stt)**2) + num/den - h - r


def minus_residual(t):
    return -residual(t)


def binary_search(f, domain, MAX = 1000):
    start, end = domain
    if start >= end:
        raise ValueError("Domain is empty")
    mid = (start + end) / 2
    fmid = f(mid)
    step = 0
    while abs(fmid) > err and step < MAX:
        if fmid < 0:
            start = mid
        else:
            end = mid
        mid = (start + end) / 2
        fmid = f(mid)
        step += 1
    return round(mid, 7)


y0 = []

for x in x0:
    y0.append(residual(x))

plt.xlim([0,np.pi])
plt.ylim([min(y0),max(y0)])
plt.plot(x0,y0,'-',color = 'black')
plt.plot(x0,n*[0])

t1 = binary_search(residual,[1,1.7])
t2 = binary_search(minus_residual,[2,3])

print(t1, expanded_ellipse(t1))
print(t2, expanded_ellipse(t2))

NOTE

The binary_search procedure needs that the function be increasing in the searching domain: so we have residual and minus_residual.

Answer 2

Assuming the vector $r = (x_1, x_2)$ is the tangency point between the ellipse and the circle, and vector $ C = (x_3, h + R)$ is the center of the tangent circle, and $R$ is the given radius of the circle, then we want:

$(r - r_0)^T Q_e (r - r_0) = 1 $

$Q_e (r - r_0) = K (r - C) $

$(r - C)^T (r - C) = R^2 $

$Q_e$ is the matrix of the ellipse, and it is computed as follows:

If $a$ the semi-major axis, and $b$ is the semi-minor axis, then define the diagonal matrix $D$ as follows

$D = \begin{bmatrix} \dfrac{1}{a^2} && 0 \\ 0 && \dfrac{1}{b^2} \end{bmatrix} $

And let $R$ be the $2 \times 2 $ rotation matrix corresponding to the rotation of the major axis counter clockwise by an angle $\theta$, then,

$R = \begin{bmatrix} \cos \theta && - \sin \theta \\ \sin \theta && \cos \theta \end{bmatrix} $

Then,

$ Q_e = R D R^T $

We can eliminate the constant $K$ from the second equation as follows:

$ V_x (x_2 - (h + R)) - V_y (x_1 - x_3) = 0 $

where $V = Q_e (r - r_0)$

To solve these equations, I used the Newton-Raphson multivariate root finder. This method is summarized by its recursion

$ \mathbf{x}_{k+1} = \mathbf{x}_k - J^{-1}(\mathbf{x}_k) \mathbf{f}_k , \ k = 0, 1, 2, \dots$

where vector $\mathbf{x} = [ x_1, x_2, x_3 ]^T $, and vector $\mathbf{f}$ is the vector of nonlinear functions whose simultaneous root is to be found. $\mathbf{f} = [f_1, f_2, f_3]^T$, where

$f_1 = (r - r_0)^T Q_e (r - r_0) - 1$

$f_2 = V_x ( x_2 - (h + R) ) - V_y (x_1 - x_3) $

$f_3 = (r - C)^T (r - C) - R^2 $

The Jacobian of this vector is the matrix $J = [J_{ij}]$ , where $J_{ij} = \dfrac{\partial f_i }{ \partial x_j } $

Depending on the initial guess of $(x_1, x_2, x_3)$, we'll get a different solution.

There are 4 possible solutions, two circles tangent to the ellipse from the outside, and two from the inside.

Here is a sample run for an ellipse with semi-major axis length of $a =3$ and semi-minor axis length of $b=1$, and an angle of rotation of the major axis from the positive x-axis of $120^\circ$. In addition, $h = 1$ and $R = 0.5$. The Newton-Raphson method converged to the solution in just $5$ iterations (this is because the initial guess of the solution was a good one and the fact that Newton-Raphson is a robust method).

And here is the full script (Microsoft Excel Visual Basic for Applications) that I wrote to implement this solution.

Public Sub place_circle_tangent_to_given_ellipse_and_line()
Dim diag(2, 2), qe(2, 2), r(2, 2), rt(2, 2) As Double
Dim r0(2), c(2) As Double
Dim np(2), rp(2) As Double
Dim x(3), y(3), jac(3, 4) As Double
Dim v(2, 2) As Double

r0(1) = 0
r0(2) = 0

r(1, 1) = Cos(2 * p / 3)
r(2, 2) = r(1, 1)
r(2, 1) = Sin(2 * p / 3)
r(1, 2) = -r(2, 1)

Call transpose(2, 2, r, rt)

diag(1, 1) = 1 / 3 ^ 2
diag(2, 2) = 1

Call multiply(2, 2, 2, r, diag, qe)
Call multiply(2, 2, 2, qe, rt, qe)

h = 1

xrad = 0.5

x(1) = 1
x(2) = 1.5
x(3) = 2

' This initial guess should produce the solution circle tangent from the 
outside on the right.
success = 0
For ic = 1 To 50

Call find_f_ellipse_circle(r0, qe, h, xrad, x, y)

If norm(y, 3) < 0.0000001 Then

MsgBox ("Converged in " + Str(ic - 1))
    success = 1
    Exit For

End If

Call find_jac_ellipse_circle(r0, qe, h, xrad, x, jac)

For i = 1 To 3
    jac(i, 4) = -y(i)
Next i

Call reduce_matrix(jac, 3, 4, 0.00000001, det)

For i = 1 To 3
    x(i) = x(i) + jac(i, 4)
Next i

Next ic

If success = 0 Then
   MsgBox ("Newton-Raphson did not converge")
   Exit Sub
End If

c(1) = x(3)
c(2) = h + xrad

For i = 1 To 2
    diag(i, i) = 1 / Sqr(diag(i, i))
Next i

Call multiply(2, 2, 2, r, diag, v)


k = 1
narcs = 0
ipass = 2

Call draw_line_segment(0, -2, h, 2, h, 4)

Call draw_ellipse(r0, v, 5)

v(1, 1) = xrad
v(2, 1) = 0

v(1, 2) = 0
v(2, 2) = xrad

Call draw_ellipse(c, v, 46)

Call plot_curve(0, 1, 0, 0, 0, "", "")

End Sub

---

Public Sub find_f_ellipse_circle(r0, qe, h, xrad, x, y)
Dim c(2), u(2), v(2), w(2) As Double
Dim e1(2), e2(2), M(2, 2), N(2, 2) As Double

c(1) = x(3)
c(2) = h + xrad

For i = 1 To 2
    u(i) = x(i) - r0(i)
    w(i) = x(i) - c(i)
Next i

Call multiplyv(2, 2, qe, u, v)

y(1) = xtqy(2, qe, u, u) - 1

e1(1) = 1
e2(2) = 1

' Matrix M is there to implement the cross product

For i = 1 To 2
For j = 1 To 2
    M(i, j) = e1(i) * e2(j) - e1(j) * e1(i)
Next j
Next i

Call multiply(2, 2, 2, qe, M, N)

'(r - r0)^T Q M ( r - C ) = u^T N w

y(2) = xtqy(2, N, u, w)

y(3) = dotn(w, w, 2) - xrad ^ 2

End Sub

---

Public Sub find_jac_ellipse_circle(r0, qe, h, xrad, x, jac)
Dim c(2), u(2), v(2), w(2) As Double
Dim e1(2), e2(2), M(2, 2), N(2, 2), Nt(2, 2) As Double
Dim Nsym(2, 2) As Double

Dim d0(2) As Double, d1(2) As Double
Dim u1(2), u2(2) As Double

For i = 1 To 3
For j = 1 To 3
    jac(i, j) = 0
Next j
Next i

e1(1) = 1
e2(2) = 1

For i = 1 To 2
For j = 1 To 2
    M(i, j) = e1(i) * e2(j) - e1(j) * e1(i)
Next j
Next i

Call multiply(2, 2, 2, qe, M, N)
Call transpose(2, 2, N, Nt)

Call madd(2, 2, N, Nt, Nsym)
Call mscale(2, 2, Nsym, 0.5, Nsym)


c(1) = x(3)
c(2) = h + xrad

For i = 1 To 2
    u(i) = x(i) - r0(i)
    w(i) = x(i) - c(i)
Next i

Call multiplyv(2, 2, qe, u, v)

For j = 1 To 2
    jac(1, j) = 2 * v(j)
Next j


Call multiplyv(2, 2, N, c, u1)
Call multiplyv(2, 2, Nt, r0, u2)

Call multiplyv(2, 2, Nsym, x, d1)
Call vscale(2, d1, 2, d1)

For j = 1 To 2
    jac(2, j) = d1(j) - u1(j) - u2(j)
Next j

' now, with respect to C :  
' (r - r0)^T N (r - C)  = (r - r0)^T N r - (r - r0)^T N C

' so the gradient with respect to C(1) is the x-entry of  - (r - r0)^T N

Call multiplyv(2, 2, Nt, u, u1)

jac(2, 3) = -u1(1)


For j = 1 To 2
    jac(3, j) = 2 * w(j)
Next j

jac(3, 3) = -2 * dotn(w, e1, 2)

End Sub