/**
 * Mandelbulb.pbk
 * Last update: 30 November 2009
 * 
 * Changelog:
 *		1.0		- Initial release
 *
 * 
 * Copyright (c) 2009 Tom Beddard
 * http://www.subblue.com
 * Spherical polar reformulation: Chris King 2009
 * 
 * For more Flash and PixelBender based generative graphics experiments see:
 * http://www.subblue.com/blog
 * 
 * Licensed under the MIT License: 
 * http://www.opensource.org/licenses/mit-license.php
 *
 *
 * Credits and references
 * ======================
 * For the story behind the 3D Mandelbrot see the following page:
 * http://www.skytopia.com/project/fractal/mandelbulb.html
 * 
 * The original forum disussion with many implementation details can be found here:
 * http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/
 * 
 * This implementation references the 4D Quaternion GPU Raytracer by Keenan Crane
 * http://www.devmaster.net/forums/showthread.php?t=4448
 * and the NVIDIA CUDA/OptiX implementation by cbuchner1:
 * http://forums.nvidia.com/index.php?showtopic=150985
 * 
 */

#define BAILOUT 4.0
#define BOUNDING_RADIUS 1.4
#define PI 3.141592653

<languageVersion : 1.0;>

kernel Mandelbulb
<	namespace : "com.subblue.filters";
	vendor : "Tom Beddard";
	version : 1;
	description : "Mandelbulb Fractal Ray Tracer";
>


{
	parameter int antialiasing
	<
		minValue:1;
		maxValue:3;
		defaultValue:1;
		description:"Super sampling quality. Number of samples squared per pixel.";
	>;
	
	parameter bool phong
	<
		defaultValue:true;
		description: "Enable phong shading.";
	>;
	
	parameter bool julia
	<
		defaultValue:false;
		description: "Enable Julia set version.";
	>;
	

    parameter float shadows

	<

		minValue:0.0;

		maxValue:1.0;

		defaultValue:0.0;

	>;

    
	parameter float ambientOcclusion
	<
		minValue:0.0;
		maxValue:1.0;
		defaultValue:1.0;
		description: "Enable fake ambient occlusion factor based on the orbit trap.";
	>;
	
	parameter float ambientOcclusionEmphasis
	<
		minValue:0.0;
		maxValue:2.0;
		defaultValue:1.0;
		description: "Emphasise the structure edges based on the number of steps it takes to reach a point in the fractal.";
	>;
	
	parameter float epsilon
	<
		minValue:1e-8;
		maxValue:1e-5;
		defaultValue:1e-5;
		description: "The starting ray marching intersection threshold.";
	>;
	
	parameter float bounding
	<
		minValue:1.0;
		maxValue:4.0;
		defaultValue:1.3;
	>;
	
	parameter int power
	<
		minValue:-16;
		maxValue:16;
		defaultValue:8;
		description: "The power of the fractal.";
	>;
	
	
	parameter float3 julia_c
	<
		minValue:float3(-2, -2, -2);
		maxValue:float3(2, 2, 2);
		defaultValue:float3(0.56, 0.8, 0.32);
		description: "The c constant for Julia set fractals";
	>;
	
	parameter float3 cameraPosition
	<
		minValue:float3(-4, -4, -4);
		maxValue:float3(4, 4, 4);
		defaultValue:float3(0, -3, 0);
		description: "Camera position.";
	>;
	
	parameter float3 cameraPositionFine
	<
		minValue:float3(-0.1, -0.1, -0.1);
		maxValue:float3(0.1, 0.1, 0.1);
		defaultValue:float3(0, 0, 0);
		description: "Fine tune position.";
	>;
	
	parameter float3 cameraRotation
	<
		minValue:float3(-180, -180, -180);
		maxValue:float3(180, 180, 180);
		defaultValue:float3(0, 0, -90);
		description: "Pointing angle in each axis of the camera.";
	>;
	
	parameter float cameraZoom
	<
		minValue:0.0;
		maxValue:10.0;
		defaultValue:0.0;
		description: "Zoom the camera.";
	>;
	
	parameter float3 light
	<
		minValue:float3(-50, -50, -50);
		maxValue:float3(50, 50, 50);
		defaultValue:float3(33, -28, 40);
		description: "Position of point light.";
	>;
	
	parameter pixel4 colorBackground
	<
		minValue:pixel4(0, 0, 0, 0);
		maxValue:pixel4(1, 1, 1, 1);
		defaultValue:pixel4(0.13, 0.15, 0.17, 1);
		description: "Background colour.";
	>;
	
	parameter float3 colorDiffuse
	<
		minValue:float3(0, 0, 0);
		maxValue:float3(1, 1, 1);
		defaultValue:float3(0.5, 0.69, 1);
		description: "Diffuse/base colour.";
	>;
	
	parameter float3 colorAmbient
	<
		minValue:float3(0, 0, 0);
		maxValue:float3(1, 1, 1);
		defaultValue:float3(0.38, 0.22, 0.11);
		description: "Ambient light colour.";
	>;
	
	parameter float3 colorLight
	<
		minValue:float3(0, 0, 0);
		maxValue:float3(1, 1, 1);
		defaultValue:float3(0.85, 1, 1);
		description: "Light colour.";
	>;
	
	parameter float colorSpread
	<
		minValue:0.0;
		maxValue:1.0;
		defaultValue:0.0;
		description: "Vary the colour based on the normal direction.";
	>;
	
	parameter float specularity
	<
		minValue:0.0;
		maxValue:1.0;
		defaultValue:0.18;
		description: "Phone specularity";
	>;
	
	parameter float specularExponent
	<
		minValue:0.1;
		maxValue:50.0;
		defaultValue:5.0;
		description: "Phong shininess";
	>;
	

    parameter float3 rotation

	<

		minValue:float3(-180, -180, -180);

		maxValue:float3(180, 180, 180);

		defaultValue:float3(0, 10.8, 21.6);

		description: "Rotate the Mandelbulb in each axis.";

	>;

    
	parameter int maxIterations
	<
		minValue:1;
		maxValue:14;
		defaultValue:5;
		description: "More iterations reveal more detail in the fractal surface but takes longer to calculate.";
	>;
	
	parameter int stepLimit
	<
		minValue:10;
		maxValue:200;
		defaultValue:110;
		description: "The maximum number of steps a ray should take.";
	>;
	
	parameter float epsilonStep
	<
		minValue:0.1;
		maxValue:1.0;
		defaultValue:1.0;
		description: "Scale the epsilon step distance. Smaller values are slower but will generate smoother results.";
	>;
	
	parameter float epsilonScale
	<
		minValue:0.1;
		maxValue:1.0;
		defaultValue:0.775;
		description: "Scale the intersection threshold, espilon, with distance.";
	>;
	
	parameter int2 size
	<
		minValue:int2(100, 100);
		maxValue:int2(2048, 2048);
		defaultValue:int2(640, 480);
		description: "The output size in pixels.";
	>;
	
	//input	 image4 envmap;
	
	region generated()
	{
		return region(float4(0, 0, size.x, size.y));
	}
	
	dependent float bailout_sr, aspectRatio, sampleStep, sampleContribution, pixel_scale;
	dependent float3 eye, lightSource;
	dependent float3x3 viewRotation, objRotation;
	output pixel4 dst;
	
	
	// The fractal calculation
	//
	// Calculate the closest distance to the fractal boundary and use this
	// distance as the size of the step to take in the ray marching.
	// 
	// Fractal formula:
	//	  z' = z^p + c
	// 
	// For each iteration we also calculate the derivative so we can estimate 
	// the distance to the nearest point in the fractal set, which then sets the
	// maxiumum step we can move the ray forward before having to repeat the calculation.
	//
	//	 dz' = p * z^(p-1)

    // 
	float distanceEstimation(float3 z, inout float min_dist)
	{
		float3 c = julia ? julia_c : z; // Julia set has fixed c, Mandelbrot c changes with location
		float p = float(power);			// power
		float pd = p - 1.0;				// power for derivative

		// Convert z to polar coordinates
		float R	 = length(z);
		float th = atan(z.y, z.x);
		float ph = acos(z.z / R);
		
		// Record z orbit distance for ambient occulsion shading
		if (R < min_dist) min_dist = R;
		
		float3 dz;
		float ph_dz = 0.0;
		float th_dz = 0.0;
		float R_dz	= 1.0;
		
		// Iterate to compute the distance estimator.
		for (int i = 0; i < maxIterations; i++) {
			// Calculate derivative of 
			float powR = p * pow(R, pd);
			float powRsin = powR * R_dz * sin(ph_dz + pd*ph);
			dz.x = powRsin * cos(th_dz + pd*th) + 1.0;
			dz.y = powRsin * sin(th_dz + pd*th);
			dz.z = powR * R_dz * cos(ph_dz + pd*ph);
			
			// polar coordinates of derivative dz
			R_dz  = length(dz);
			th_dz = atan(dz.y, dz.x);
			ph_dz = acos(dz.z / R_dz);
			
			// z iteration
			powR = pow(R, p);
			powRsin = sin(p*ph);
			z.x = powR * powRsin * cos(p*th);
			z.y = powR * powRsin * sin(p*th);
			z.z = powR * cos(p*ph);
			z += c;
			
			R  = length(z);
			th = atan(z.y, z.x);
			ph = acos(z.z / R);
		    
			if (R < min_dist) min_dist = R;
			if (R > BAILOUT) break;
		}
		
		// Return the distance estimation value which determines the next raytracing
		// step size, or if whether we are within the threshold of the surface.
		return 0.5 * R * log(R)/R_dz;
	}
	
	
	// Intersect bounding sphere
	// 
	// If we intersect then set the tmin and tmax values to set the start and
	// end distances the ray should traverse.
	bool intersectBoundingSphere(float3 origin, 
								 float3 direction,
								 out float tmin,
								 out float tmax)
	{
		bool hit = false;
		
		float b = dot(origin, direction);
		float c = dot(origin, origin) - bounding;
		float disc = b*b - c;			// discriminant
		tmin = tmax = 0.0;
		
		if (disc > 0.0) {
			// Real root of disc, so intersection
			float sdisc = sqrt(disc);
			float t0 = -b - sdisc;			// closest intersection distance
			float t1 = -b + sdisc;			// furthest intersection distance
			
			if (t0 >= 0.0) {
				// Ray intersects front of sphere

                float min_dist;
				tmin = distanceEstimation(origin + t0 * direction, min_dist);
				tmax = t1;
			} else if (t0 < 0.0) {
				// Ray starts inside sphere
				float min_dist;
				tmin = distanceEstimation(origin, min_dist);
				tmax = t1;
			}
			hit = true;
		}
		
		return hit;
	}
	
	
	// Calculate the gradient in each dimension from the intersection point
	float3 estimate_normal(float3 z, float e)
	{
		float min_dst;	// Not actually used in this particular case
		float dx = distanceEstimation(z + float3(e, 0, 0), min_dst) - distanceEstimation(z - float3(e, 0, 0), min_dst);
		float dy = distanceEstimation(z + float3(0, e, 0), min_dst) - distanceEstimation(z - float3(0, e, 0), min_dst);
		float dz = distanceEstimation(z + float3(0, 0, e), min_dst) - distanceEstimation(z - float3(0, 0, e), min_dst);
		
		return normalize(float3(dx, dy, dz) / (2.0*e));
	}
	
	
	// Computes the direct illumination for point pt with normal N due to
	// a point light at light and a viewer at eye.
	float3 Phong(float3 pt, float3 N)
	{
		float3 diffuse	= float3(0);			// Diffuse contribution
		float3 specular = float3(0);			// Specular contribution
		float3 color	= float3(0);
		
		float3 L = normalize(light * objRotation - pt);	// find the vector to the light
		float  NdotL = dot(N, L);			// find the cosine of the angle between light and normal
		
		if (NdotL > 0.0) {
			// Diffuse shading
			diffuse = colorDiffuse * NdotL * colorLight;
            
			// Phong highlight
			float3 E = normalize(eye - pt);		// find the vector to the eye
			float3 R = L - 2.0 * NdotL * N;		// find the reflected vector
			float  RdE = dot(R,E);
			
			if (RdE <= 0.0) {
				specular = specularity * colorLight * pow(abs(RdE), specularExponent);
			}
		}
		
		color = (colorAmbient + diffuse + specular);

        // Add some of the normal to the color to make it more interesting
		color += abs(N) * colorSpread;
		
		return clamp(color, 0.0, 1.0);
	}
	
	
	// Define the ray direction from the pixel coordinates
	float3 rayDirection(float2 p)
	{
		float3 direction = float3( 2.0 * aspectRatio * p.x / float(size.x) - aspectRatio, 
								  -2.0 * p.y / float(size.y) + 1.0, 
								  -2.0 * exp(cameraZoom));
		return normalize(direction * viewRotation * objRotation);
	}
	
	
	// Calculate the output colour for each input pixel
	pixel4 renderPixel(float2 pixel)
	{
		float tmin, tmax;
		float3 ray_direction = rayDirection(pixel);
		pixel4 pixel_color = colorBackground;
		
		if (intersectBoundingSphere(eye, ray_direction, tmin, tmax)) {
			float3 ray = eye + tmin * ray_direction;
			
			float dist, ao;
			float min_dist = 4.0;
			float ray_length = tmin;
			float eps = epsilon;

			// number of raymarching steps scales inversely with factor
			int max_steps = int(float(stepLimit) / epsilonStep);
			int i;
			float f;
			
			for (i = 0; i < max_steps; ++i) {
				dist = distanceEstimation(ray, min_dist);

				// March ray forward
				f = epsilonStep * dist;
				ray += f * ray_direction;
				ray_length += f * dist;

				// Set the intersection threshold as a function of the ray length away from the camera
				eps = max(epsilon, pixel_scale * pow(ray_length, epsilonScale));
				
				// Are we within the intersection threshold or completely missed the fractal
				if (dist < eps || ray_length > tmax) break;
			}
			
			
			// Found intersection?
			if (dist < eps) {
				if (phong) {
					float3 normal = estimate_normal(ray, eps/2.0);
					pixel_color.rgb = Phong(ray, normal);

                    

                    if (shadows > 0.0) {

                        // The shadow ray will start at the intersection point and go

                        // towards the point light. We initially move the ray origin

                        // a little bit along this direction so that we don't mistakenly

                        // find an intersection with the same point again.

                        float3 light_direction = normalize((lightSource - ray) * objRotation);

                        ray += normal * eps * 2.0;

                        

                        float min_dist2;

                        dist = 4.0;

                        

                        for (int j = 0; j < max_steps; ++j) {

                            dist = distanceEstimation(ray, min_dist2);



                            // March ray forward

                            f = epsilonStep * dist;

                            ray += f * light_direction;

                            

                            // Are we within the intersection threshold or completely missed the fractal

                            if (dist < eps || dot(ray, ray) > bounding * bounding) break;

                        }

                        

                        // Again, if our estimate of the distance to the set is small, we say

                        // that there was a hit and so the source point must be in shadow.

                        if (dist < eps) pixel_color.rgb *= 1.0 - shadows;

                    }
				} else {

                    // Just use the base colour
					pixel_color.rgb = colorDiffuse;
				}
				
				ao	= 1.0 - clamp(1.0 - min_dist * min_dist, 0.0, 1.0) * ambientOcclusion;
				ao *= 1.0 - (float(i) / float(max_steps)) * ambientOcclusionEmphasis;
				pixel_color.rgb *= ao;
			}
		} else {
			//pixel_color = backgroundEnvironement(ray_direction);
		}
		
		return pixel_color;
	}
	
	
	// Common values used by all pixels
	void evaluateDependents()
	{
		aspectRatio = float(size.x) / float(size.y);
		
		// Camera orientation
		float c1 = cos(radians(-cameraRotation.x));
		float s1 = sin(radians(-cameraRotation.x));
		float3x3 viewRotationY = float3x3( c1, 0, s1, 
                                            0, 1, 0, 
                                          -s1, 0, c1);
		
		float c2 = cos(radians(-cameraRotation.y));
		float s2 = sin(radians(-cameraRotation.y));
		float3x3 viewRotationZ = float3x3( c2, -s2, 0,
                                           s2, c2, 0, 
                                            0, 0, 1);
		
		float c3 = cos(radians(-cameraRotation.z));
		float s3 = sin(radians(-cameraRotation.z));
		float3x3 viewRotationX = float3x3( 1, 0, 0,
                                           0, c3, -s3, 
                                           0, s3, c3);
		
		viewRotation = viewRotationX * viewRotationY * viewRotationZ;
		

        // Object rotation

        c1 = cos(radians(-rotation.x));

		s1 = sin(radians(-rotation.x));

		float3x3 objRotationY = float3x3( c1, 0, s1, 

                                            0, 1, 0, 

                                          -s1, 0, c1);

		

		c2 = cos(radians(-rotation.y));

		s2 = sin(radians(-rotation.y));

		float3x3 objRotationZ = float3x3( c2, -s2, 0,

                                           s2, c2, 0, 

                                            0, 0, 1);

		

		c3 = cos(radians(-rotation.z));

		s3 = sin(radians(-rotation.z));

		float3x3 objRotationX = float3x3( 1, 0, 0,

                                           0, c3, -s3, 

                                           0, s3, c3);

		

		objRotation = objRotationX * objRotationY * objRotationZ;

        

        
		//eye = float3(0, 0, camera.w) * viewRotation;
		//lightSource = light * viewRotation * 100.0;
		eye = (cameraPosition + cameraPositionFine) * objRotation;
		lightSource = light;
		
		// Super sampling
		sampleStep = 1.0 / float(antialiasing);
		sampleContribution = 1.0 / pow(float(antialiasing), 2.0);
		pixel_scale = 1.0 / max(float(size.x), float(size.y));
	}
	
	
	// The main loop
	void evaluatePixel()
	{
		pixel4 dst_color = pixel4(0, 0, 0, 0);
		
		if (antialiasing > 1) {
			// Average detailSuperSample^2 points per pixel
			for (float i = 0.0; i < 1.0; i += sampleStep)
				for (float j = 0.0; j < 1.0; j += sampleStep)
					dst_color += sampleContribution * renderPixel(float2(outCoord().x + i, outCoord().y + j));
		} else {
			dst_color = renderPixel(outCoord());
		}
		
		// Return the final color which is still the background color if we didn't hit anything.
		dst = dst_color;
	}
}