Update: Here and here I've found an explanation of the problem
As far as my understanding goes, by finding the intersection of two lines and calculating the angle, then dividing it by two, a line can be obtained which passes through the center of the circle.
As an example I have 2 lines, and a circle:
Line 1
x1 = 7
y1 = 10
x2 = 11
y2 = 6
Line 2
x3 = 12
y3 = 3
x4 = 5
y4 = 0
Circle
rad = 2
xc = 2
yc = 8
And I want to calculate the centerpoint and radius of the circle touching those (tangent).
The result should be something like:
radius = 3.754
xcenter = 6.825
ycenter = 4.866
The final C# function made out of the accepted answer
public static bool GetCircleFromLineLineCircle()
{
//The given values
double P1_x = 7;
double P1_y = 10;
double P2_x = 11;
double P2_y = 6;
double P3_x = 12;
double P3_y = 3;
double P4_x = 5;
double P4_y = 0;
double C_x = 2;
double C_y = 8;
double C_radius = 2;
//1. Horizontal and vertical distance between line points
double d1_x = P2_x - P1_x;
double d1_y = P2_y - P1_y;
double d2_x = P4_x - P3_x;
double d2_y = P4_y - P3_y;
//2. Find the intersection of the two lines containing the two line segments.
double denom = d1_y * d2_x - d1_x * d2_y;
//Return false if lines are parallel or coincident
if (denom == 0) { return false; }
//Intersection point
double T = ((P1_x - P3_x) * d2_y + (P3_y - P1_y) * d2_x) / denom;
double I_x = P1_x + d1_x * T;
double I_y = P1_y + d1_y * T;
//3. Calculate the vector extending from intersection to center
double V_x = C_x - I_x;
double V_y = C_y - I_y;
//4. Calculate the dot product of vector with d1 and d2:
double dotp1 = V_x * d1_x + V_y * d1_y;
double dotp2 = V_x * d2_x + V_y * d2_y;
//If sign is negative, negative line points
if (Math.Sign(dotp1) == -1) { d1_x = -d1_x; d1_y = -d1_y; }
if (Math.Sign(dotp2) == -1) { d2_x = -d2_x; d2_y = -d2_y; }
//5. Calculate the direction vector of the bisector
double dist1 = Math.Sqrt(d1_x * d1_x + d1_y * d1_y);
double dist2 = Math.Sqrt(d2_x * d2_x + d2_y * d2_y);
double U_x = d1_x / dist1 + d2_x / dist2;
double U_y = d1_y / dist1 + d2_y / dist2;
//6. A point on the bisector
double dotp_U_d1 = U_x * d1_x + U_y * d1_y;
double dotp_U_U = U_x * U_x + U_y * U_y;
double dotp_d1_d1 = d1_x * d1_x + d1_y * d1_y;
double theta = Math.Acos(dotp_U_d1 / (Math.Sqrt(dotp_U_U) * Math.Sqrt(dotp_d1_d1)));
//7. The equation we want to solve
double dotp_V_U = V_x * U_x + V_y * U_y;
double dotp_V_V = V_x * V_x + V_y * V_y;
double cos2_theta = Math.Pow(Math.Cos(theta), 2);
double a = dotp_U_U * cos2_theta;
double b = -2 * dotp_V_U - 2 * C_radius * Math.Sqrt(dotp_U_U) * Math.Sin(theta);
double c = dotp_V_V - C_radius * C_radius;
double discriminant = b * b - 4 * a * c;
double t1 = (-b - Math.Sqrt(discriminant)) / (2 * a);
double t2 = (-b + Math.Sqrt(discriminant)) / (2 * a);
//8. Calculate the center and C_radius
//Solution 1
double C1_x = I_x + t1 * U_x;
double C1_y = I_y + t1 * U_y;
double C1_radius = t1 * Math.Abs(Math.Sqrt(dotp_U_U)) * Math.Sin(theta);
//Solution 2
double C2_x = I_x + t2 * U_x;
double C2_y = I_y + t2 * U_y;
double C2_radius = t2 * Math.Abs(Math.Sqrt(dotp_U_U)) * Math.Sin(theta);
Console.WriteLine($"C1_x: {C1_x:0.0##} C1_y: {C1_y:0.0##} C1_radius: {C1_radius:0.0##}");
//C1_x: 6.825 C1_y: 4.866 C1_radius: 3.754
return true;
}
Best Answer
As you said in the question, the center of the unknown circle will lie on the angle bisector of the lines that contain the given line segments. But there are two such bisectors.
To know which bisector to use, we have to locate the center of the given circle within the four segments of the plane defined by the two lines. This is easy to do, as all it takes is finding the dot product of the vector from the intersection of the two lines to the center with the direction vectors of the two lines.
So using the given values, we have
$P_1 = (7,10) , P_2 = (11, 6) , P_3 = (12, 3) , P_4= (5,0) $
Given circle : $C = (2, 8) $ , $ r = 2 $
$ d_1 = P_2 - P_1 = (4, -4) $
$ d_2 = P_4 - P_3 = (-7, -3) $
$ I = P_1 + t d_1 = P_3 + s d_2 $
$ (7, 10) + t ( 4, -4) = (12, 3) + s (-7, -3) $
Therefore,
$4 t + 7 s = 5$
$-4 t + 3 s = -7 $
Adding gives $ 10 s = -2 $, so $ s = - \dfrac{1}{5} $
Hence, the intersection points is
$I = P_3 + s d_2 = (12, 3) - \dfrac{1}{5} (-7, -3) = \dfrac{1}{5}(67 , 18 ) $
$ V = IC = C - I = (2, 8) - \dfrac{1}{5} (67, 18) = \dfrac{1}{5} (-57, 22) $
$ V \cdot d_1 = \dfrac{1}{5} (-57, 22) \cdot (4, -4) \lt 0 $
$ V \cdot d_2 = \dfrac{1}{5}(-57, 22) \cdot (-7, -3) \gt 0 $
Using this, our new direction vectors will be
$ d_1 = (-4, 4) , d_2 = (-7, -3) $
$ U = \dfrac{d_1}{\| d_1 \| } + \dfrac{d_2} { \| d_2 \| } $
This gives
$ U = \dfrac{ (-4, 4) }{ 4 \sqrt{2} } + \dfrac{ (-7, -3) } { \sqrt{ 58 } } $
Numerically, this evaluates to
$ U = (-1.6262518, 0.3131875 ) $
$ P = I + t U $ where $ t \gt 0 $
We want the distance $(\overline{PC} - r )$ to be equal to the perpendicular distance between $P$ and the two lines, i.e. we want
$ (\| PC \| - r) = \| I + t U - C \| = \| PI \| \sin \theta $
where $\theta $ is the angle between the direction of the bisector $U$ and either $d_1$ or $d_2$, so
$ \theta = \cos^{-1} \left(\dfrac{ U \cdot d_1 }{\sqrt{U \cdot U} \sqrt{d_1 \cdot d_1}} \right)$
Numerically, this evaluates to
$ \theta = \cos^{-1} \left(\dfrac{7.7577572} { (1.65613445)(4 \sqrt{2}) }\right) = \cos^{-1}( 0.828067237) = 0.595145 \text{ radians} $
$ \| PC \| - r = \| I + t U - C \| - r = \| PI \| \sin \theta $
Re-arranging,
$ \| I - C + t U \| = t \| U \| \sin \theta + r $
Squaring, we obtain
$ (I - C) \cdot (I - C) + 2 t (I - C) \cdot U + t^2 ( U \cdot U ) = t^2 \| U \|^2 \sin^2 \theta + 2 r t \| U \| \sin \theta + r^2 $
Recall that $ V = C - I $
Therefore the above equation reduces to
$ V \cdot V + t (-2 V \cdot U - 2 r \|U \| \sin \theta ) + t^2 (U \cdot U) \cos^2 \theta - r^2 = 0 $
Numerically, the above equation evaluates to
$ 145.32 - 43.5485 t + 1.8807123 t^2 = 0 $
The roots of this quadratic equation in $t$ are
$ t_1 = 4.04283 $
$ t_2 = 19.1125 $
Since both $t$'s are positive, they're both valid solutions.
Using the first value of $t$, we can obtain the center of the radius of the unknown circle
The center is: $ P = I + t U = \dfrac{1}{5}(67 , 18 ) + 4.04283 (-1.6262518, 0.3131875 ) = (6.82534, 4.8661) $
And the radius given by $ t \|U\| \sin \theta = 4.04283 (1.65613445 ) \sin 0.595145 = 3.75367 $