/*
   Copyright (c) 2015 Techmeology

   Permission is hereby granted, free of charge, to any person obtaining a copy
   of this software and associated documentation files (the "Software"), to deal
   in the Software without restriction, including without limitation the rights
   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   copies of the Software, and to permit persons to whom the Software is
   furnished to do so, subject to the following conditions:

   The above copyright notice and this permission notice shall be included in
   all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   SOFTWARE.
*/

/*
 * Calculate something funny based on the Mandelbrot set
 * g++ --std=c++11 -fopenmp -O3 -s -o mandelbrot mandelbrot.cpp
 * ./mandelbrot <width> <height> [oversample] [maximum_iteration] [x0 x1 y0 y1]
 */

#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>

namespace
{
	const double PI = 3.141592654;

	double ColourFunctionInterpolate(double v, double min, double max)
	{
		return (v - min)/(max - min);
	}

	void ColourFunction(double iteration, double argument, double &red, double &green, double &blue)
	{
		red =
			(argument < 1.0/6.0) ? 1.0 :
			(argument < 2.0/6.0) ? ColourFunctionInterpolate(argument, 2.0/6.0, 1.0/6.0) :
			(argument < 4.0/6.0) ? 0.0 :
			(argument < 5.0/6.0) ? ColourFunctionInterpolate(argument, 4.0/6.0, 5.0/6.0) :
			1.0;
		green =
			(argument < 1.0/6.0) ? ColourFunctionInterpolate(argument, 0.0/6.0, 1.0/6.0) :
			(argument < 3.0/6.0) ? 1.0 :
			(argument < 4.0/6.0) ? ColourFunctionInterpolate(argument, 4.0/6.0, 3.0/6.0) :
			0.0;
		blue =
			(argument < 2.0/6.0) ? 0.0 :
			(argument < 3.0/6.0) ? ColourFunctionInterpolate(argument, 2.0/6.0, 3.0/6.0) :
			(argument < 5.0/6.0) ? 1.0 :
			ColourFunctionInterpolate(argument, 6.0/6.0, 5.0/6.0);
	
		red *= 1.0 - iteration;
		green *= 1.0 - iteration;
		blue *= 1.0 - iteration;
	}

	void CalculatePoint(double x0, double y0, unsigned int max_iteration, double &iteration, double &argument) //returns iteration,argument
	{
		//Calculate x, y, and i in the standard way (see http://en.wikipedia.org/wiki/Mandelbrot_set#Escape_time_algorithm)
		double x = 0.0;
		double y = 0.0;
		unsigned int i = 0;
		while ((x*x + y*y < 4.0) && (i < max_iteration)){
			double t = x*x - y*y + x0;
			y = 2.0*x*y + y0;
			x = t;
			i++;
		}
	
		//Calculate iteration in the range 0 to 1
		iteration = (double)i/(double)max_iteration;
	
		//Calculate the argument
		argument = (atan2(y, x) + PI)/(2.0 * PI);
	}

	void CalculateImage(double *image, unsigned int width, unsigned int height, unsigned int max_iteration, double min_x, double max_x, double min_y, double max_y)
	{
		#pragma omp parallel for
		for (unsigned int yi = 0; yi < height; yi++){
			for (unsigned int xi = 0; xi < width; xi++){
				//Calculate x and y in mathematical space rather than image space
				double x = min_x + (max_x - min_x) * (double)xi/(double)width;
				double y = min_y + (max_y - min_y) * (double)yi/(double)height;
			
				//Calculate the pixel
				double iteration, argument, red, green, blue;
				CalculatePoint(x, y, max_iteration, iteration, argument); //Calculate the information from which the pixel will be coloured
				ColourFunction(iteration, argument, red, green, blue); //Calculate the actual colour
			
				//Write the pixel
				image[(yi * width + xi) * 3] = red;
				image[(yi * width + xi) * 3 + 1] = green;
				image[(yi * width + xi) * 3 + 2] = blue;
			}
		}
	}

	void WriteImage(double *image, unsigned int width, unsigned int height, unsigned int oversample)
	{
		//Write the header
		printf("P6\n#A thing based on the mandelbrot set\n%u %u\n255\n", width, height);
	
		//Iterate over each pixel
		for (unsigned int y0 = 0; y0 < height; y0++){
			for (unsigned int x0 = 0; x0 < width; x0++){
				double red = 0.0;
				double green = 0.0;
				double blue = 0.0;
			
				//Iterate over each sub pixel
				for (unsigned int y1 = 0; y1 < oversample; y1++){
					for (unsigned int x1 = 0; x1 < oversample; x1++){
						red += image[((y0 * oversample + y1) * width * oversample + x0 * oversample + x1) * 3];
						green += image[((y0 * oversample + y1) * width * oversample + x0 * oversample + x1) * 3 + 1];
						blue += image[((y0 * oversample + y1) * width * oversample + x0 * oversample + x1) * 3 + 2];
					}
				}
			
				//Scale between 0 and 1
				red /= oversample * oversample;
				green /= oversample * oversample;
				blue /= oversample * oversample;
			
				//Write the pixels
				uint8_t red_i = red * 255.0;
				uint8_t green_i = green * 255.0;
				uint8_t blue_i = blue * 255.0;
				fwrite(&red_i, 1, 1, stdout);
				fwrite(&green_i, 1, 1, stdout);
				fwrite(&blue_i, 1, 1, stdout);
			}
		}
	}

	void PrintUsage(const char *name)
	{
		fprintf(stderr, "Usage is: %s <width> <height> [oversample] [maximum_iteration] [x0 x1 y0 y1]\n", name);
	}
}

int main(int argc, char **argv)
{
	//Check the basic argument requirements
	if (argc < 3){
		PrintUsage(argv[0]);
		return 1;
	}
	
	//Options
	unsigned int width = atoi(argv[1]);
	unsigned int height = atoi(argv[2]);
	unsigned int oversample = 1; //Default oversample is 1
	unsigned int max_iteration = 255; //Default maximum iteration is 255
	
	//Image window size
	double min_x = -2.5;
	double max_x = 1.0;
	double min_y = -1.0;
	double max_y = 1.0;
	
	//User specified an oversample option
	if (argc >= 4){
		oversample = atoi(argv[3]);
	}
	
	//User specified maximum iteration option
	if (argc >= 5){
		max_iteration = atoi(argv[4]);
	}
	
	if (argc == 9){
		min_x = atof(argv[5]);
		max_x = atof(argv[6]);
		min_y = atof(argv[7]);
		max_y = atof(argv[8]);
	}
	else if (argc > 5){
		PrintUsage(argv[0]);
		return 1;
	}
	
	//Check the validity of the arguments
	if (!width || !height || !oversample || !max_iteration){
		PrintUsage(argv[0]);
		return 1;
	}
	
	//Print what we're using
	fprintf(stderr,
		"[*] Wdith = %u\n"
		"[*] Height = %u\n"
		"[*] Oversample = %u\n"
		"[*] Maximum iteration = %u\n",
		width, height, oversample, max_iteration
	);
	
	//Actually do stuff
	double *image = new double[width * height * oversample * oversample * 3]; //Allocate memory
	CalculateImage(image, width * oversample, height * oversample, max_iteration, min_x, max_x, min_y, max_y); //Calculate the mandelbrot set image
	WriteImage(image, width, height, oversample); //Write the image
	delete []image; //Deallocate the memory
	
	return 0; //End of program
}