Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
/// 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 })
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment