centroid.rs

type Point = (f64, f64);

const POLY: &[Point] = &[
    (40547359.805348977, 20991407.4519994),
    (40547406.539216213, 20991408.345532853),
    (40547405.607148655, 20991456.331498079),
    (40547358.891922772, 20991455.404874228),
];

fn naive_centroid(poly: &[Point]) -> Point {
    let mut area = 0.0;
    let mut x = 0.0;
    let mut y = 0.0;
    
    for i in 0..poly.len() {
        let a = &poly[i];
        let b = &poly[(i + 1) % poly.len()];

        let cross = a.0 * b.1 - b.0 * a.1;

        x += (a.0 + b.0) * cross;
        y += (a.1 + b.1) * cross;
        area += 3.0 * cross;
    }

    (x / area, y / area)
}

fn move_to_origin_centroid(poly: &[Point]) -> Point {
    let (x0, y0) = poly[0];
    let moved_to_origin = poly.iter().map(|(x, y)| (x - x0, y - y0)).collect::<Vec<_>>();
    println!("{:?}", moved_to_origin);

    let (x, y) = naive_centroid(&moved_to_origin);
    (x + x0, y + y0)
}

fn stable_centroid(poly: &[Point]) -> Point {
    let mut area = 0.0;
    let mut x = 0.0;
    let mut y = 0.0;
    
    for i in 0..poly.len() {
        let a = &poly[i];
        let b = &poly[(i + 1) % poly.len()];
        
        let w = b.0 * a.1;
        let e = b.0.mul_add(-a.1, w);
        let f = a.0.mul_add(b.1, -w);
        let cross = f + e;

        x += (a.0 + b.0) * cross;
        y += (a.1 + b.1) * cross;
        area += 3.0 * cross;
    }

    (x / area, y / area)
}

fn main() {
    for (idx, p) in POLY.iter().enumerate() {
        println!("POINT #{} = {}, {}", idx, p.0, p.1);
    }
    println!("===");

    let naive = naive_centroid(POLY);
    println!("NAIVE = {}, {}", naive.0, naive.1);

    let centered = move_to_origin_centroid(POLY);
    println!("CENTERED = {}, {}", centered.0, centered.1);

    let stable = stable_centroid(POLY);
    println!("STABLE = {}, {}", stable.0, stable.1);
}