#include "noises.h"
#include "filterwidth.h"
#include "raysphere.h"

//! Macro
/**
* A macro to normalize a given value to a specified range.
* \param mi - the minimum value of the range
* \param mx - the maximum value of the range
* \param value - the value to be normalized
* \return the normalized value
*/
#define normalizeRange(mi,mx,value) (value-mi)/(mx-mi)

//! The Bark2 Shader
/**
* This file contains the Bark2 displacement shader. 
* @author Hannes Ricklefs - NCCA 2005
*/

//! A Function
/**
* Calculate the voronoi pattern. This function has been taken
* from the supporting material of the Advanced Renderman Book.
* It has been included to be more obvious to the author.
* \param P - position to test
* \param jitter - how much to jitter the cell center positions (1 is typical,
*             	  smaller values make a more regular pattern, larger values
*                 make a more jagged pattern; use jitter >1 at your risk!).
* \param f1, f2 - distance to the nearest feature (f1 is closest, f2 is the
*				  distance to the 2nd closest)
* \param pos1, pos2 - the position of the nearest feature
*/
void
voronoi(point P;
		 float jitter;
		 output float f1;  output point pos1;
		 output float f2;  output point pos2;
    )
{
	//cellcenter
    point thiscell = point (floor(xcomp(P))+0.5, floor(ycomp(P))+0.5,
			    floor(zcomp(P))+0.5);
	f1 = f2 = 1000;
    uniform float i, j, k;
    for (i = -1;  i <= 1;  i += 1) {
        for (j = -1;  j <= 1;  j += 1) {
            for (k = -1;  k <= 1;  k += 1) {
				
				point testcell = thiscell + vector(i,j,k);
			
				//featurepos
                point pos = testcell + 
		            jitter * (vector cellnoise (testcell) -0.5);
				
				vector offset = pos - P;
                
				float dist = offset . offset; /* actually dist^2 */
                
				if (dist < f1) {
                    f2 = f1;  pos2 = pos1;
                    f1 = dist;  pos1 = pos;
                } else if (dist < f2) {
                    f2 = dist;  pos2 = pos;
				}
    		}
		}
    }
    f1 = sqrt(f1);  f2 = sqrt(f2);
}

//! The displacement shader
/**
* This displacement shader provides the main functionality to create photo realitc looking
* tree bark. The main aspect is the use of a fractal voronoi pattern.
* The best approach in using this shader is to test different values for the crack_stretch vector
* as this vector provides the way of controlling the voronoi pattern. In order to stretch the pattern
* the user should apply a vector where the value of the direction to stretch in is less than the other 
* two directions. Depending on the given geometry these values could range between 0.1 - 1 or 1 - 20. 
* Once the desired displacement is achieved the user should use the provided normalrange.pl script to
* identify the range of displacement this pattern generate. The resulting min and max values should be 
* passed into lowNormaliseRange and topNormaliseRange respectively.
*
* Inputs: 
* \param Km 			- a scale value for the overall texture. Should be mostly in the range 0.05 - 0.2
* \param fBmoctaves 	- the octaves of the Fractional Brownian Motion in order to affect the cracks.
* \param fBmlacunarity 	- the lacunarity of the Fractional Brownian Motion in order to affect the cracks.
* \param fBmgain 		- the gain of the Fractional Brownian Motion in order to affect the cracks.
* \param turboctaves 	- the octaves of the Turbulence in order to add turbulence to the overall texture.
* \param turblcunarity 	- the lacunarity of the Turbulence in order to add turbulence to the overall texture.
* \param turbgain	 	- the gain of the Turbulence in order to add turbulence to the overall texture.
* \param barksmoothness - scales the Fractional Motion applied to the cracks
* \param squarness_of_bark_patches - will be passed to the jitter of the voronoi function jitter> 1 at your own risk
* \param crack_detail	- the amount of times the voronoi pattern should be applied, increase in high increases your render times 
* \param lowNormaliseRange - the min value calculated by the displacement inorder to normalise the output hump.
*							.The normlrange.pl script should be used to calculate the min value.
* \param topNormaliseRange - the max value calculated by the displacement inorder to normalise the output hump.
*							.The normlrange.pl script should be used to calculate the max value.
* \param space 			- the space in which the shader should be be applied
* \param crack_strech	- stretches the voronoi pattern by this vector
* 
* Outputs:
* \param hump			- the displacement normalise according to lowNormaliseRange and topNormaliseRange
* \param cellid			- the cellid of the voronoi pattern at this point
*/
displacement Bark2(		uniform float Km = 0.1; 
						float fBmoctaves=6,
						fBmlac=2,
						fBmgain=0.5,
						turboctaves=4,
						turblac=2,
						turbgain=0.5,
						barksmoothness=0.07,
						squareness_of_bark_patches=1,
						crack_detail=10,
						lowNormaliseRange=0.712768, 
						topNormaliseRange=2.883801;
						string space="shader";
						vector crack_stretch=(0.5,0.25,0.5);
						output varying float hump = 0;
						output varying point cellid = 0) {
						
	//Declare the positions for the nearest and second nearest feature.
	point pos1, pos2;
	
	//Declare the distance to the nearest and second nearest feature as 
	//well as the differance between the two.
	float f1, f2, diff;
	
	//Standard shader initialisation calls
	normal n = normalize(N); 
	normal nf = faceforward(n,I);
	point PP = transform(space,P);
	
	//Transform the vector into the relevant space inorder to stretch the voronoi pattern
	vector cs = vtransform(space,crack_stretch);	
	
	// Noise for the cracks 
	PP = PP + (barksmoothness*vfBm(2*PP,filterwidthp(PP),fBmoctaves,fBmlac,fBmgain));
	
	//stretch the voronoi pattern
	PP = point (comp(cs,0) * comp(PP,0), comp(cs,1) * comp(PP,1), comp(cs,2) * comp(PP,2));
	
	//apply the voronoi function
	voronoi(PP,squareness_of_bark_patches,f1,pos1,f2,pos2);
	
	
	/* The voronoi cell */					
	cellid = pos1 * 10;

	/* The difference */
	diff = 1-(f2-f1);
	
	//Values for the fractal voronoi pattern
	float i = 0;
	//The amplitude
	float amp = 1.0;	
	
	/* fractal higher frequnecy and lower amplitude */
	for(i = 0; i < crack_detail; i+=1) {
		
		PP = transform(space,P);
		/* Noise for the cracks */
		PP = PP + (barksmoothness*vfBm(2*PP,filterwidthp(PP),fBmoctaves,fBmlac,fBmgain));
		PP = point (comp(cs,0) * comp(PP,0), comp(cs,1) * comp(PP,1), comp(cs,2) * comp(PP,2));
		voronoi(PP,squareness_of_bark_patches,f1,pos1,f2,pos2);

		/* The difference */
		diff = 1-(f2-f1);
	
		/* turbulence for the humps */
		hump += amp*  (1-(smoothstep(0.77,0.95,diff)) + turbulence(PP,filterwidthp(PP),turboctaves,turblac,turbgain));
		
		cs = cs * 1.56;
		amp *= 0.6;
	}
	
	/* Recalculte the surface */
	P = P + n * hump * Km ;
	
	/* Normalize the output of hump for the color */
	//printf ("%f\n",hump);
	hump = normalizeRange(lowNormaliseRange, 
							topNormaliseRange,
							hump);
	//Recalculate the values
	N = calculatenormal(P);
	n = normalize(N);
	nf = faceforward(n,I);
}