/// Construct a circle from three tangent points. /// /// Returns None if there is no unique solution. This happens when /// the triangle formed by three points has zero area. /// /// The algorithm goes as follows: /// Construct two segment bisectors, defined parametrically: /// /// Bisector for segment ab: /// mid_ab + t * norm_ab = p, where t \in Real /// /// Bisector for segment ac: /// mid_ac + t * norm_ac = p, where t \in Real /// /// Where mid_xy is a midpoint between x and y, /// norm_xy is a vector normal to (y - x) /// /// To calculate intersection between these lines, let's write down /// two criteria of the intersection point: /// /// 1) the intersection point lies on the bisector of ab /// center = mid_ab + t * norm_ab /// /// 2) the intersection point lies on the bisector of ac, written down in /// an alternative way -- that mid_ac - center is collinear to norm_ac /// (mid_ac - center).dot(norm(norm_ac)) == 0 /// (mid_ac - center).dot(dir_ac) == 0 /// /// Combine these equations by substituting center in (2): /// (mid_ac - (mid_ab + t * norm_ab)).dot(dir_ac) == 0 /// /// This is a linear equation that has the following solution: /// (mid_ac - mid_ab - t * norm_ab).dot(dir_ac) == 0 /// (mid_ac - mid_ab).dot(dir_ac) - t * norm_ab.dot(dir_ac) == 0 /// t * norm_ab.dot(dir_ac) == (mid_ac - mid_ab).dot(dir_ac) /// t == (mid_ac - mid_ab).dot(dir_ac) / norm_ab.dot(dir_ac) /// /// The only thing that needs to be checked is the denominator, /// norm_ab.dot(dir_ac).abs() > \epsilon /// /// If it's small, then this is a case of "zero" area. fn circle_from_three(a: Point, b: Point, c: Point) -> Option { const EPS: f32 = 1e-6; let dir_ab = b - a; let dir_ac = c - a; // A vector normal to (x, y) is (-y, x) let norm_ab = Point::new(-dir_ab.y, dir_ab.x); // The denominator for the solution to the linear equation let dir_sin = norm_ab.dot(dir_ac); // Triangle area is too small, there's no unique solution if dir_sin.abs() < EPS { return None; } // Calculate midpoints for ab and ac let mid_ab = (a + b) * 0.5; let mid_ac = (a + c) * 0.5; // Solution to the linear equation let t = (mid_ac - mid_ab).dot(dir_ac) / dir_sin; // Use the t in parametric definition for bisector of ab let center = mid_ab + t * norm_ab; // Radius is equal to distance to any of three points let radius = (a - center).len(); Some(Circle { center, radius }) }