Rug Cellular Automata (demo)

rug.cpp

///// Rug Cellular Automata (demo)
///// Copyright (c) 2013, Sergey Shishmintzev
///// License: Apache/GPLv3
///// Idea: http://www.fourmilab.ch/cellab/manual/rules.html#Rug
///// Example: http://youtu.be/OQQ8zq-xjXI
///// Compiler: gcc-4.7.3 (cygwin)
///// Compile command: g++ -O3  rug.cpp -lpng
///// P.S. Sorry for my English. It isn't my mother tongue.

#include <assert.h>
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <png.h>
#include <iostream>
#include <string>
#include <vector>
#include <boost/numeric/ublas/matrix.hpp>

using namespace boost::numeric::ublas;

inline int mod(int a, int b) { // http://stackoverflow.com/questions/4003232/how-to-code-a-modulo-operator-in-c-c-obj-c-that-handles-negative-numbers
	if(b < 0) //you can check for b == 0 separately and do what you want
		return mod(-a, -b);   
	int ret = a % b;
	if (ret < 0)
		ret += b;
	return ret;
}


inline int make_rgb(int r, int g, int b) {
	return (b&0xFF) | ((g&0xFF)<<8) | ((r&0xFF)<<16);
}


inline int get_r(int rgb) {
	return (rgb>>16) & 0xFF;
}


inline int get_g(int rgb) {
	return (rgb>>8) & 0xFF;
}


inline int get_b(int rgb) {
	return (rgb) & 0xFF;
}


inline void setRGB(png_byte *ptr, int val) {
	ptr[0] = get_r(val);
	ptr[1] = get_g(val);
	ptr[2] = get_b(val);
}


int save_matrix(const matrix<int> &m, const std::string &filename, const std::vector<int> &palette) {
	// http://www.labbookpages.co.uk/software/imgProc/libPNG.html
	int code = 0;
	FILE *fp = NULL;
	png_structp png_ptr = NULL;
	png_infop info_ptr = NULL;
	png_bytep row = NULL;
	
	// Open file for writing (binary mode)
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
		fprintf(stderr, "Could not open file %s for writing\n", filename.c_str());
		code = 1;
		goto finalise;
	}

	// Initialize write structure
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
		fprintf(stderr, "Could not allocate write struct\n");
		code = 1;
		goto finalise;
	}

	// Initialize info structure
	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
		fprintf(stderr, "Could not allocate info struct\n");
		code = 1;
		goto finalise;
	}

	// Setup Exception handling
	if (setjmp(png_jmpbuf(png_ptr))) {
		fprintf(stderr, "Error during png creation\n");
		code = 1;
		goto finalise;
	}

	png_init_io(png_ptr, fp);

	// Write header (8 bit colour depth)
	png_set_IHDR(png_ptr, info_ptr, m.size2(), m.size1(),
			8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
			PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

	png_write_info(png_ptr, info_ptr);

	// Allocate memory for one row (3 bytes per pixel - RGB)
	row = (png_bytep) malloc(3 * m.size2() * sizeof(png_byte));

	// Write image data
	unsigned int i, j;
	for (i = 0; i < m.size1(); i++) {
		for (j = 0; j < m.size2(); j++) {
			setRGB(&(row[j*3]), palette[m(i, j)]);
		}
		png_write_row(png_ptr, row);
	}

	// End write
	png_write_end(png_ptr, NULL);

finalise:
	if (fp != NULL) fclose(fp);
	if (info_ptr != NULL) png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
	if (png_ptr != NULL) png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
	if (row != NULL) free(row);

	return code;
}

// Gradient palette routines
// http://jtauber.com/blog/2008/05/18/creating_gradients_programmatically_in_python/
int linear_gradient(int start_value, int stop_value, float start_offset, float stop_offset, float offset) {
	return floor(start_value + ((offset - start_offset) / (stop_offset - start_offset) * (stop_value - start_value)));
}

int gradient(int start, int end, float initial_offset, float stop_offset, float offset) {
	// make gradient color	
	int r = linear_gradient(get_r(start), get_r(end), initial_offset, stop_offset, offset);
	int g = linear_gradient(get_g(start), get_g(end), initial_offset, stop_offset, offset);
	int b = linear_gradient(get_b(start), get_b(end), initial_offset, stop_offset, offset);
	return make_rgb(r, g, b);
}


int main() {
	using std::cout;
	using std::endl;
	using std::min;
	using std::max;
	using std::vector;

	// random generators seed
	const int modulo = pow(2, 11); // arbitrary power of two number
	// output image size
	const int width = 1920;
	const int height = 1080;

	cout << "modulo = " << modulo << endl;
	cout << "size = " << width << "x" << height << endl;

	// generate multiband palette
	vector<int> palette(modulo);

	const int color1 = make_rgb(250, 100, 30); 
	const int color2 = make_rgb(99, 222, 22);
	const int color3 = make_rgb(100, 10, 233);
	const int color4 = make_rgb(190, 50, 11);
	const int color5 = make_rgb(211, 200, 1);
	const int color6 = make_rgb(111, 20, 144);

	for(int i = 0; i < modulo/2; i++)
		palette[i] = gradient(color1, color2, 0.0, 1.0, i*1.0/(modulo/2));
	for(int i = 0; i < modulo/4; i++)
		palette[modulo/2+i] = gradient(color2, color3, 0.0, 1.0, i*1.0/(modulo/4));
	for(int i = 0; i < modulo/8; i++)
		palette[modulo/2+modulo/4+i] = gradient(color3, color4, 0.0, 1.0, i*1.0/(modulo/8));
	for(int i = 0; i < modulo/16; i++)
		palette[modulo/2+modulo/4+modulo/8+i] = gradient(color4, color5, 0.0, 1.0, i*1.0/(modulo/16));
	for(int i = 0; i < modulo/16; i++)
		palette[modulo/2+modulo/4+modulo/8+modulo/16+i] = gradient(color5, color6, 0.0, 1.0, i*1.0/(modulo/16));
	

	const int frames = 50000; // count of the images to produce

	boost::array<matrix<int>, 2> _state;
	_state.assign(matrix<int>(height, width));
	int cur_state_index = 0;
	int next_state_index = 1;

	// seed
	_state[0](height/2, width/4 + width/8) = modulo-1;
	_state[0](height/2, width/2 + width/4 - width/8) = modulo-1;

	for (int frame = 0; frame < frames; frame++) {
		cout << "frame #" << frame << endl;
		matrix<int> &next_state(_state[next_state_index]);
		matrix<int> &cur_state(_state[cur_state_index]);
		next_state.clear();

		int i0, i1, i2, j0, j1, j2;
		int s00, s01, s02, s10, s11, s12, s20, s21, s22;
		int val;
		for (i1 = 0; i1 < height; i1++) {
			i0 = mod(i1 - 1, height);
			i2 = mod(i1 + 1, height);
			for (j1 = 0; j1 < width; j1++) {
				j0 = mod(j1 - 1, width);
				j2 = mod(j1 + 1, width);
				s00 = cur_state(i0, j0);
				s01 = cur_state(i0, j1);
				s02 = cur_state(i0, j2);
				s10 = cur_state(i1, j0);
				s11 = cur_state(i1, j1);
				s12 = cur_state(i1, j2);
				s20 = cur_state(i2, j0);
				s21 = cur_state(i2, j1);
				s22 = cur_state(i2, j2);
				// Arithmetic mean rule
				val = (  s00 +   s01 +   s02 + 
				         s10 +   s11 +   s12 + 
				         s20 +   s21 +   s22) / 9;
				if (val) 
					next_state(i1, j1) = (val + 1) % modulo;
			}
		}

		// produce image
		char fn[100];
		sprintf(fn, "rug%d-%dx%d-frame%05d.png", modulo, width, height, frame);
		if ((frame > 4050) || frame % 50 == 0)  // skip non-interesting frames
			save_matrix(next_state, fn, palette);

		// alter state
		cur_state_index ^= 1;
		next_state_index ^= 1;
	}
	return 0;
}