Cargo.toml

[package]
name = "planet"
version = "0.0.1"

[dependencies]
sdl2 = "0.19"
sdl2_image = "0.19"
cgmath = "*"
rand = "*"
rayon = "*"

planet.rs

extern crate sdl2;
extern crate sdl2_image;
extern crate cgmath;
extern crate rand;
extern crate rayon;

use rand::{XorShiftRng, Rand, Rng, SeedableRng};
use sdl2::surface::{Surface, SurfaceRef};
use sdl2::pixels::PixelFormatEnum::ARGB8888;
use sdl2_image::SaveSurface;

type Vec2u = cgmath::Vector2<u32>;
type Vec2f = cgmath::Vector2<f64>;
type Vec3f = cgmath::Vector3<f64>;

#[derive(Copy, Clone, Debug)]
struct Altitude(f64);

#[derive(Copy, Clone, Debug)]
struct Seed([u32; 4]);

#[derive(Copy, Clone, Debug)]
struct Vertex {
    alt: Altitude,
    pos: Vec3f,
    seed: Seed,
}

#[derive(Copy, Clone, Debug)]
struct Tetra(Vertex, Vertex, Vertex, Vertex);

#[derive(Copy, Clone, Debug)]
struct Plane(Vec3f, Vec3f, Vec3f);

#[derive(Clone, Debug)]
struct Image {
    width: u32,
    height: u32,
    pixels: Vec<Option<Altitude>>,
    altmax: Altitude,
}

#[derive(Copy, Clone, Debug)]
struct ImageSize(u32, u32);

#[derive(Copy, Clone, Debug)]
struct Rgb(u8, u8, u8);

#[derive(Clone, Debug)]
struct ColorTable(Vec<(u16, Rgb)>);

#[derive(Copy, Clone, Debug)]
enum Projection {
    Ortho,
    Square,
}

fn main() {
    //let seed = [111111, 413, 612, 1025];
    let seed = Seed([11123249, 849348, 279347923, 849833048]);
    let default_altitude = Altitude(-0.01);

    let t = Tetra::new(seed, default_altitude);
    let img_size = ImageSize::default();
    let ImageSize(width, height) = img_size;

    let scale = 0.9;

    let sdl = sdl2::init().unwrap();
    let vid = sdl.video().unwrap();
    let mut window = vid.window("Planet", width, height)
        .position_centered()
        .build()
        .unwrap();

    let mut event_pump = sdl.event_pump().unwrap();

    let mut angle = 0.0;

    let colors = ColorTable::default();

    'render: loop {
        for event in event_pump.poll_iter() {
            use sdl2::event::Event;

            match event {
                Event::Quit { .. } => break 'render,
                Event::KeyDown { keycode: Some(key), .. } => {
                    use sdl2::keyboard::Keycode;

                    match key {
                        Keycode::Q => break 'render,
                        Keycode::Escape => break 'render,
                        _ => (),
                    }
                }
                _ => (),
            }
        }

        let img = Image::render(img_size, scale, angle, t);

        match window.surface_mut(&mut event_pump) {
            Ok(surf) => img.draw_to_surface(surf, &colors),
            Err(e) => { println!("{}", e); break 'render },
        }

        window.update_surface().unwrap();
        angle += 10.0;
    }

    let _img = { use sdl2_image::*; init(INIT_PNG) }.unwrap();
    let img = Image::render(img_size, scale, angle, t);
    img.write("test.png", &colors).unwrap();
}

fn sample(t: Tetra, p: Vec3f, l: f64, depth: u8) -> Altitude {
    use cgmath::MetricSpace;

    let Tetra(a, b, c, d) = t.sort();
    let dist = a.pos.distance(b.pos);
    if dist < l {
        return t.mean_altitude();
    }

    let m = {
        let seed = random(Seed::mean(a.seed, b.seed));

        Vertex {
            seed: seed,
            pos: (a.pos + b.pos) / 2.0,
            alt: Altitude::offset(seed, dist, a.alt, b.alt),
        }
    };

    // The paper calls for checking whether p is inside Tetra(m, b, c, d),
    // but the only alternative is that it's inside Tetra(a, m, c, d), so
    // at this step we only check the bisecting plane.
    if Plane(m.pos, c.pos, d.pos).contains_both(b.pos, p) {
        sample(Tetra(m, b, c, d), p, l, depth + 1)
    } else {
        sample(Tetra(a, m, c, d), p, l, depth + 1)
    }
}

fn random(seed: Seed) -> Seed {
    seed.into_rng().gen()
}

#[test]
fn random_should_alter_seed() {
    let old_seed = Seed([123249, 849348, 279347923, 849833048]);
    let new_seed = random(old_seed);
    assert!(old_seed.0 != new_seed.0);
}

impl Image {
    fn render(size: ImageSize, scale: f64, angle: f64, t: Tetra) -> Self {
        let ImageSize(width, height) = size;

        //let square_projection = |x, y| {
        //};

        let projection = |x, y| {
            let half_w = width as f64 / 2.0;
            let wx = (x as f64 - half_w) / (half_w * scale);
            let half_h = height as f64 / 2.0;
            let wy = (half_h - y as f64) / (half_h * scale);
            let r2 = wx*wx + wy*wy;
            if r2 > 1.0 { return None; }
            let wz = (1.0 - r2).sqrt();
            Some(Vec3f::new(wx, wy, wz))
        };

        let l = 1.0 / (height as f64 * 8.0);

        use cgmath::{Rotation, Rotation3};

        let angle = cgmath::Rad((-angle).to_radians());
        let rot = cgmath::Basis3::from_angle_y(angle);

        use rayon::prelude::*;

        let mut pixels = vec![None; (width * height) as usize];
        pixels.par_chunks_mut(width as usize)
            .enumerate()
            .for_each(|(y, row)| {
                for x in 0 .. width {
                    row[x as usize] = projection(x, y).map(|v| {
                        let v = rot.rotate_vector(v);
                        sample(t, v, l, 0)
                    })
                }
            });

        let mut altmax = 0.0;
        for pixel in pixels.iter() {
            if let &Some(Altitude(ref a)) = pixel {
                let mag = a.abs();
                if mag > altmax {
                    altmax = mag;
                }
            }
        }

        Image {
            width: width,
            height: height,
            pixels: pixels,
            altmax: Altitude(altmax),
        }
    }

    fn draw_to_surface(&self, surf: &mut SurfaceRef, colors: &ColorTable) {
        surf.with_lock_mut(|data| {
            for (i, p) in self.pixels.iter().cloned().enumerate() {
                let alt = match p {
                    Some(Altitude(a)) => a,
                    None => continue,
                };

                let (r, g, b, a) = match colors.lookup(alt, self.altmax) {
                    Some(Rgb(r, g, b)) => (r, g, b, 255),
                    None => (0, 0, 0, 0),
                };

                let j = i * 4;

                // SDL seems to have some funny notions about channel order on
                // little-endian architectures. This might be OS-dependent, so
                // watch your back.
                data[j + 0] = b;
                data[j + 1] = g;
                data[j + 2] = r;
                data[j + 3] = a;
            }
        });
    }

    fn write(self, path: &str, colors: &ColorTable) -> Result<(), String> {
        let mut surf = try!(Surface::new(self.width, self.height, ARGB8888));
        self.draw_to_surface(&mut surf, colors);
        surf.save(path.as_ref())
    }
}

//impl Projection {
//    fn project(self, x: f64, y: f64) -> Option<Vec3f> {
//        match self {
//            Projection::Ortho => {
//            },
//            Projection::Square => {
//            },
//        }
//    }
//}

impl ColorTable {
    fn lookup(&self, a: f64, altmax: Altitude) -> Option<Rgb> {
        let a: u16 = {
            let (low, high) = match self.range() {
                None => return None,
                Some(r) => r,
            };

            let Altitude(altmax) = altmax;
            let span = (high - low) as f64;

            ((a + altmax) * (span / (altmax * 2.0))) as u16
        };

        fn lerp(low: u8, high: u8, amt: f64) -> u8 {
            if low > high { return lerp(high, low, 1.0 - amt); }
            ((high - low) as f64 * amt) as u8 + low
        }

        let mut prev = *self.0.first().unwrap();
        for (i, rgb) in self.0.iter().cloned() {
            if a == i { return Some(rgb) }
            let (prev_i, prev_rgb) = prev;
            if a < i {
                let amt = (a - prev_i) as f64 / (i - prev_i) as f64;
                let r = lerp(prev_rgb.0, rgb.0, amt);
                let g = lerp(prev_rgb.1, rgb.1, amt);
                let b = lerp(prev_rgb.2, rgb.2, amt);
                return Some(Rgb(r, g, b))
            }
            prev = (i, rgb);
        }

        Some(prev.1)
    }

    fn range(&self) -> Option<(u16, u16)> {
        match (self.0.first(), self.0.last()) {
            (Some(low), Some(high)) => Some((low.0, high.0)),
            _ => None,
        }
    }
}

impl Plane {
    fn contains_both(self, a: Vec3f, b: Vec3f) -> bool {
        use cgmath::InnerSpace;
        let Plane(p1, p2, p3) = self;
        let normal = (p2 - p1).cross(p3 - p1);
        let dot_a = normal.dot(a - p1);
        let dot_b = normal.dot(b - p1);
        dot_a.signum() == dot_b.signum()
    }
}

impl Tetra {
    fn new(seed: Seed, initial_altitude: Altitude) -> Self {
        let mut rng = XorShiftRng::from_seed(seed.0);

        let q = 3.0f64.sqrt();

        let a = Vertex {
            seed: rng.gen(),
            alt: initial_altitude,
            pos: Vec3f::new(-q - 0.20, -q - 0.22, -q - 0.23),
        };

        let b = Vertex {
            seed: rng.gen(),
            alt: initial_altitude,
            pos: Vec3f::new(-q - 0.19, q + 0.18, q + 0.17),
        };

        let c = Vertex {
            seed: rng.gen(),
            alt: initial_altitude,
            pos: Vec3f::new(q + 0.21, -q - 0.24, q + 0.15),
        };

        let d = Vertex {
            seed: rng.gen(),
            alt: initial_altitude,
            pos: Vec3f::new(q + 0.24, q + 0.22, -q - 0.25),
        };

        Tetra(a, b, c, d)
    }

    /// Rearrange the vertices so AB is the longest edge.
    fn sort(self) -> Self {
        use cgmath::MetricSpace;
        let Tetra(a, b, c, d) = self;

        enum Edge { AB, AC, AD, BC, BD, CD, }

        // TODO: See if we can skip some of this.

        let edge_lengths = vec![
            (Edge::AC, a.pos.distance2(c.pos)),
            (Edge::AD, a.pos.distance2(d.pos)),
            (Edge::BC, b.pos.distance2(c.pos)),
            (Edge::BD, b.pos.distance2(d.pos)),
            (Edge::CD, c.pos.distance2(d.pos)),
        ];

        let mut max = (Edge::AB, a.pos.distance2(b.pos));
        for (name, len) in edge_lengths.into_iter() {
            if len > max.1 { max = (name, len); }
        }

        match max.0 {
            Edge::AB => self,
            Edge::AC => Tetra(a, c, b, d),
            Edge::AD => Tetra(a, d, c, b),
            Edge::BC => Tetra(c, b, a, d),
            Edge::BD => Tetra(d, b, c, a),
            Edge::CD => Tetra(d, c, b, a),
        }
    }

    fn mean_altitude(self) -> Altitude {
        let Tetra(a, b, c, d) = self;
        Altitude((a.alt.0 + b.alt.0 + c.alt.0 + d.alt.0) / 4.0)
    }
}

impl Altitude {
    /// Calculate a new random altitude using the given seed.
    fn offset(seed: Seed, mut dist: f64, a1: Self, a2: Self) -> Self {
        // I swiped these from the reference implementation after having
        // trouble with rendering artifacts. Ought to be configurable.
        const WEIGHT_ALT: f64 = 0.45;
        const POW_ALT: f64 = 1.0;
        const WEIGHT_DIST: f64 = 0.035;
        const POW_DIST: f64 = 0.47;

        let midpoint = (a1.0 + a2.0) / 2.0;

        let mut rng = seed.into_rng();
        let es = rng.gen::<f64>() - 0.5;
        let es1 = rng.gen::<f64>() - 0.5;

        let alt_mod = es * WEIGHT_ALT * (a1.0 - a2.0).abs().powf(POW_ALT);

        if dist > 1.0 { dist = dist.powf(0.5); }
        let dist_mod = es1 * WEIGHT_DIST * dist.powf(POW_DIST);

        Altitude(midpoint + alt_mod + dist_mod)
    }
}

impl Seed {
    fn mean(lhs: Self, rhs: Self) -> Self {
        let Seed(mut lhs) = lhs;
        let Seed(rhs) = rhs;
        for i in 0 .. 4 {
            lhs[i] = lhs[i].wrapping_add(rhs[i]) / 2;
        }
        Seed(lhs)
    }

    fn into_rng(self) -> XorShiftRng {
        let Seed(s) = self;
        XorShiftRng::from_seed(s)
    }
}

impl Rand for Seed {
    fn rand<R: Rng>(rng: &mut R) -> Self {
        Seed(rng.gen())
    }
}

impl Rand for Altitude {
    fn rand<R: Rng>(rng: &mut R) -> Self {
        let a: f64 = rng.gen();
        let b = a - 0.5;
        Altitude(b * 0.1 - 0.01)
    }
}

impl Default for ImageSize {
    fn default() -> Self {
        ImageSize(512, 512)
    }
}

impl Default for ColorTable {
    fn default() -> Self {
        ColorTable(vec![
            (6, Rgb(55, 28, 75)),
            (130, Rgb(136, 92, 80)),
            (131, Rgb(148, 100, 82)),
            (180, Rgb(232, 164, 90)),
            (230, Rgb(248, 180, 92)),
            (255, Rgb(165, 128, 108)),
        ])

        //ColorTable(vec![
        //    (1, Rgb(0, 0, 0)),
        //    (16384, Rgb(255, 255, 255)),
        //])
    }
}