Calculate the inscribed circle between one line and two circles

Question

I am currently working on finding a calculation that will help me determine the values of $C_x$, $C_y$ and $C_r$ as shown in the image below.

The goal is for the unknown red circle to touch the line and be tangent to the two green circles.

This website explains the situation as pointed out by Intelligenti pauca

For example, the given values are:
C1r = 2
C1x = 2
C1y = 8
C2r = 2.5
C2x = 11
C2y = 9
P1x = 5
P1y = 0
P2x = 12
P2y = 3
The outcome should be somewhere around:
Cr = 3.671
Cx = 6.596
Cy = 4.678

The worked out result in C# code (based on the accepted answer)

    public static bool GetCircleFromLCC(out double C_x, out double C_y, out double C_radius)
    {
        // -- Input values --
        double C1_radius = 2;
        double C1_x = 2;
        double C1_y = 8;
        double C2_radius = 2.5;
        double C2_x = 11;
        double C2_y = 9;
        double P1_x = 5;
        double P1_y = 0;
        double P2_x = 12;
        double P2_y = 3;
    // -- Calculations --
    double u_x = P2_x - P1_x;
    double u_y = P2_y - P1_y;
    double norm = Math.Sqrt(u_x * u_x + u_y * u_y);
    u_x /= norm;
    u_y /= norm;
    double u1_x = -(u_x * P1_x + u_y * P1_y);
    double u1_y = -(-u_y * P1_x + u_x * P1_y);
    double u2_x = u_x * C1_x + u_y * C1_y;
    double u2_y = -u_y * C1_x + u_x * C1_y;
    double c1p_x = u2_x + u1_x;
    double c1p_y = u2_y + u1_y;
    u2_x = u_x * C2_x + u_y * C2_y;
    u2_y = -u_y * C2_x + u_x * C2_y;
    double c2p_x = u2_x + u1_x;
    double c2p_y = u2_y + u1_y;
    double a = 2.0 * (c1p_x - c2p_x);
    double b = 2.0 * (c1p_y - c2p_y + C1_radius - C2_radius);
    double c = Math.Pow(C1_radius, 2) - Math.Pow(C2_radius, 2) + Math.Pow(c2p_x, 2) - Math.Pow(c1p_x, 2) + Math.Pow(c2p_y, 2) - Math.Pow(c1p_y, 2);
    double a0 = Math.Pow(c1p_x, 2) + Math.Pow(c1p_y, 2) - Math.Pow(C1_radius, 2) + 2.0 * c / b * (c1p_y + C1_radius);
    double a1 = -2.0 * c1p_x + 2.0 * a / b * (C1_radius + c1p_y);
    double xdisc = Math.Pow(a1, 2) - 4.0 * a0;

    // Find quadratic roots
    if (xdisc &gt; 0.00000001)
    {
        double root1 = (-a1 - Math.Sqrt(xdisc)) / 2;
        C_radius = 1 / b * (-c - a * root1);
        if (Math.Sign(C_radius) == Math.Sign(c1p_y) &amp;&amp; Math.Sign(C_radius) == Math.Sign(c2p_y))
        {
            // Valid root, calculate the corresponding x,y coordinates
            u2_x = u_x * root1 + -u_y * C_radius;
            u2_y = u_y * root1 + u_x * C_radius;
            C_x = u2_x + P1_x;
            C_y = u2_y + P1_y;

            Console.WriteLine($&quot;Solution 1   C_x: {C_x:0.0##}   y: {C_y:0.0##}   radius: {C_radius:0.0##}&quot;);
            // Solution 1   x: 6.596   y: 4.678   radius: 3.671
        }
        double root2 = (-a1 + Math.Sqrt(xdisc)) / 2;
        C_radius = 1 / b * (-c - a * root2);
        if (Math.Sign(C_radius) == Math.Sign(c1p_y) &amp;&amp; Math.Sign(C_radius) == Math.Sign(c2p_y))
        {
            // Valid root, calculate the corresponding x,y coordinates
            u2_x = u_x * root2 + -u_y * C_radius;
            u2_y = u_y * root2 + u_x * C_radius;
            C_x = u2_x + P1_x;
            C_y = u2_y + P1_y;

            Console.WriteLine($&quot;Solution 2   C_x: {C_x:0.0##}   y: {C_y:0.0##}   radius: {C_radius:0.0##}&quot;);
            // Solution 2   x: -48.359   y: 336.121   radius: 329.963
            return true;
        }
    }

    C_x = 0;
    C_y = 0;
    C_radius = 0;
    return false;

}

This is a particular case of the Problem of Apollonius: https://www.cut-the-knot.org/Curriculum/Geometry/GeoGebra/LCC.shtml — Intelligenti pauca, Mar 03 '23 at 20:26
Thanks! Didn’t know what it was called, looks clear with the illustrations — Sander, Mar 04 '23 at 11:14

score 6 · Answer 1 · answered Mar 03 '23 at 16:13

For convenience, rotate the picture so the line is parallel to the $x$ axis, say $y = y_0$. Now you want the centre of your third circle to be a point $C = (x,y)$ such that $\text{dist}(C, C_1) = r_1 + y - y_1$ and $\text{dist}(C, C_2) = r_2 + y - y_2$, where $C_1 = (x_1, y_1)$, $C_2 = (x_2, y_2)$, and the given circles have radii $r_1$ and $r_2$. Expanding $(x - x_1)^2 + (y - y_1)^2 = (r_1 + y - y_1)^2$ and simplifying gives $x^2 - 2 x_1 x - 2 r_1 y + a_1 = 0$ where $a_1 = x_1^2 + 2 r_1 y_1 - r_1^2$. Similarly for the second circle you get $x^2 - 2 x_2 x - 2 r_2 y + a_2 = 0$. You can then eliminate $y$ and solve a quadratic equation to get $x$.

Do the calculations become more complex when dealing with a line that is not horizontal? — Sander, Mar 08 '23 at 08:00

Hosam Hajeer · Accepted Answer · 2023-03-08T15:58:06.493

Similar to Robert Israel solution, I am going to define a new reference frame with its origin at $P_1 = (5, 0) $, and its $u$ axis extending along $P_1 P_2 $

So the first coordinate vector is

$ \hat{u} = \dfrac{ P_2 - P_1 }{\| P_2 - P_1 \| } $

And the second coordinate vector is perpendicular to that

$ \hat{v} = [ - u_y , u_x ] $

Having defined the new coordinate system, express all the points with respect to this frame, this can be done as follows:

The relation between the coordinates in the $xy$ plane and $ uv $ plane is

$ P = P_1 + R Q $

So

$Q = R^T (P - P_1) $

where $R = [\hat{u} , \hat{v}] $

Apply this equation to $C_1, C_2, P_1, P_2 $ and name the new coordinates $C'_1, C'_2, P'_1, P'_2 $, so that

$ C'_1 = [u_1, v_1 ] $

$ C'_2 = [u_2, v_2 ] $

$P'_1 = [0, 0 ] $

$P'_2 = [ u_4, 0 ] $

Note: It is assumed that the circles are on the same side of the line segment. If $v_1 \lt 0 $ and $v_2 \lt 0 $ then: reverse the direction of vectors $\hat{u}$ and $\hat{v}$ and redefine the $R$ matrix, then redo the calculations for the computation of the coordinates in the $uv$ coordinate frame.

Next, let our center be expressed in the new reference frame be $ C' = (u, v) $, then it must satisfy

$ (u - u_1)^2 + (v - v_1)^2 = (r_1 + r)^2 $

$ (u - u_2)^2 + (v - v_2)^2 = (r_2 + r)^2 $

but the radius $r$ is equal to $v$, therefore,

$ (u - u_1)^2 + (v- v_1)^2 = (r_1 + v)^2 $

$ (u - u_2)^2 + (v - v_2)^2 = (r_2 + v)^2 $

Subtracting these two equations gives us

$ 2 u ( u_1 - u_2 ) + 2 v (v_1 - v_2 + r_1 - r_2) + r_1^2 - r_2^2 + u_2^2 - u_1^2 + v_2^2 - v_1^2 = 0 $

Solving for $v$ and substituting this equation $(3)$, one can determine the two solutions of $u$, and the corresponding values of $v$. The sign of $v$ must be the same as the sign of $v_1$ and $v_2$.

Having found the $(u, v)$ pair we can now calculate the $xy$ coordinates of the center from

$ P = P_1 + R Q $

I've implemented the above method into a VBA program on Excel.

Here's a listing of the subroutine.

Public Sub two_circles_and_a_line()
Dim r1, r2, c1(2), c2(2), p1(2), p2(2) As Double
Dim u(3), v(2) As Double
Dim r(2, 2) As Double, rt(2, 2) As Double
Dim c1p(2), c2p(2), p1p(2), p2p(2) As Double
Dim u1(2), u2(2) As Double
Dim roots(2) As Double
r1 = 2
r2 = 2.5
c1(1) = 2
c1(2) = 8
c2(1) = 11
c2(2) = 9
p1(1) = 5
p1(2) = 0
p2(1) = 12
p2(2) = 3
For i = 1 To 2
u(i) = p2(i) - p1(i)
Next i
Call normalize1(u)
v(1) = -u(2)
v(2) = u(1)
For i = 1 To 2
r(i, 1) = u(i)
r(i, 2) = v(i)
Next i
Call transpose(2, 2, r, rt)
' calculate coordinates in uv frame
p1p(1) = 0
p1p(2) = 0
Call multiplyv(2, 2, rt, p1, u1)
' negate u1
Call vscale(2, u1, -1, u1)
Call multiplyv(2, 2, rt, p2, u2)
Call vadd(2, u2, u1, p2p)
Call multiplyv(2, 2, rt, c1, u2)
Call vadd(2, u2, u1, c1p)
Call multiplyv(2, 2, rt, c2, u2)
Call vadd(2, u2, u1, c2p)
a = 2 * (c1p(1) - c2p(1))
b = 2 * (c1p(2) - c2p(2) + r1 - r2)
c = r1 ^ 2 - r2 ^ 2 + c2p(1) ^ 2 - c1p(1) ^ 2 + c2p(2) ^ 2 - c1p(2) 
^ 2
' now we have
' a u + b v + c = 0
' solve for v
'  v = 1/b ( - c - a u)
' substitute this into the distance equation
' (u - u1)^2 + (v - v1)^2 = (r1 + v )^2
' expanding and cancelling v^2 on both sides
' u^2 - 2 u1 u + u1^2 - 2 v v1 + v1^2 = r1^2 + 2 r1 v
' substitute  v = -1/b (a u + c )
' u^2 + u ( - 2 u1 + 2 a v1 / b + 2 a r1 / b ) + u1^2 + v1^2 - r1^2

2 v1 c / b + 2 r1 c / b = 0

a2 = 1
a1 = -2 * c1p(1) + 2 * a / b * (r1 + c1p(2))
a0 = c1p(1) ^ 2 + c1p(2) ^ 2 - r1 ^ 2 + 2 * c / b * (c1p(2) + r1)
Call find_quadratic_roots(a2, a1, a0, roots)
isol = 0
For i = 1 To roots(0)
u_ = roots(i)
v_ = 1 / b * (-c - a * u_)
If Sgn(v_) = Sgn(c1p(2)) And Sgn(v_) = Sgn(c2p(2)) Then
isol = isol + 1
' valid root
' calculate the corresponding x,y coordinates
u1(1) = u_
u1(2) = v_
Call multiplyv(2, 2, r, u1, u2)
Call vadd(2, u2, p1, u3)
ActiveSheet.Cells(isol, 1) = u3(1)
ActiveSheet.Cells(isol, 2) = u3(2)
ActiveSheet.Cells(isol, 4) = v_
End If
Next i
End Sub

Here are the user defined subroutines that the above subroutines calls

Public Sub normalize1(x)
Dim xnorm As Double
xnorm = Sqr(x(1) ^ 2 + x(2) ^ 2 + x(3) ^ 2)
For j = 1 To 3
x(j) = x(j) / xnorm
Next j
End Sub
Public Sub transpose(ByVal m As Integer, 
ByVal n As Integer, a, At)
Dim c1(1500, 1500) As Double
' m  no. of rows of a
' n  no. of columns of a
For i = 1 To n
For j = 1 To m
c1(i, j) = a(j, i)
Next j
Next i
For i = 1 To n
For j = 1 To m
At(i, j) = c1(i, j)
Next j
Next i
End Sub
Public Sub multiplyv(n, m, a, b, c)
Dim i, j, k As Integer
Dim c1(200) As Double
For i = 1 To n
c1(i) = 0
For k = 1 To m
c1(i) = c1(i) + a(i, k) * b(k)
Next k
Next i
For i = 1 To n
c(i) = c1(i)
Next i
End Sub
Public Sub vadd(ByVal n As Integer, x, y, z)
For i = 1 To n
z(i) = x(i) + y(i)
Next i
End Sub
Public Sub vdiff(ByVal n As Integer, x, y, z)
For i = 1 To n
z(i) = x(i) - y(i)
Next i
End Sub
Public Sub vscale(ByVal n, x, a, y)
For i = 1 To n
y(i) = a * x(i)
Next i
End Sub
Public Sub find_quadratic_roots(a2, a1, a0, r)
If Abs(a2) > 0.0000001 Then
xdisc = a1 ^ 2 - 4 * a2 * a0
If xdisc < 0 And Abs(xdisc) < 0.00000001 Then xdisc = 0
If xdisc < 0 Then
r(0) = 0
Else
r(0) = 2
  r(1) = (-a1 - Sqr(xdisc)) / (2 * a2)
  r(2) = (-a1 + Sqr(xdisc)) / (2 * a2)
End If
Else
If Abs(a1) > 0.0000001 Then
  r(0) = 1
  r(1) = -a0 / a1
Else
MsgBox ("Invalid Input in  Find_Quadratic_Roots")
  Exit Sub
End If
End If
End Sub

The output of the program is the following, where column $A$ is for the $x$ coordinate of the center and the $B$ column is for the $y$ coordinate. Column $D$ lists the corresponding radius of the solution circle. There were two valid solutions.

Thanks that's amazing :) Did you wrote sub functions for the normalize1, transpose, multiplyv, vscale, ... your self in Excel or is this some kind of add-in? — Sander, Mar 08 '23 at 13:28
No, it is not an add-in. I wrote all these subroutines in VBA. — Hosam Hajeer, Mar 08 '23 at 13:35
Could you provide those subroutines, so I can see what is happening inside? I'll try to make a C# version and post it in my question above. — Sander, Mar 08 '23 at 15:17

Calculate the inscribed circle between one line and two circles

2 Answers2