// Noise class for generating pseudorandom noise fields
//based on noise lattice a la Peachey/Perlin
#include "points.h"
#include "aux.h"
#include <assert.h>
#include <iostream>

Noise::Noise()//construct a noise object
{		
	int i;
	index = new unsigned char[256]; assert(index);
	for(i = 0; i < 256; i++) index[i] = i;//fill array with indices
	for(i = 0; i < 256; i++) // shuffle it
	{
		int which = rand() % 256; // choose random place in array
		unsigned char tmp = index[which]; // swap them
		index[which] = index[i];
		index[i] = tmp;
	}
	noiseTable = new float[256]; assert(noiseTable);
	for(i = 0; i < 256; i++) noiseTable[i] = rand()/32767.99;
} // end of constructor

float Noise::noise(float scale, Point3& p)
{ // linearly interpolated lattice noise
	#define lerp(f, A, B) A + f * (B - A)
	float d[2][2][2];
	Point3 pp;
	pp.x = p.x * scale + 10000; //offset avoids negative values
	pp.y = p.y * scale + 10000; 
	pp.z = p.z * scale + 10000;
	long ix = (long)pp.x; long iy = (long)pp.y; long iz = (long)pp.z;
	float tx,ty,tz, x0,x1,x2,x3, y0,y1;
	tx = pp.x - ix; ty = pp.y - iy; tz = pp.z - iz; // fractional parts
	float mtx = 1.0 - tx, mty = 1.0 - ty, mtz = 1.0 - tz;

	for(int k = 0; k <= 1; k++) // get noise at 8 lattice points
	for(int j = 0; j <= 1; j++)
	for(int i = 0; i <= 1; i++)
	d[k][j][i] = latticeNoise(ix + i, iy + j,iz + k);

	x0 = lerp(tx, d[0][0][0],d[0][0][1]);
	x1 = lerp(tx, d[0][1][0],d[0][1][1]);
	x2 = lerp(tx, d[1][0][0],d[1][0][1]);
	x3 = lerp(tx, d[1][1][0],d[1][1][1]);
	y0 = lerp(ty, x0, x1);
	y1 = lerp(ty, x2, x3);
	return lerp(tz, y0, y1);
}

float Noise::turbulence(float s, Point3& p)
{
	float val = noise(s    , p) / 2 +
		   noise(s * 2, p) / 4 +
		   noise(s * 4, p) / 8 +
		   noise(s * 8, p) / 16;
	return val;
}

float Noise::marble(float strength,Point3& p)
{
	float turbul = turbulence(10, p);
	float val = sin(6 * p.z + strength * turbul);
	return mySpline(val);
}

float Noise::mySpline(float x) // used for marble
{
	//yap: replace SQR by sqrt:
	if(x <-0.4) return 0.15 +  2.857 * sqrt(x + 0.75);
	else if(x < 0.4) return 0.95 - 2.8125 * sqrt(x);
	else return 0.26 + 2.666 * sqrt(x - 0.7);
}

float Noise::latticeNoise(int i, int j, int k)
{ // return PR noise value on an integer lattice
	#define PERM(x) index[(x) & 255]
	#define INDEX(ix, iy, iz) PERM( (ix) + PERM((iy) + PERM(iz)) )
	return noiseTable[INDEX(i,j,k)];
}

// uses several of the basic classes developed above
// uses classes: DefUnit, DefUnitStack 
// uses Affine4 and associated classes
// uses Shapes class as well

//@@@@@@@@@@@@@  Scene class @@@@@@@@@@@@@@@@@@@@

// @@@@@@@@@@@@@@@@@@@@@ Affine4 class @@@@@@@@@@@@@@@@
Affine4::Affine4(){ // make identity transform
	m[0] = m[5]  = m[10] = m[15] = 1.0;
	m[1] = m[2]  = m[3]  = m[4]  = 0.0;
	m[6] = m[7]  = m[8]  = m[9]  = 0.0;
	m[11]= m[12] = m[13] = m[14] = 0.0;
}

void Affine4::setIdentityMatrix(){ // make identity transform
	m[0] = m[5]  = m[10] = m[15] = 1.0;
	m[1] = m[2]  = m[3]  = m[4]  = 0.0;
	m[6] = m[7]  = m[8]  = m[9]  = 0.0;
	m[11]= m[12] = m[13] = m[14] = 0.0;
}

void Affine4::set(Affine4 a)// set this matrix to a
{
	for(int i = 0; i < 16; i++)
		m[i]=a.m[i];
}

	//<<<<<<<<<<<< postMult >>>>>>>>>>>
void Affine4::postMult(Affine4 n){// postmultiplies this with n
	float sum;
	Affine4 tmp; 
	tmp.set(*this); // tmp copy
	for(int c = 0; c < 4; c++)// form this = tmp * n
		for(int r = 0; r <4 ; r++)
		{
			sum = 0;
			for(int k = 0; k < 4; k++)
				sum += tmp.m[4 * k + r]* n.m[4 * c + k];
			m[4 * c + r] = sum;
		}// end of for loops
}

//@@@@@@@@@@ AffineNode class @@@@@@@@@@@
AffineNode::AffineNode()
{
	next = NULL;
	affn = new Affine4; // new affine with identity in it
	invAffn = new Affine4; // and for the inverse
}

AffineNode::~AffineNode() //destructor
{
	delete affn;
	delete invAffn;
}

//@@@@@@@@@@@@@@@@ AffineStack class @@@@@@@@@@@@
AffineStack::AffineStack()//default constructor;puts identity on top
{
	tos = new AffineNode; // node with identity in it
	tos->next = NULL;
}

void AffineStack::dup()
{
	AffineNode* tmp = new AffineNode;			
	tmp->affn = new Affine4(*(tos->affn));
	tmp->invAffn = new Affine4(*(tos->invAffn));
	tmp->next = tos;
	tos = tmp;
}

void AffineStack::setIdentity() // make top item the identity matrix
{ 
	assert(tos != NULL);
	tos->affn->setIdentityMatrix();
	tos->invAffn->setIdentityMatrix();
}

void AffineStack::popAndDrop()
{
	if(tos == NULL) return; // do nothing
	AffineNode *tmp = tos;
	tos = tos->next;
	delete tmp; // should call destructor, which deletes trices
}	

void AffineStack::releaseAffines()
{ // pop and drop all remaining items 
	while(tos) popAndDrop();
}

void AffineStack::rotate(float angle, Vector3 u)
{
	Affine4 rm; // make identity matrix
	Affine4 invRm;
	u.normalize(); // make the rotation axis unit length
	float ang = angle * 3.14159265/ 180; // deg to  
	float c = cos(ang), s = sin(ang);
	float mc = 1.0 - c;
 	//fill the 3x3 upper left matrix - Chap.5 p. 29
	rm.m[0] = c + mc * u.x * u.x;
	rm.m[1] = mc * u.x * u.y + s * u.z;
	rm.m[2] = mc * u.x * u.z - s * u.y;
	rm.m[4] = mc * u.y * u.x - s * u.z;
	rm.m[5] = c + mc * u.y * u.y;
	rm.m[6] = mc * u.y * u.z + s * u.x;
	rm.m[8] = mc * u.z * u.x + s * u.y;
	rm.m[9] = mc * u.z * u.y - s * u.x;
	rm.m[10] = c + mc * u.z * u.z;
	// same for inverse : just sign of s is changed
	invRm.m[0] = c + mc * u.x * u.x;
	invRm.m[1] = mc * u.x * u.y - s * u.z;
	invRm.m[2] = mc * u.x * u.z + s * u.y;
	invRm.m[4] = mc * u.y * u.x + s * u.z;
	invRm.m[5] = c + mc * u.y * u.y;
	invRm.m[6] = mc * u.y * u.z - s * u.x;
	invRm.m[8] = mc * u.z * u.x - s * u.y;
	invRm.m[9] = mc * u.z * u.y + s * u.x;
	invRm.m[10] = c + mc * u.z * u.z;
	tos->affn->postMult(rm);
	// yap: just comment out to avoid compiler error
	// tos->invAffn->preMult(invRm);
}

void AffineStack::scale(float sx, float sy, float sz)
{ // post-multiply top item by scaling
	#define sEps 0.00001
	Affine4 scl;// make an identity
	Affine4 invScl;
	scl.m[0]  = sx; 
	scl.m[5]  = sy; 
	scl.m[10] = sz;// adjust it to a scaling matrix
// yap: syntax error here -- so I just fixed it w/o regard to logic
	if (fabs(sx) < sEps || fabs(sy) < sEps || fabs(sz) < sEps) 
	{
		cerr << "degenerate scaling transformation!\n"; 
	};
	invScl.m[0]  = 1/sx; invScl.m[5]  = 1/sy; invScl.m[10] = 1/sz;
	tos->affn->postMult(scl); //
	// tos->invAffn->preMult(invScl);
}		

void AffineStack::translate(Vector3 d)
{
	Affine4 tr; // make identity matrix
	Affine4 invTr;
	tr.m[12] = d.x; tr.m[13] = d.y;	tr.m[14] = d.z;
	invTr.m[12] = -d.x;	invTr.m[13] = -d.y; invTr.m[14] = -d.z;
	tos->affn->postMult(tr);
	// tos->invAffn->preMult(invTr);
}



//@@@@@@@@@@@@@@@@@@@ PointCluster class @@@@@@@@@@@@@@@@@@
PointCluster::PointCluster()
{num = 0; pt = NULL;}

PointCluster::PointCluster(int n)// make a cluster of n points
{
	pt = new Point3[n]; assert(pt);
	num = n;
}

//@@@@@@@@@@@@@@@@@@@@@@@@@@ SphereInfo @@@@@@@@@@@@@@@@@@@@@
void SphereInfo::set(float x, float y, float z, float rsq)
{
	center.set(x,y,z);
	radSq = rsq;
}

//@@@@@@@@@@@@@@@@@@@@@@@@ Cuboid @@@@@@@@@@@@@@@@@@@@@@@@
void Cuboid::set(float l, float t, float r, float b, float f, float bk)
{
	left = l; top = t; right = r; bott = b; front = f; back = bk;
}

void Cuboid::set(Cuboid& c)
{
	left = c.left;top = c.top; right = c.right; bott = c.bott;
	front = c.front; back = c.back;
}

void Cuboid::print(){
	cout << "\n(l,t,r,b,f,bk) = (" << left << "," << top << "," << right << "," << bott << "," << front << "," << back << ")";
}

//@@@@@@@@@@@@@@@@@@@@ Ray @@@@@@@@@@@@@@@@@@@@@@@@@