Apollonian Gasket - My Saturday night project

[ZiCog]

This is an image of an Apollonian Gasket:

It's basically what you get when 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.

[danjl1100]

I feel like this would benefit from a collection that's aware of the circle centers as a unique sortable key, so the duplicates check is not a linear search each time in the innermost function.

Not sure if BTreeSet with a suitable (Partial)Ord impl would be enough, but something to try.

Or maybe a HashMap<CircleCenter, CircleOtherData>?

Though I guess Hash/Eq is tricky with floats...

[ZiCog]

Well of course better than detecting duplicates is not creating duplicates in the first place.

Turns out our slow run time problem lies in the fact that as we recursively generate ever smaller circles from circles we are effectively creating big bushy tree, and that many of the nodes of that tree represent circles of the same radius and position as other nodes on of the tree. Only the tree does not know that, it's just crating more nodes on different branches.

This is compounded by the fact that the "fix" to that problem is a duplicate detector that runs in N squared time or some such.

So we end up wasting time creating circles we don't need and then more time trying to detect them.

The revised Apollonian Gasket generator below removes the duplicate detector and adds logic to not create duplicate circles in the first place.

The result is stunning. With a recursion depth of 12 the run time is down from over 6 minutes to 0.13 seconds on my Mac M1. A speed up of over 2500 times!

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
        }
    }
}

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)
}

/// Finds the single correct fourth circle tangent to c1, c2, c3
/// while explicitly excluding the 'exclude' circle.
fn find_fourth_circle(
    c1: &Circle,
    c2: &Circle,
    c3: &Circle,
    exclude: &Circle,
    next_depth: usize,
) -> Option<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();

    // Descartes' Circle Theorem yields two curvature roots
    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 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 {
            // Step 1: Skip if it matches the parent circle we want to avoid backtracking to
            if (z4 - exclude.center).norm() < epsilon && (r4 - exclude.radius).abs() < epsilon {
                continue;
            }

            // Step 2: Verify geometry layout configurations match up safely
            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 is_valid {
                return Some(Circle {
                    center: z4,
                    radius: r4,
                    is_enclosing: k4_is_enclosing,
                    depth: next_depth,
                });
            }
        }
    }
    None
}

/// Structured recursion without global arrays or search lookups
fn generate_gasket_recursive(
    c1: Circle,
    c2: Circle,
    c3: Circle,
    exclude: Circle,
    depth: usize,
    max_depth: usize,
    all_circles: &mut Vec<Circle>,
) {
    if depth > max_depth {
        return;
    }

    if let Some(c4) = find_fourth_circle(&c1, &c2, &c3, &exclude, depth) {
        all_circles.push(c4);

        // Branch out recursively to the three new tangent combinations
        generate_gasket_recursive(c1, c2, c4, c3, depth + 1, max_depth, all_circles);
        generate_gasket_recursive(c1, c3, c4, c2, depth + 1, max_depth, all_circles);
        generate_gasket_recursive(c2, c3, c4, c1, depth + 1, max_depth, all_circles);
    }
}

pub 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];

    // Initial step: find the inner root circles between the starting trio setup
    let k1 = c1.curvature();
    let k2 = c2.curvature();
    let k3 = c3.curvature();
    let sum_k = k1 + k2 + k3;
    let radical_k = 2.0 * (k1 * k2 + k2 * k3 + k3 * k1).max(0.0).sqrt();

    // The two initial solutions within the bounding outer circle frame
    let k4_roots = [sum_k + radical_k, sum_k - radical_k];
    let term1 = k1 * c1.center + k2 * c2.center + k3 * c3.center;
    let radical_z = 2.0
        * (k1 * k2 * c1.center * c2.center
            + k2 * k3 * c2.center * c3.center
            + k3 * k1 * c3.center * c1.center)
            .sqrt();

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

    for &k4 in &k4_roots {
        let r4 = 1.0 / k4.abs();
        let possible_centers = [(term1 + radical_z) / k4, (term1 - radical_z) / k4];
        for z4 in possible_centers {
            let mut valid = true;
            for c in &[c1, c2, c3] {
                let actual = (z4 - c.center).norm();
                let expected = match (c.is_enclosing, k4 < 0.0) {
                    (true, true) => (c.radius - r4).abs(),
                    (true, false) => c.radius - r4,
                    (false, true) => r4 - c.radius,
                    (false, false) => c.radius + r4,
                };
                if (actual - expected).abs() > epsilon {
                    valid = false;
                }
            }
            if valid
                && !initial_inner_circles
                    .iter()
                    .any(|c: &Circle| (c.center - z4).norm() < epsilon)
            {
                initial_inner_circles.push(Circle::new(z4.re, z4.im, r4, k4 < 0.0, 1));
            }
        }
    }

    // Begin generation across all primary open gaps
    for c4 in initial_inner_circles {
        all_circles.push(c4);
        generate_gasket_recursive(c1, c2, c4, c3, 2, 11, &mut all_circles);
        generate_gasket_recursive(c1, c3, c4, c2, 2, 11, &mut all_circles);
        generate_gasket_recursive(c2, c3, c4, c1, 2, 11, &mut all_circles);
    }

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

    let mut file = File::create("gasket.svg")?;
    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 r < 0.25 {
            continue;
        }

        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=\"{:.2}\" cy=\"{:.2}\" r=\"{:.2}\" 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(())
}

[StefanSalewski]

@ZiCog

The result is stunning.

Was an AI able to find that improvement?

Actually, when you posted the original code some days ago, I already considered asking an AI for improvements, but finally I forgot about it

I think such a task is not a bad test, at least more significant than the popular Rubic cube simulator often used by a few YouTubers.

But of course the solution might be just contained in the training data as well.

[MOCKBA]

The next project should be -Apollo finds Venus.

[ZiCog]

@MOCKBA

Apollo finds Venus

What's that? I don't get the reference.