This is an image of an Apollonian Gasket:

It's basically what you get you pack circles into a circle.

Now some YouTube video hit me with the fact that if the radius of the outer circle is 1, then the radius of the two biggest inner circles is 1/2, the the radius of the next smallest circles is 1/3. and so on. But the curvature of a circle is 1/radius so the curvatures go 1, 2, 3, .... Amazingly every circle in the packing has an exactly integer curvature. WTF? How could that be?

That leads to a lot of googling around and of course culminates in a Rust program to generate such an image.

use num_complex::Complex64;
use std::fs::File;
use std::io::Write;

#[derive(Debug, Clone, Copy)]
struct Circle {
    center: Complex64,
    radius: f64,
    is_enclosing: bool,
    depth: usize,
}

impl Circle {
    fn new(x: f64, y: f64, radius: f64, is_enclosing: bool, depth: usize) -> Self {
        Circle {
            center: Complex64::new(x, y),
            radius,
            is_enclosing,
            depth,
        }
    }

    fn curvature(&self) -> f64 {
        if self.is_enclosing { -1.0 / self.radius } else { 1.0 / self.radius }
    }
}

// Converts HSL colors directly to hex strings so Inkscape can parse them without crashing
fn hsl_to_hex(h: f64, s: f64, l: f64) -> String {
    let s = s / 100.0;
    let l = l / 100.0;

    let c = (1.0 - (2.0 * l - 1.0).abs()) * s;
    let x = c * (1.0 - (((h / 60.0) % 2.0) - 1.0).abs());
    let m = l - c / 2.0;

    let (r_prim, g_prim, b_prim) = if h < 60.0 {
        (c, x, 0.0)
    } else if h < 120.0 {
        (x, c, 0.0)
    } else if h < 180.0 {
        (0.0, c, x)
    } else if h < 240.0 {
        (0.0, x, c)
    } else if h < 300.0 {
        (x, 0.0, c)
    } else {
        (c, 0.0, x)
    };

    let r = ((r_prim + m) * 255.0).round() as u8;
    let g = ((g_prim + m) * 255.0).round() as u8;
    let b = ((b_prim + m) * 255.0).round() as u8;

    format!("#{:02x}{:02x}{:02x}", r, g, b)
}

fn find_fourth_circle_excluding(c1: &Circle, c2: &Circle, c3: &Circle, exclude: Option<&Circle>, next_depth: usize) -> Vec<Circle> {
    let k1 = c1.curvature();
    let k2 = c2.curvature();
    let k3 = c3.curvature();

    let sum_k = k1 + k2 + k3;
    let product_sum_k = k1 * k2 + k2 * k3 + k3 * k1;
    let radical_k = 2.0 * product_sum_k.max(0.0).sqrt();

    let k4_solutions = [sum_k + radical_k, sum_k - radical_k];

    let z1 = c1.center;
    let z2 = c2.center;
    let z3 = c3.center;

    let term1 = k1 * z1 + k2 * z2 + k3 * z3;
    let term2 = k1 * k2 * z1 * z2 + k2 * k3 * z2 * z3 + k3 * k1 * z3 * z1;
    let radical_z = 2.0 * term2.sqrt();

    let mut valid_circles = Vec::new();
    let epsilon = 1e-4;

    for &k4 in &k4_solutions {
        if k4.abs() < epsilon { continue; }

        let r4 = 1.0 / k4.abs();
        let k4_is_enclosing = k4 < 0.0;

        let possible_centers = [
            (term1 + radical_z) / k4,
            (term1 - radical_z) / k4,
        ];

        for z4 in possible_centers {
            let mut is_valid = true;
            for &c in &[c1, c2, c3] {
                let actual_dist = (z4 - c.center).norm();
                let expected_dist = match (c.is_enclosing, k4_is_enclosing) {
                    (true, true)   => (c.radius - r4).abs(),
                    (true, false)  => c.radius - r4,
                    (false, true)  => r4 - c.radius,
                    (false, false) => c.radius + r4,
                };

                if (actual_dist - expected_dist).abs() > epsilon {
                    is_valid = false;
                    break;
                }
            }

            if let Some(ex) = exclude {
                if (z4 - ex.center).norm() < epsilon && (r4 - ex.radius).abs() < epsilon {
                    is_valid = false;
                }
            }

            if is_valid {
                if !valid_circles.iter().any(|c: &Circle| (c.center - z4).norm() < epsilon) {
                    valid_circles.push(Circle { 
                        center: z4, 
                        radius: r4, 
                        is_enclosing: k4_is_enclosing,
                        depth: next_depth 
                    });
                }
            }
        }
    }
    valid_circles
}

fn generate_gasket(c1: Circle, c2: Circle, c3: Circle, depth: usize, max_depth: usize, all_circles: &mut Vec<Circle>) {
    if depth > max_depth { return; }

    let next_circles = find_fourth_circle_excluding(&c1, &c2, &c3, None, depth + 1);

    for c4 in next_circles {
        let duplicate = all_circles.iter().any(|c| (c.center - c4.center).norm() < 1e-4 && (c.radius - c4.radius).abs() < 1e-4);
        if duplicate { continue; }

        all_circles.push(c4);

        generate_gasket(c1, c2, c4, depth + 1, max_depth, all_circles);
        generate_gasket(c1, c3, c4, depth + 1, max_depth, all_circles);
        //generate_gasket(c2, c3, c4, depth + 1, max_depth, all_circles);
    }
}

fn main() -> std::io::Result<()> {
    let c1 = Circle::new(0.0, 0.0, 2.0, true, 0);
    let c2 = Circle::new(-1.0, 0.0, 1.0, false, 0);
    let c3 = Circle::new(1.0, 0.0, 1.0, false, 0);

    let mut all_circles = vec![c1, c2, c3];
    generate_gasket(c1, c2, c3, 0, 8, &mut all_circles);

    println!("Total filled circles generated: {}", all_circles.len());

    let mut file = File::create("gasket.svg")?;
    let size = 600;
    let scale = 140.0;
    let offset = 300.0;

    writeln!(
        file, 
        "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"600px\" height=\"600px\" viewBox=\"0 0 600 600\">"
    )?;
    writeln!(file, "  <rect width=\"600\" height=\"600\" fill=\"#0b0c10\"/>")?;

    for circle in all_circles {
        let cx = circle.center.re * scale + offset;
        let cy = circle.center.im * scale + offset;
        let r = circle.radius * scale;

        if cx.is_nan() || cy.is_nan() || r.is_nan() || r <= 0.0 {
            continue; 
        }

        if circle.is_enclosing {
            writeln!(file, "  <circle cx=\"{:.2}\" cy=\"{:.2}\" r=\"{:.2}\" fill=\"#1f2833\" stroke=\"#45f3ff\" stroke-width=\"0\"/>", cx, cy, r)?;
        } else {
            let hue = (180 + (circle.depth * 35)) % 360;
            let hex_color = hsl_to_hex(hue as f64, 85.0, 50.0);
            
            writeln!(
                file, 
                "  <circle cx=\"{:.3}\" cy=\"{:.3}\" r=\"{:.3}\" fill=\"{}\" fill-opacity=\"0.85\" stroke=\"#000000\" stroke-width=\"0\"/>", 
                cx, cy, r, hex_color
            )?;
        }
    }

    writeln!(file, "</svg>")?;
    println!("Successfully exported clean vector file to 'gasket.svg'!");

    Ok(())
}

What else is a lonely boy supposes to do on a Saturday night?

Now that is great but it takes 6 minutes to generate that image on my MacBook M1. If anyone has any suggestions to turbo charge it that would be great.

1 post - 1 participant

Read full topic