Calculate the inscribed circle between one line and two circles

geometry

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 > 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) && 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($"Solution 1   C_x: {C_x:0.0##}   y: {C_y:0.0##}   radius: {C_radius:0.0##}");
                // 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) && 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($"Solution 2   C_x: {C_x:0.0##}   y: {C_y:0.0##}   radius: {C_radius:0.0##}");
                // Solution 2   x: -48.359   y: 336.121   radius: 329.963
                return true;
            }
        }

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

    }

Best Answer

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$.

Related Solutions

[Math] Cylinder in 3D from five points

Interesting question!

(Is this called a generic conic equation?)

To me, a generic conic would be something planar. You essentially describe the 3d analogon here, so I'd call this a quadric.

I know it would take nine sets of points $[x,y,z]=p_i$ to find the solution to the nine coefficients

And yet you only use 8 variables $a\ldots h$. Well, as you found out and wrote in a comment, that's because you forgot the linear term in $z$. So yes, 9 points would match 9 variables. You assume (without explicitely stating that assumption) a constant coefficient of $1$ for the $x^2$ term. That's permissible except if that coefficient would turn out to be zero. Even more general, you'd take a tenth variable for the $x^2$ term, then use 9 points to determine said variables up to a common scalar factor which is irrelevant as the right hand side is zero. The equation is homogeneous.

Anyone got some thoughts on how to proceed? Or know of a simpler solution?

Abbreviating $x_i:=p_i \cdot u_{0x}$ and $y_i:=p_i \cdot u_{0y}$ as the projected in-plane coordinates of the point, I'd write the circle as

$$(x_i^2+y_i^2)+ax_i+by_i+c=0$$

That way you have 6 unknowns: $u_x,u_y,s,a,b,c$. This description feels easier than your parametrization using $c_x,c_y,r$. But the relations between these are still non-linear, which makes the whole problem quite tricky. You probably should aim to make all conditions polynomial, i.e. rewrite your condition for $s$ to

$$s^2(1+u_y^2)=1+u_x^2\cdot u_y^2+u_x^2$$

You can of course try to let some computer algebra system solve this for you. Or you can look in to elimination theory, and start eliminating some of the unknowns from your equations using Gröbner bases or resultants. Which is probably what your computer algebra system would do internally.

I've actually tried to do the CAS computation with a bunch of random points. Wolfram Alpha doesn't accept as much input as I'd need for this, so I'll be using Sage, and in particular its method for computing the variety of an ideal.

At first I had made a mistake when I tried to use your $u_{0x}$ and $u_{0y}$, see later section for details. So looking for an alternate solution I defined $u_{0x}=(x_1,y_1,0)$ and $u_{0y}=(x_2,y_2,z_2)$ and then added conditions for $u_{0x}\cdot u_{0x}=1$, $u_{0y}\cdot u_{0y}=1$ and $u_{0x}\cdot u_{0y}=0$. Or in other words, both vectors should be unit length and orthogonal to one another. I've assumed that the $z$ coordinate of the axis is non-zero, while you chose $x\neq 0$, but that shouldn't matter. (If you automate this computation, handling such special cases automatically can be troublesome.)

A specific example

With this formulation I tried to solve a specific instance in Sage. I used five points, chosen at random from $\{0,1,2,3,4\}^3$.

$$ p_1=(1, 4, 4)\quad p_2=(4, 3, 4)\quad p_3=(2, 0, 1)\quad p_4=(0, 0, 4)\quad p_5=(2, 4, 0) $$

For these I got $24$ distinct solutions. Actually I usually got $24$ solutions over $\mathbb C$ but only some multiple of $8$ of them real, so the above example was specifically chosen from a number of random examples as one where all the solutions are real.

Here is the Sage code to compute these solutions:

# pts = random_matrix(ZZ, 5, 3, x=0, y=5).rows() # for new random instance
pts = matrix([(1,4,4), (4,3,4), (2,0,1), (0,0,4), (2,4,0)]).rows() # use my case
PR.<a,b,c,x1,y1,z1,x2,y2,z2> = QQ[] # multi-variate polynomial ring
u1 = vector((x1,y1,z1))
u2 = vector((x2,y2,z2))
id1 = ideal([(p*u1)^2 + (p*u2)^2 + a*p*u1 + b*p*u2 + c for p in pts]
          + [z1-0, u1*u2, u1*u1-1, u2*u2-1])
id1.dimension() # must be zero
len(id1.variety(CC)) # will be 24 for most examples
len(id1.variety(RR)) # will be 24 for my example, and usually a multiple of 8

These solutions have to occur in groups of four, since one can negate all coordinates of $u_{0x}$ and/or $u_{0y}$ and still describe essentially the same conic. But this still leaves $24/4=6$ really different cylinders for the situation above. Their axis directions according to my computation using the example above are approximately

$$ (0.678, 0.729, 0.0989) \quad (-0.605, 0.489, -0.628) \quad (-0.408, 0.179, -0.895) \\ (-0.422, -0.104, 0.900) \quad (-0.0774, -0.964, 0.253) \quad (-0.963, -0.108, 0.245) $$

These six coordinates also correspond to six different values, e.g. for the unknown $c$. These six values are the six roots of the polynomial

\begin{align*} 16947006302936317345&\cdot c^6\\ - 423964877845474796312&\cdot c^5\\ + 3414707935712257584272&\cdot c^4\\ - 8163333268350816357120&\cdot c^3\\ - 9034354158004982553600&\cdot c^2\\ + 35770219295155494912000&\cdot c\\ + 13673321316256573440000&= 0 \end{align*}

This I found by exact computation using algebraic numbers (QQbar in Sage, but note that the computation can take considerable time). The Galois group of this polynomial is $S_6$, which is not solvable. So according to Galois theory there can be no way to express the value of $c$ for any of these solutions using radicals. Which in turn means that I'd strongly advise against solving this kind of problem without help from a computer to do the heavy number crunching.

Back to your set of variables

When I first used your vhoice of $u_{0x}$ and $u_{0y}$, I must have forgotten some extra variables, which made me believe that there were additional real degrees of freedom which I couldn't solve for. Once I realized and fixed that mistake, I could apply the above approach to your parametrization as well. I found $12$ solutions, instead of my $24$, which is to be expected as you fixed one coordinate to $1$ which fixed the sign of that vector. So my choice of $u_{0x}$ and $u_{0y}$ is probably not better than yours, while the parametrization of the circle using just three linear unknown variables should still be a good idea.

Feeding my example into your description, the variables that define your vectors are characterized by these polynomials:

\begin{align*} \small 1478610113817600&\cdot u_x^6 & \small 266500&\cdot u_y^6 & \small 2186287868683696026446069760000&\cdot s^{12} \\ \small - 3863371808155392&\cdot u_x^5 & \small + 539357&\cdot u_y^5 & \small -33034531469518995706122669195264&\cdot s^{10} \\ \small - 3541905515941632&\cdot u_x^4 & \small - 8742480&\cdot u_y^4 & \small +96487595035479832160430206853120&\cdot s^8 \\ \small + 565191711017856&\cdot u_x^3 & \small + 12151234&\cdot u_y^3 & \small -106608759197577962918341435697664&\cdot s^6 \\ \small + 691347653169852&\cdot u_x^2 & \small - 5732364&\cdot u_y^2 & \small +44504844774915504691341808859664&\cdot s^4 \\ \small + 59575310863341&\cdot u_x & \small + 1024565&\cdot u_y & \small -3621892801242489875603841840033&\cdot s^2 \\ \small - 1445045015375 &= 0 & \small - 59000 &= 0 & \small + 75225907399919929034745390625 &= 0 \end{align*}

So you “only” need to compute the roots and check which combinations fit together, then you know how to project the points and continue from there. Perhaps this is of use at some point, if only to illustrate how hard an exact solution can be to describe even for very simple input.

[Math] inscribed circles inside circle

"Inscribed" is the key word. You necessarily have despite $d_1$ and $d_2$ can vary a lot $$D=d_1+d_2$$ where $D,d_1,d_2$ are the diameters. Thus the equivalent equality $$D\pi=d_1\pi+d_2\pi$$ Consequently the sums are equal.