turing_morphogenesis.rs

/* This is a naive implementation of Jonathan McCabe's elaboration of
 * Alan Turing's model of morphogenesis, as described here:
 * http://www.jonathanmccabe.com/Cyclic_Symmetric_Multi-Scale_Turing_Patterns.pdf
 */

extern crate rand;

use rand::Rng;
use std::env;
use std::io::Write;
use std::ops::{Index,IndexMut};
use std::fs;
use std::path::Path;

struct Config {
    inner_size: usize,
    outer_size: usize,
    incr_amt: f32,
    decr_amt: f32,
    steps: usize,
    img_size: usize,
}

fn clamp(low: f32, high: f32, num: f32) -> f32 {
    if num < low {
        low
    } else if num > high {
        high
    } else {
        num
    }
}

/* We're going to, for simplicity, assume that images are always
 * square, so an image is a size along one dimension plus a vector
 * of floats. */
struct Image {
    sz: usize,
    px: Vec<f32>,
}

impl Image {
    /* Our constructors are straightforward */
    fn new(sz: usize) -> Image {
        let mut v = Vec::new();
        for _ in 0..(sz*sz) {
            v.push(0.0);
        }
        Image { sz: sz, px: v }
    }

    fn new_rand(sz: usize) -> Image {
        let mut rng = rand::thread_rng();
        let mut v = Vec::new();
        for _ in 0..(sz*sz) {
            v.push(rng.gen());
        }
        Image { sz: sz, px: v }
    }

    /* This probably isn't super useful here, but it would be
     * if we kept the Image representation appropriately
     * private and encapsulated. */
    fn dim(&self) -> usize {
        self.sz
    }

    fn in_bounds(&self, x: usize, y: usize) -> bool {
        x < self.sz && y < self.sz
    }

    /* We're printing as a pgm file, which is straightforward,
     * especially as this ignores errors entirely: */
    fn print(&self, mut f: fs::File) {
        /* A header, which is just the text P2 */
        let _ = writeln!(f, "P2");
        /* The width and height in pixels: */
        let _ = writeln!(f, "{} {}", self.dim(), self.dim());
        /* and the maximum greyscale value. I'm just using 128
         * for no real reason. */
        let _ = writeln!(f, "128");
        /* And then we just convert all the floats to ints in the
         * range [0,128] and we're good. We don't need to use a
         * newlines like we do, but it's nice for debugging. */
        for x in 0..self.dim() {
            for y in 0..self.dim() {
                let _ = writeln!(f, "{} ", (self[(x,y)] * 128.0).floor());
            }
            let _ = writeln!(f, "");
        }
    }
}

/* We also want to index our images using pairs, so we can implement
 * the indexing traits to make our code nicer. */
impl Index<(usize, usize)> for Image {
    type Output = f32;

    fn index<'a>(&'a self, (x, y): (usize, usize)) -> &'a f32 {
        &self.px[x * self.sz + y]
    }
}

impl IndexMut<(usize, usize)> for Image {
    fn index_mut<'a>(&'a mut self, (x, y): (usize, usize)) -> &'a mut f32 {
        &mut self.px[x * self.sz + y]
    }
}

/* The actual step-running is pretty simple: */
fn run_step(old: &Image, new: &mut Image, conf: &Config) {
    /* We loop over every pixel in the image... */
    for x in 0..conf.img_size {
        for y in 0..conf.img_size {
            /* And compute the average value of two neighborhoods: a
             * smaller one... */
            let near = gather_neighbors(old, (x, y), conf.inner_size);
            /* And a larger one. */
            let far = gather_neighbors(old, (x, y), conf.outer_size);
            /* Depending on which one is bigger, we either add or subtract
             * from the existing value of the pixel. */
            let dx = if near < far { conf.incr_amt } else { conf.decr_amt };
            /* And we clamp the range to [0.0,1.0] for good measure. */
            new[(x, y)] = clamp(0.0, 1.0, old[(x, y)] + dx);
        }
    }
}

/* Our neighborhood is based on Manhattan distance. We could futz with this
 * to create different, interesting patterns, too. */
fn gather_neighbors(img: &Image, (x, y): (usize, usize), n: usize) -> f32 {
    /* We keep a running average... */
    let mut amt = 0.0;
    let mut tot = 0.0;
    /* This is for going back and forth between signed and unsigned types.
     * There is almost certainly a better way of doing this and I don't
     * really care. */
    let ns = n as isize;
    for i in -ns..ns {
        for j in -ns..ns {
            let xn = (i + x as isize) as usize;
            let yn = (j + y as isize) as usize;
            if img.in_bounds(xn, yn) {
                amt += img[(xn,yn)];
                tot += 1.0;
            }
        }
    }
    amt / tot
}

fn main() {
    /* We can vary these parameters here if we want. */
    let conf = Config {
        inner_size: 18,
        outer_size: 12,
        incr_amt: 0.05,
        decr_amt: -0.05,
        steps: 100,
        img_size: 256,
    };

    /* Take the filename to write to. */
    let filename = match env::args().nth(1) {
        Some(n) => n,
        None => panic!("Usage: [target]"),
    };
    println!("Printing to {:?}", filename);

    /* Run the above step for some number of times */
    let final_image = {
        let mut old = Image::new_rand(conf.img_size);
        let mut new = Image::new(conf.img_size);
        for _ in 0..conf.steps {
            run_step(&old, &mut new, &conf);
            std::mem::swap(&mut new, &mut old);
        }
        new
    };

    /* ...and print to the specified file */
    match fs::File::create(Path::new(&filename)) {
        Ok(file) => final_image.print(file),
        _ => panic!("Unable to open file."),
    };
}