C++ Mandelbrot with gradiented colors Multithreaded generation

mandelbrot.cpp

// Mandelbrot.cpp : Defines the entry point for the console application.

//#include "stdafx.h" // if you are using visual studio
#include <fstream>
#include <iostream>
#include <vector>
#include <string> 
#include <sstream>
#include <thread>
using namespace std;

template <class T> struct RGB { T r, g, b; };

template <class T>
class Matrix {
public:
	Matrix(const size_t rows, const size_t cols) : _rows(rows), _cols(cols) {
		_matrix = new T*[rows];
		for (size_t i = 0; i < rows; ++i) {
			_matrix[i] = new T[cols];
		}
	}
	Matrix(const Matrix &m) : _rows(m._rows), _cols(m._cols) {
		_matrix = new T*[m._rows];
		for (size_t i = 0; i < m._rows; ++i) {
			_matrix[i] = new T[m._cols];
			for (size_t j = 0; j < m._cols; ++j) {
				_matrix[i][j] = m._matrix[i][j];
			}
		}
	}
	~Matrix() {
		for (size_t i = 0; i < _rows; ++i) {
			delete[] _matrix[i];
		}
		delete[] _matrix;
	}
	T *operator[] (const size_t nIndex)
	{
		return _matrix[nIndex];
	}
	size_t width() const { return _cols; }
	size_t height() const { return _rows; }
protected:
	size_t _rows, _cols;
	T **_matrix;
};

// Portable PixMap image
class PPMImage : public Matrix<RGB<unsigned char>>
{
public:
	PPMImage(const size_t height, const size_t width) : Matrix(height, width) { }
	void save(const std::string &filename)
	{
		std::ofstream out(filename, std::ios_base::binary);
		out << "P6" << std::endl << _cols << " " << _rows << std::endl << 255 << std::endl;
		for (size_t y = 0; y<_rows; y++)
			for (size_t x = 0; x<_cols; x++)
				out << _matrix[y][x].r << _matrix[y][x].g << _matrix[y][x].b;
	}
};

// transforms the given RGB collection by HSV (hue-saturation shifting)
RGB<unsigned char> transformHSV(
	const RGB<unsigned char> &in,  // color to transform
	float H,          // hue shift (in degrees)
	float S,          // saturation multiplier (scalar)
	float V           // value multiplier (scalar)
)
{
	float VSU = V * S*cos(H*acos(-1) / 180);
	float VSW = V * S*sin(H*acos(-1) / 180);

	RGB<unsigned char> ret;
	ret.r = (.299*V + .701*VSU + .168*VSW)*in.r
		+ (.587*V - .587*VSU + .330*VSW)*in.g
		+ (.114*V - .114*VSU - .497*VSW)*in.b;
	ret.g = (.299*V - .299*VSU - .328*VSW)*in.r
		+ (.587*V + .413*VSU + .035*VSW)*in.g
		+ (.114*V - .114*VSU + .292*VSW)*in.b;
	ret.b = (.299*V - .3*VSU + 1.25*VSW)*in.r
		+ (.587*V - .588*VSU - 1.05*VSW)*in.g
		+ (.114*V + .886*VSU - .203*VSW)*in.b;
	return ret;
}

// Calculate whether c belongs to the Mandelbrot set or not
int mandelbrot(double c_im, double c_re)
{
	// configurable
	const unsigned MaxIterations = 255;

	// Set Z = c
	double
		Z_re = c_re,
		Z_im = c_im;

	// iteration
	unsigned n;

	for (n = 0; n<MaxIterations; ++n)
	{
		double
			Z_re2 = pow(Z_re, 2),
			Z_im2 = pow(Z_im, 2);

		/* Z = Z2 + c */
		// addition of the two complex numbers
		Z_im = 2 * Z_re*Z_im + c_im;
		Z_re = Z_re2 - Z_im2 + c_re;

		// absolute value of Z is > 2
		// (are we within the imaginary magnitude boundary?
		// if (sqrt(Z_re*Z_re + Z_im * Z_im) > 2)
		// simplify to:
		if (Z_re2 + Z_im2 >= 4.0)
			break;
	}

	return n;
}

// interpolates two RGB collections by a value
RGB<unsigned char> interpolate(const RGB<unsigned char> &a, const RGB<unsigned char> &b, double value) {
	// (b - a) * value + a
	auto red = static_cast<unsigned char>((b.r - a.r) * value + a.r);
	auto green = static_cast<unsigned char>((b.g - a.g) * value + a.g);
	auto blue = static_cast<unsigned char>((b.b - a.b) * value + a.b);

	return { red, green, blue };
}

// generates a gradiented RGB by positions on interpolated gradiented colors
RGB<unsigned char> mandelbrot_color(double value) {

	// you could also try different positions here, but these work well
	const double positions[] = 
	{
		0.0,
		0.16,
		0.42,
		0.6425,
		0.8575,
	};

	// try different gradients here!
	const std::vector<RGB<unsigned char>> gradients =
	{
		//{ 0,   7,   100 },
		//{ 32,  107, 203 },
		//{ 237, 255, 255 },
		//{ 255, 170, 0 },
		//{ 0,   2,   0 },
		{3,0,30},
		{115,3,192},
		{236,56,188},
		{253,239,249},
		{255,255,255},
	};
	const int num_grads = sizeof(positions) / sizeof(double);

	double value_d = value * num_grads;
	auto value_i = static_cast<int>(value_d);

	if (value_i >= num_grads - 1) {
		return gradients[num_grads - 1];
	}
	else {
		return interpolate(gradients[value_i], gradients[value_i + 1], value_d - value_i);
	}
}

unsigned CUR_THREAD = 0;

// handles a mandelbrot point on the image and sets the pixel color appropiately on a certain 'region'
// region specifying a place x*y in size, with optional offsets
void mandelbrot_point(PPMImage *image, 
	unsigned int BLOCK_OFFSETY, unsigned int BLOCK_HEIGHT,
	unsigned int BLOCK_OFFSETX, unsigned int BLOCK_WIDTH,
	double MaxIm, double Im_factor, double MinRe, double Re_factor ) {
	
	unsigned MY_THREAD = ++CUR_THREAD;

	chrono::milliseconds start = chrono::duration_cast<chrono::milliseconds>(
		chrono::system_clock::now().time_since_epoch()
		);

	//loop y: height
	for (unsigned y = BLOCK_OFFSETY; y < BLOCK_OFFSETY + BLOCK_HEIGHT; ++y)
	{
		// get the imaginary c
		double c_im = MaxIm - y * Im_factor;

		// loop x
		for (unsigned x = BLOCK_OFFSETX; x < BLOCK_OFFSETX + BLOCK_WIDTH; ++x)
		{
			// get the real c
			double c_re = MinRe + x * Re_factor;

			// get the RGBA color value, then set the appropiate color
			int byteVal = mandelbrot(c_im, c_re);
			(*image)[y][x] = mandelbrot_color(1 - pow(0.99, byteVal));
		}
	}

	chrono::milliseconds end = chrono::duration_cast<chrono::milliseconds>(
		chrono::system_clock::now().time_since_epoch()
		);

	stringstream msg;
	msg << "Thread " << this_thread::get_id() << "(" << MY_THREAD << ") finished after " << (end-start).count() << " ms" << endl;
	cout << msg.str();
}

int main()
{
	const unsigned width = 1600;
	const unsigned height = 1600;

	PPMImage image(height, width);

	// Area definition:
	double MinRe = -2.0; // min real
	double MaxRe = 1.0; // max real
	double MinIm = -1.5; // min imaginary
	double MaxIm = MinIm + (MaxRe - MinRe)*height / width; // max imaginary
	// imaginary is calculated dynamically from other boundaries

	// Complex value for pixels inbetween:
	// double c_re = MinRe + x * (MaxRe - MinRe) / (width - 1);
	// double c_im = MaxIm - y * (MaxIm - MinIm) / (height - 1);

	// Define constants from that:
	double Re_factor = (MaxRe - MinRe) / (width - 1);
	double Im_factor = (MaxIm - MinIm) / (height - 1);
	// because these values, do not change: they are constant. so we turn them into factors
	// Now we can access the values easier:
	// double c_re = MinRe + x * Re_factor;
	// double c_im = MaxIm - y * Im_factor;

	// Here, you can set any number of threads.
	// Note: one thread is the same as no threading!
	// On my i5 4670k cpu, 3 threads perform just as well as 64
	// One thread (no multithreading) is more than twice as slow as multithreaded
	// To test out remaining pixels, you can try 8,24,etc. threads (square root has decimals!)
	const unsigned NUM_THREADS = 24;

	// the magical square root of our threads, used in calculations
	unsigned SQRT_THREADS = ceil(sqrt(NUM_THREADS));
	vector<thread> threads;

	// every thread handles a 'block' of the image, block meaning a specific section x*y section of the image
	unsigned BLOCK_WIDTH = (floor)(width / SQRT_THREADS);
	unsigned BLOCK_HEIGHT = (floor)(height / SQRT_THREADS);

	// sometimes, with a certain number of threads, the image can't be evenly divided into blocks 
	// without leaving a few pixels behind, so a there is can be a remainder of pixels
	unsigned REM_WIDTH = width - BLOCK_WIDTH * SQRT_THREADS;
	unsigned REM_HEIGHT = height - BLOCK_HEIGHT * SQRT_THREADS;

	// we divide de image into columns and rows
	// for example, if we consider the image to be 1600x1600 and use 64 threads
	// that means the image will be created by those 64 threads,
	// and every column and every row is worked on by the square root of 64 (threads)
	// so each column and row is worked on by 8 blocks (8*8 = 64)

	// loop columns
	for (unsigned col = 0; col < SQRT_THREADS; ++col)
	{
		//  the offset of each column is the width of a block * column count
		unsigned BLOCK_OFFSETX = col * BLOCK_WIDTH;

		// loop rows
		for (unsigned row = 0; row < SQRT_THREADS; ++row)
		{
			// the offset of each row is the height of a block * row count
			// modulated by the total height of the image
			// ensuring that the offset resets by every new column
			unsigned BLOCK_OFFSETY = row * BLOCK_HEIGHT % height;

			// push the thread that handles this block
			threads.push_back(
				thread(mandelbrot_point, &image,
					BLOCK_OFFSETY, BLOCK_HEIGHT,
					BLOCK_OFFSETX, BLOCK_WIDTH,
					MaxIm, Im_factor, MinRe, Re_factor));
		}
	}

	// we have remaining height
	if (REM_HEIGHT > 0)
	{
		// for every 'strand' of remaining height, process a thread to handle it
		for (unsigned h = 0; h < REM_HEIGHT; ++h)
		{
			// push the thread to handle this 'strand'
			threads.push_back(
				thread(mandelbrot_point, &image,
					height - REM_HEIGHT, 1,
					0, width,
					MaxIm, Im_factor, MinRe, Re_factor));
		}
	}
	// we have remaining width
	if (REM_WIDTH > 0)
	{
		// for every 'strand' of remaining width, process a thread to handle it
		for (unsigned w = 0; w < REM_WIDTH; ++w)
		{
			// push the thread to handle this 'strand'
			threads.push_back(
				thread(mandelbrot_point, &image,
					0, height,
					width - REM_WIDTH , 1,
					MaxIm, Im_factor, MinRe, Re_factor));
		}
	}

	// join all the threads with the main thread
	for (auto &thread : threads)
		thread.join();

	// save the image
	image.save("mandelbrot.ppm");

	// if you are on windows, or for some reason can not view thread times
	// uncomment this code to accept an input at the end to stop the process from exiting immediately
	//std::string str;
	//std::getline(std::cin, str);

	return 0;
}

Output: