circle_from_tangent_points.rs

/// 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<Circle> {
    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 })
}