/* Copyright (c) 1999, TNCCI Inc.
 * 
 * musgrave.c
 *
 * a procedural texture plugin for blender that generates (and bumps !)
 * Kent Musgrave's fractal noise patterns.
 *
 * (source written/compiled by rob haarsma, between dec 95 and april 1999)
 *
 * Contact:      rob@captainvideo.nl   
 * Information:  http://www.captainvideo.nl/blender
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 * 
 */

#include "math.h"
#include "plugin.h"

#define NR_TYPES	1

#define N 0x1000
#define B 0x100
#define BM 0xff
#define NP 12   /* 2^N */
#define NM 0xfff
#define lerp(t, a, b) ( a + t * (b - a) )
#define s_curve(t) ( t * t * (3. - 2. * t) )
#define setup(i,b0,b1,r0,r1)\
	t = vec[i] + N;\
	b0 = ((int)t) & BM;\
	b1 = (b0+1) & BM;\
	r0 = t - (int)t;\
	r1 = r0 - 1.;

static float p[B + B + 2];
static float g3[B + B + 2][3];
static float g2[B + B + 2][2];
static float g1[B + B + 2];

float p_jit[8][2] = {	0.968813, 0.018570,
			0.344245, 0.147180,
			0.719306, 0.261499, 
			0.106714, 0.389051,
			0.475301, 0.513908,
			0.843805, 0.647431,
			0.227145, 0.769918, 
			0.595909, 0.888710 };



float result[8];
float cfra;

static	int rseed = 1;
static	int start = 1;
float	x, y, z, x1;

 /* timeconsuming oversampled rendering is disabled,
   (i didn't notice much difference) */
int	do_osa = FALSE;


extern float hnoise(float noisesize, float x, float y, float z);

float Noise3(float vec[]);


/* plugin menu and variables******************************************* */

char name[]= "Fractal";
char stnames[NR_TYPES][16]= {"Musgrave" };

VarStruct varstr[]= {
	/* butcode,	name,		default,min,	max	tooltip comment*/
	{ NUM|FLO,	"H",	    	0.0,	-1.0,	1.0,	
	"Musgrave plugin: Play with it..."},
	{ NUM|FLO,	"lacun",	1.4,    0.0,	10.0,	
	"Musgrave plugin: Play with it..."},  
	{ NUM|FLO,	"octav",	1.05,	0.0,	5.0,	
	"Musgrave plugin: Play with it..."}, 
	{ NUM|FLO,	"offset",	0.6,	-1.0,	2.0,	
	"Musgrave plugin: Play with it..."}, 
	{ NUM|FLO,	"gain",	    	5.0,	0.0,	10.0,	
	"Musgrave plugin: Play with it..."}, 
	{ LABEL,		"",		0,	 0,	0,	""},
	{ NUM|INT,	"seed",	    	1.0,	1.0,	1000.0,	
	"Musgrave plugin: Random seed number for the pattern"},
	{ NUM|FLO,	"scale",	1.3,	0.0,	10.0,	
	"Musgrave plugin: Scale pattern"}, 
	{ TOG|INT,	"bump it!",	-0.5,	-1.0,	1.0,	
	"Musgrave plugin: Activate bump mapping"}, 
	{ NUM|FLO,	"bump size",	0.5,	-5.0,	5.0,	
	"Musgrave plugin: Bump size"}, 
	{ NUM|FLO,	"bump factor",	0.009,	0.0,	0.05,	
		"Musgrave plugin: Bump sharpness"}, 
	{ LABEL,		"phaseIV  '99",	0,	 0,	0,	""},
};

typedef struct Cast {
	float H, lacun, octav, offst, gain, dum1;
	int seed;
	float scale;
	int bump_it;
	float  bsize, bfac, dum2;
} Cast;

/* ******************************************************************* */

int plugin_tex_doit(int, Cast*, float*, float*, float*);


int plugin_tex_getversion(void) 
{	
	return B_PLUGIN_VERSION;
}


void plugin_but_changed(int but) 
{
}


void plugin_init(void)
{
}


void plugin_getinfo(PluginInfo *info)
{
	info->name= name;
	info->stypes= NR_TYPES;
	info->nvars= sizeof(varstr)/sizeof(VarStruct);
	
	info->snames= stnames[0];
	info->result= result;
	info->cfra= &cfra;
	info->varstr= varstr;

	info->init= plugin_init;
	info->tex_doit=  (TexDoit) plugin_tex_doit;
	info->callback= plugin_but_changed;
}


float ridged(Cast *cast, float *texvec, float xafw, float yafw, float zafw)
{
      float            frequency, Noise3(), point[3], octaves, Hh;
      float            result, signal, weight;
      int               i;
      static float     Hfirst = -3.0;
      static float      exponent_array[10];

	point[0] = (texvec[0]+xafw) * cast->scale;
	point[1] = (texvec[1]+yafw) * cast->scale;
	point[2] = (texvec[2]+zafw) * cast->scale;

	octaves = cast->octav;
	Hh = cast->H;
		frequency = 1.0;
	if(Hfirst != Hh){
		for (i=0; i<=octaves+1; i++) {
			  /* compute weight for each frequency */
			  exponent_array[i] = pow( frequency, -Hh );
			  frequency *= cast->lacun;
		}
		Hfirst = Hh;
	}

	/* get first octave */
	signal = Noise3( point );
	/* get absolute value of signal (this creates the ridges) */
	if ( signal < 0.0 ) signal = -signal;
	/* invert and translate (note that "cast->offst" should be ~= 1.0) */
	signal = cast->offst - signal;
	/* square the signal, to increase "sharpness" of ridges */
	signal *= signal;
	/* assign initial values */
	result = signal;
	weight = 1.0;
	
	for( i=1; i<octaves + 1; i++ ) {
		/* increase the frequency */
				point[0] *= cast->lacun;
				point[1] *= cast->lacun;
				point[2] *= cast->lacun;
		
		/* weight successive contributions by previous signal */
		weight = signal * cast->gain;
		if ( weight > 1.0 ) weight = 1.0;
		if ( weight < 0.0 ) weight = 0.0;
		signal = Noise3( point );
		if ( signal < 0.0 ) signal = -signal;
		signal = cast->offst - signal;
		signal *= signal;
		/* weight the contribution */
		signal *= weight;
		result += signal * exponent_array[i];
	}
	return result;
}


static void normalize2(float v[2])
{
	float s;

	s = sqrt(v[0] * v[0] + v[1] * v[1]);
	v[0] = v[0] / s;
	v[1] = v[1] / s;
}


static void normalize3(float v[3])
{
	float s;

	s = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
	v[0] = v[0] / s;
	v[1] = v[1] / s;
	v[2] = v[2] / s;
}


static void init(void)
{
	int i, j, k;

	srand((unsigned int)rseed);
	for (i = 0 ; i < B ; i++) {
		p[i] = i;

		g1[i] = (float)((rand() % (B + B)) - B) / B;

		for (j = 0 ; j < 2 ; j++){
			g2[i][j] = (float)((rand() % (B + B)) - B) / B;
		}
		normalize2(g2[i]);

		for (j = 0 ; j < 3 ; j++){
			g3[i][j] = (float)((rand() % (B + B)) - B) / B;
		}
		normalize3(g3[i]);
	}

	while (--i) {
		k = p[i];
		p[i] = p[j = rand() % B];
		p[j] = k;
	}

	for (i = 0 ; i < B + 2 ; i++) {
		p[B + i] = p[i];
		g1[B + i] = g1[i];
		for (j = 0 ; j < 2 ; j++)
			g2[B + i][j] = g2[i][j];
		for (j = 0 ; j < 3 ; j++)
			g3[B + i][j] = g3[i][j];
	}
}


float Noise3(float vec[3])
{
	int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
	float rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
	register i, j;

	if (start) {
		start = 0;
		init();
	}

	setup(0, bx0,bx1, rx0,rx1);
	setup(1, by0,by1, ry0,ry1);
	setup(2, bz0,bz1, rz0,rz1);

	i = p[ bx0 ];
	j = p[ bx1 ];

	b00 = p[ i + by0 ];
	b10 = p[ j + by0 ];
	b01 = p[ i + by1 ];
	b11 = p[ j + by1 ];

	t  = s_curve(rx0);
	sy = s_curve(ry0);
	sz = s_curve(rz0);

#define at3(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] )

	q = g3[ b00 + bz0 ] ; u = at3(rx0,ry0,rz0);
	q = g3[ b10 + bz0 ] ; v = at3(rx1,ry0,rz0);
	a = lerp(t, u, v);

	q = g3[ b01 + bz0 ] ; u = at3(rx0,ry1,rz0);
	q = g3[ b11 + bz0 ] ; v = at3(rx1,ry1,rz0);
	b = lerp(t, u, v);

	c = lerp(sy, a, b);

	q = g3[ b00 + bz1 ] ; u = at3(rx0,ry0,rz1);
	q = g3[ b10 + bz1 ] ; v = at3(rx1,ry0,rz1);
	a = lerp(t, u, v);

	q = g3[ b01 + bz1 ] ; u = at3(rx0,ry1,rz1);
	q = g3[ b11 + bz1 ] ; v = at3(rx1,ry1,rz1);
	b = lerp(t, u, v);

	d = lerp(sy, a, b);

	return lerp(sz, c, d);
}


int plugin_tex_doit(int stype, Cast *cast, float *texvec, float *dxt, float *dyt)
{
	int  i;
	float base, osabase = 0.0, bump[3], osavec[3];
	
	if(rseed != cast->seed){
		rseed = cast->seed;
		start = 1;
	}

	bump[0] = 0.0;
	bump[1] = 0.0;
	bump[2] = 0.0;
	if((dxt || dyt) && do_osa){
		for(i = 0; i < 8; i++){
			osavec[0] = texvec[0] + (p_jit[i][0] * dxt[0] + p_jit[i][1] * dyt[0]);
			osavec[1] = texvec[1] + (p_jit[i][0] * dxt[1] + p_jit[i][1] * dyt[1]);
			osavec[2] = texvec[2] + (p_jit[i][0] * dxt[2] + p_jit[i][1] * dyt[2]);
			base= ridged(cast, osavec, 0.0, 0.0, 0.0);
			osabase += base;
			if(cast->bump_it){
				bump[0] += (base - ridged(cast, osavec, cast->bfac, 0.0, 0.0) )* -cast->bsize;
				bump[2] += (- base  + ridged(cast, osavec, 0.0, cast->bfac, 0.0) )* -cast->bsize;
				bump[1] += (base - ridged(cast, osavec, 0.0, 0.0, cast->bfac) )* -cast->bsize;
			}
		}
		base = osabase / 8.0;
		result[5] = bump[0] / 8.0;
		result[6] = bump[1] / 8.0;
		result[7] = bump[2] / 8.0;
	} else {
		base= ridged(cast, texvec, 0.0, 0.0, 0.0);
		if(cast->bump_it){
			result[5] = bump[0]= (base - ridged(cast, texvec, cast->bfac, 0.0, 0.0) )* -cast->bsize;
			result[6] = bump[2]= (- base  + ridged(cast, texvec, 0.0, cast->bfac, 0.0) )* -cast->bsize;
			result[7] = bump[1]= (base - ridged(cast, texvec, 0.0, 0.0, cast->bfac) )* -cast->bsize;
			/* VecAddf(result+5, result+5, bump); */
			/* Normalise(result+5); */
		}
	}
	result[0]= base;
	result[1]= result[0];
	result[2]= result[0];
	result[3]= result[0];
	result[4]= 1.0;

	if(cast->bump_it) {
		return 2;
	} else
		return 1;
}