// file: app3.cc
//
// NOTE:
// I have put together this file for my Visualization Course, Fall 2001.
// It is mostly from F.S.Hill's book (appendix 3), plus some
// additional routines (e.g. writeBMPFile) and extensions (e.g. draw(x,y)).
// --Chee Yap (yap@cs.nyu.edu)
//
//Appendix 3. Some useful Classes and Utility Routines.
//This appendix defines some basic types and algorithms 
//that may prove useful as you develop graphics applications as 
//described throughout the book. The basic types are developed through 
//classes which can be used at will.

//Some of the classes are fully developed here. Others have a number
//of methods declared but not defined: It is up to the student to 
//flesh these out. Some classes are given here in very skeletal form,
//to suggest what might be developed in an actual application.

//In most classes a relaxed approach to encapsulation is taken:
//most data fields are declared public rather than private, as a
//matter of brevity and to avoid the need to define a large number
//of accessor and mutator functions.

//Another excellent source of classes and utilities may be found in
//the Graphics Gems series, whose on-line repository is
//http://www.acm.org/tog/GraphicsGems/index.html.

//Classes for 2D graphics.
// definition of simple support classes:
#include <string>
#include <iostream>
#include <fstream>
#include <strstream>
using namespace std;

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#ifdef _WINDOZ
	#include <windows.h> //change if using xWindows
	#include <gl/glut.h>
#else
	#include <GL/glut.h>
#endif

//@@@@@@@@@@@@@@@@@@ IntPoint class @@@@@@@@@@@@@@@@
class IntPoint{ // for 2D points with integer coordinates
public: 
	int x,y;
	void set(int dx, int dy){x = dx; y = dy;}
	void set(IntPoint& p){ x = p.x; y = p.y;} 
	IntPoint(int xx, int yy){x = xx; y = yy;}
	IntPoint(){ x = y = 0;}
};
//@@@@@@@@@@@@@@@@@@ Point2 class @@@@@@@@@@@@@@@@
class Point2{ // for 2D points with real coordinates
public: 
	float x,y;
	void set(float dx, float dy){x = dx; y = dy;}
	void set(Point2& p){ x = p.x; y = p.y;} 
	Point2(float xx, float yy){x = xx; y = yy;}
	Point2(){x = y = 0;}
};
//<<<<<<<<<<<<<<<<<<<<<< PolyLine >>>>>>>>>>>>>>>>>>>>>>>>>
class PolyLine{ // a polyline is a num plus an array of points
public:
	int num;
	Point2 pt[80]; //may need larger arrays in some circumstances
	PolyLine(){num = 0;}
};
// @@@@@@@@@@@@@@@@@@@@@@@@ IntRect class @@@@@@@@@@@@@@@@@@@@
class IntRect{ // a rectangle with integer border values
public: 
  int left, top, right, bott;
  IntRect(){left = top = right = bott = 0;}
  IntRect(int l, int t, int r, int b)
  {left = l; top = t; right = r; bott = b;}
  void set(int l, int t, int r, int b)
  {left = l; top = t; right = r; bott = b;}
  void set(IntRect& r)
  {left = r.left; top = r.top; right = r.right; bott = r.bott;}
};
//@@@@@@@@@@@@@@@@@@ Vector2 class @@@@@@@@@@@@@@@@
class Vector2{ 
public: 
	float x,y;
	void set(float dx, float dy){ x = dx; y = dy; } 
	void set(Vector2& v){ x = v.x; y = v.y;}
	void setDiff(Point2& a, Point2& b)//set to difference a - b
	{x = a.x - b.x; y = a.y - b.y;}
	void normalize()//adjust this vector to unit length
	{		double sizeSq = x * x + y * y;
		if(sizeSq < 0.0000001)
		{
			cerr << "\nnormalize() sees vector (0,0)!";
			return; // does nothing to zero vectors;
		}
		float scaleFactor = 1.0/(float)sqrt(sizeSq);
		x *= scaleFactor; y *= scaleFactor;
	}
	Vector2(float xx, float yy){x = xx; y = yy; }
	Vector2(Vector2& v){x = v.x; y = v.y; }
	Vector2(){x = y = 0;} //default constructor
	float dot(Vector2 b) // return this dotted with b
 	{return x * b.x + y * b.y;}
	void perp() // perp this vector
	{float tmp = x; x = -y; y = tmp;}
	float perpDot(Vector2& v) // return perp of this dotted with v
	{return x *v.x - y * v.y;}
};
//<<<<<<<<<<<<<<<<<<<< Canvas class >>>>>>>>>>>
// a global Canvas object (described in Chapter 3) knows how 
//to draw lines in world coordinates and to perform turtlegraphics
class Canvas {
 private:
	  Point2 CP; 	// current position in world
	  float CD; // current direction in degrees
 public:
	  float windowAspect;
	  Canvas(int width, int height, char* title)
	  {
		char* list; //dummy list for glutInit
		int numArgs = 1;//dummy value for glutInit
		glutInit(&numArgs, &list);
		glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
		glutInitWindowSize(width, height);
		glutInitWindowPosition(100, 100);
		glutCreateWindow(title);
		CP.x = CP.y = 0.0;
		windowAspect = 1.0;
	}
	void setWindow(float l, float r, float b, float t)
	{
		glMatrixMode(GL_PROJECTION);
		glLoadIdentity();
		gluOrtho2D((GLdouble)l, (GLdouble)r, (GLdouble)b, (GLdouble)t);
		if(t == b) return;
		windowAspect = (r - l)/(t - b);
	}
	void setViewport(int l, int r, int b, int t)
		{glViewport((GLint)l, (GLint)b, (GLint)(r-l), (GLint)(t-b));}	
	float getWindowAspect(void) { return windowAspect;}
	void lineTo(float x, float y)
	{
		glBegin(GL_LINES);
		glVertex2f((GLfloat)CP.x, (GLfloat)CP.y);
		CP.x = x; CP.y = y;
		glVertex2f((GLfloat)CP.x, (GLfloat)CP.y);
		glEnd(); 	glFlush();
	}
	void moveTo(float x, float y){CP.x = x; CP.y = y;}
	void turn(float ang) {CD += ang;}
	void turnTo(float ang) {CD = ang;}
	void forward(float dist, int vis)	
	{
		#define RadPerDeg 0.017453393 //radians per degree
		float x = CP.x + dist * cos(RadPerDeg * CD);
		float y = CP.y + dist * sin(RadPerDeg * CD);
		if(vis) lineTo(x, y);
		else moveTo(x, y);
		CP.x = x; CP.y = y;
	}
	void initCT() // initialize the CT (model view matrix)
	{
		glMatrixMode(GL_MODELVIEW); glLoadIdentity();
	}
	void rotate2D(double angle)
	{
		glMatrixMode(GL_MODELVIEW); glRotated(angle, 0.0, 0.0, 1.0);	
  }
	void translate2D(double dx, double dy)
	{
		glMatrixMode(GL_MODELVIEW); glTranslated(dx, dy, 0.0);
	}
	void scale2D(double sx, double sy)
	{
		glMatrixMode(GL_MODELVIEW); glScaled(sx, sy, 1.0);
	}
	void pushCT(void)
	{
		glMatrixMode(GL_MODELVIEW); glPushMatrix();
	}
	void popCT(void)
	{
		glMatrixMode(GL_MODELVIEW); glPopMatrix();
	}
	void ngon(int n, float cx, float cy, float radius);
};
//3D Classes for Graphics

//@@@@@@@@@@@@@@@@@@ Point3 class @@@@@@@@@@@@@@@@
class Point3{ 
public: 
	float x,y,z;
	void set(float dx, float dy, float dz){x = dx; y = dy; z = dz;}
	void set(Point3& p){x = p.x; y = p.y; z = p.z;} 	
 Point3(float xx,     float yy, float zz){x = xx; y = yy; z = zz;}
	Point3(){x = y = z = 0;}
	void build4tuple(float v[])
	{// load 4-tuple with this color: v[3] = 1 for homogeneous
		v[0] = x; v[1] = y; v[2] = z; v[3] = 1.0f;
	}	
};
//@@@@@@@@@@@@@@@@@@ Vector3 class @@@@@@@@@@@@@@@@
class Vector3{ 
public: 
	float x,y,z;
	void set(float dx, float dy, float dz){ x = dx; y = dy; z = dz;} 
	void set(Vector3& v){ x = v.x; y = v.y; z = v.z;}
	void flip(){x = -x; y = -y; z = -z;} // reverse this vector
	void setDiff(Point3& a, Point3& b)//set to difference a - b
	{ x = a.x - b.x; y = a.y - b.y; z = a.z - b.z;}
	void normalize()//adjust this vector to unit length
	{		
		double sizeSq = x * x + y * y + z * z;
		if(sizeSq < 0.0000001)
		{
			cerr << "\nnormalize() sees vector (0,0,0)!";
			return; // does nothing to zero vectors;
		}
		float scaleFactor = 1.0/(float)sqrt(sizeSq);
		x *= scaleFactor; y *= scaleFactor; z *= scaleFactor;
	}
	Vector3(float xx, float yy, float zz){x = xx; y = yy; z = zz;}
	Vector3(Vector3& v){x = v.x; y = v.y; z = v.z;}
	Vector3(){x = y = z = 0;} //default constructor
	Vector3 cross(Vector3 b) //return this cross b
	{
	   Vector3 c(y*b.z - z*b.y, z*b.x - x*b.z, x*b.y - y*b.x);
	   return c;
	}
	float dot(Vector3 b) // return this dotted with b
 {return x * b.x + y * b.y + z * b.z;}
};

// @@@@@@@@@@@@@@@@@@@@@ Color3 class @@@@@@@@@@@@@@@@
class Color3 { // holds an red,green,blue 3-tuple
public:
	float red, green, blue;
	Color3(){red = green = blue = 0;}
	Color3(float r, float g, float b){red = r; green = g; blue = b;}
	Color3(Color3& c){red = c.red; green = c.green; blue = c.blue;}
	void set(float r, float g, float b){red = r; green = g; blue = b;}
	void set(Color3& c) 
      	   {red = c.red; green = c.green; blue = c.blue;}
	void add(float r, float g, float b)
      	   {red += r; green += g; blue += b;}
	void add(Color3& src, Color3& refl)
	   { // add the product of source color and reflection coefficient
		red   += src.red   * refl.red;
		green += src.green * refl.green; 
		blue  += src.blue  * refl.blue;
	}
	void add(Color3& colr)
	{ // add colr to this color
	 red += colr.red ; green += colr.green; blue += colr.blue;}
	void build4tuple(float v[])
	{// load 4-tuple with this color: v[3] = 1 for homogeneous
		v[0] = red; v[1] = green; v[2] = blue; v[3] = 1.0f;
	}
};
//@@@@@@@@@@@@@@@@@@@@ light class @@@@@@@@@@@@@@@@@@@
class Light{ // for a linked list of light sources’ color and position
public:
	Point3 pos;
	Color3 color;
	Light* next;
	void setPosition(Point3 p){pos.set(p);}
	void setColor(Color3 c){color.set(c);}
	Light(){next = NULL;}
};


// RGBpixmap.h: a class to support working with RGB pixmaps.
#ifndef _RGBPIXMAP
#define _RGBPIXMAP
#include <fstream>

typedef unsigned char uchar;

class mRGB{   // the name RGB is already used by Windows
public: uchar r,g,b;
	mRGB(){r = g = b = 0;}
	mRGB(mRGB& p){r = p.r; g = p.g; b = p.b;}
	mRGB(uchar rr, uchar gg, uchar bb){r = rr; g = gg; b = bb;}
	void set(uchar rr, uchar gg, uchar bb){r = rr; g = gg; b = bb;}
};
//$$$$$$$$$$$$$$$$$ RGBPixmap class $$$$$$$$$$$$$$$
class RGBpixmap{
private: 
public:
	int nRows, nCols; // dimensions of the pixmap
	mRGB* pixel; // array of pixels

//METHODS:
	RGBpixmap() {nRows = nCols = 0; pixel = 0;}
	RGBpixmap(int rows, int cols) //constructor
	{
		nRows = rows;
		nCols = cols;
		pixel = new mRGB[rows*cols]; 
	}
	int readBMPFile(string fname); // read BMP file into this pixmap
	int writeBMPFile(char * fname); // write pixmap into BMP file
	void freeIt() // give back memory for this pixmap
	{
		delete []pixel; nRows = nCols = 0;
	}
	//<<<<<<<<<<<<<<<<<< copy >>>>>>>>>>>>>>>>>>>
	void copy(IntPoint from, IntPoint to, int x, int y, int width, int height)
	{ // copy a region of the display back onto the display
		if(nRows == 0 || nCols == 0) return;
		glCopyPixels(x, y, width, height,GL_COLOR);
	}
	//<<<<<<<<<<<<<<<<<<< draw >>>>>>>>>>>>>>>>>
	void draw()
	{ // draw this pixmap at current raster position
	if(nRows == 0 || nCols == 0) return;
		//tell OpenGL: don't align pixels to 4 byte boundaries in memory
		glPixelStorei(GL_UNPACK_ALIGNMENT,1);
		glDrawPixels(nCols, nRows,GL_RGB, GL_UNSIGNED_BYTE,pixel);
	}
	//<<<<<<<<<<<<<<<<<<< draw >>>>>>>>>>>>>>>>>
	void draw(int x, int y)
	{	// (x,y) is the origin of the translated image to be drawn
		// 4 cases, depending on which quadrant (x,y) is in 
		if (x >= 0 && y >= 0) {	// drag left and down
		  glRasterPos2i(0, 0);
		  glPixelStorei(GL_UNPACK_SKIP_ROWS, y);
		  glPixelStorei(GL_UNPACK_SKIP_PIXELS, x);
		  glDrawPixels(nCols - x, nRows - y, GL_RGB,
					GL_UNSIGNED_BYTE,pixel);
		}
		else if (x >= 0 && y < 0) { // drag left and up
		  glRasterPos2i(0, -y);
		  glPixelStorei(GL_UNPACK_SKIP_ROWS, 0);
		  glPixelStorei(GL_UNPACK_SKIP_PIXELS, x);
		  glDrawPixels(nCols - x, nRows + y, GL_RGB,
               		                GL_UNSIGNED_BYTE,pixel);
		}
		else if (x < 0 && y >= 0) { // drag right and down
		  glRasterPos2i(-x, 0);
		  glPixelStorei(GL_UNPACK_SKIP_ROWS, y);
		  glPixelStorei(GL_UNPACK_SKIP_PIXELS, 0);
		  glDrawPixels(nCols + x, nRows - y, GL_RGB,
               		                GL_UNSIGNED_BYTE,pixel);
		}
		else { // (x < 0 && y < 0), drag right and up
		  glRasterPos2i(-x, -y);
		  glPixelStorei(GL_UNPACK_SKIP_ROWS, 0);
		  glPixelStorei(GL_UNPACK_SKIP_PIXELS, 0);
		  glDrawPixels(nCols + x, nRows + y, GL_RGB,
               		                GL_UNSIGNED_BYTE,pixel);
		}
	} // draw(x, y)
	//<<<<<<<<<<<<<<<<< read >>>>>>>>>>>>>>>>
	int read(int x, int y, int wid, int ht)
	{ // read a rectangle of pixels into this pixmap
		nRows = ht;
		nCols = wid;
		pixel = new mRGB[nRows *nCols]; if(!pixel) return -1;

		//tell OpenGL: don't align pixels to 4 byte boundaries in memory
		glPixelStorei(GL_PACK_ALIGNMENT,1);
		glReadPixels(x, y, nCols, nRows, GL_RGB,GL_UNSIGNED_BYTE,pixel);
		return 0;
	}
	//<<<<<<<<<<<<<<<<< read from IntRect >>>>>>>>>>>>>>>>
	int read(IntRect r)
	{ // read a rectangle of pixels into this pixmap
		nRows = r.top - r.bott;
		nCols = r.right - r.left;
		pixel = new mRGB[nRows *nCols]; if(!pixel) return -1;
		//tell OpenGL: don’t align pixels to 4 byte boundaries in memory
		glPixelStorei(GL_PACK_ALIGNMENT,1);
		glReadPixels(r.left,r.bott, nCols, nRows, GL_RGB, GL_UNSIGNED_BYTE, pixel);
		return 0;
	}
	//<<<<<<<<<<<<<< setPixel >>>>>>>>>>>>>
	void setPixel(int x, int y, mRGB color)
	{
		if(x>=0 && x <nCols && y >=0 && y < nRows)
			pixel[nCols * y + x] = color;
	}
	//<<<<<<<<<<<<<<<< getPixel >>>>>>>>>>>
	mRGB getPixel(int x, int y)
	{ 		
		mRGB bad(255,255,255);
		assert(x >= 0 && x < nCols);
		assert(y >= 0 && y < nRows);
		return pixel[nCols * y + x];
	}
}; //end of class RGBpixmap
#endif

// RGBpixmap.cpp - routines to read a BMP file
// #include "RGBpixmap.h"
typedef unsigned short ushort;
typedef unsigned long ulong;
fstream inf; // global in this file for convenience
fstream outf; // output file
//<<<<<<<<<<<<<<<<<<<<< getShort >>>>>>>>>>>>>>>>>>>>
ushort getShort() //helper function
{ //BMP format uses little-endian integer types
  // get a 2-byte integer stored in little-endian form
		char ic;
		ushort ip;
		inf.get(ic); ip = ic;  //first byte is little one 
		inf.get(ic);  ip |= ((ushort)ic << 8); // or in high order byte
		return ip;
}
//<<<<<<<<<<<<<<<<<<<< getLong >>>>>>>>>>>>>>>>>>>
ulong getLong() //helper function
{  //BMP format uses little-endian integer types
   // get a 4-byte integer stored in little-endian form
		ulong ip = 0;
		char ic = 0;
		unsigned char uc = ic;
		inf.get(ic); uc = ic; ip = uc;
		inf.get(ic); uc = ic; ip |=((ulong)uc << 8);
		inf.get(ic); uc = ic; ip |=((ulong)uc << 16);
		inf.get(ic); uc = ic; ip |=((ulong)uc << 24);
		return ip;
	}
//<<<<<<<<<<<<<<<<<< RGBPixmap:: readBmpFile>>>>>>>>>>>>>
int RGBpixmap:: readBMPFile(string fname) 
{  // Read into memory an mRGB image from an uncompressed BMP file.
	// return 0 on failure, 1 on success
	inf.open(fname.c_str(), ios::in|ios::binary); //read binary char's
	if(!inf){ cout << " can't open file: " << fname << endl; return 0;}
	int k, row, col, numPadBytes, nBytesInRow;
	// read the file header information
	char ch1, ch2;
	inf.get(ch1); inf.get(ch2); 		//type: always 'BM'
	ulong fileSize =     	getLong();
	ushort reserved1 =   	getShort();     // always 0
	ushort reserved2= 	getShort();     // always 0 
	ulong offBits =		getLong();	// offset to image - unreliable
	ulong headerSize = 	getLong();   	// always 40
	ulong numCols =		getLong();	// number of columns in image
	ulong numRows = 	getLong();	// number of rows in image
	ushort planes=    	getShort();     // always 1 
	ushort bitsPerPixel=	getShort();   	//8 or 24; allow 24 here
	ulong compression =  	getLong();      // must be 0 for uncompressed 
	ulong imageSize = 	getLong();      // total bytes in image 
	ulong xPels =    	getLong();   	// always 0 
	ulong yPels =    	getLong();   	// always 0 
	ulong numLUTentries =	getLong();   	// 256 for 8 bit, otherwise 0 
	ulong impColors = 	getLong();      // always 0 
	if(bitsPerPixel != 24) 
	{ // error - must be a 24 bit uncompressed image
		cout << "not a 24 bit/pixelimage, or is compressed!\n";
		inf.close(); return 0;
	} 
	//add bytes at end of each row so total # is a multiple of 4
	// round up 3*numCols to next mult. of 4
	nBytesInRow = ((3 * numCols + 3)/4) * 4;
	numPadBytes = nBytesInRow - 3 * numCols; // need this many
	nRows = numRows; // set class's data members
	nCols = numCols;
   if (pixel) delete pixel;	// in case we read over an old image
   pixel = new mRGB[nRows * nCols]; //make space for array
	if(!pixel) return 0; // out of memory!
	long count = 0;
	char dum;
	for(row = 0; row < nRows; row++) // read pixel values
	{
		for(col = 0; col < nCols; col++)
		{
			char r,g,b;
			inf.get(b); inf.get(g); inf.get(r); //read bytes
			pixel[count].r = r; //place them in colors
			pixel[count].g = g;
			pixel[count++].b = b;
		}
   		for(k = 0; k < numPadBytes ; k++) //skip pad bytes at row's end
			inf >> dum;
	}
	inf.close(); return 1; // success
}// readBMPFile

// #include "writeBMP.cc"
// 	contributed by TingJen (yentj@cs.nyu.edu).
// 	It is slow because it writes one
//	byte at a time!  It uses putLong() and putShort() utility:

void putLong(unsigned long ip)
{
  outf.put((unsigned char)(ip & 0xff));
  outf.put((unsigned char)((ip>>8) & 0xff));
  outf.put((unsigned char)((ip>>16) & 0xff));
  outf.put((unsigned char)((ip>>24) & 0xff));
}

void putShort(unsigned short ip)
{
  outf.put((unsigned char)(ip & 0xff));
  outf.put((unsigned char)((ip>>8) & 0xff));
}

int RGBpixmap::writeBMPFile(char *fname)
{
  outf.open(fname, ios::out | ios::binary);
  if (!outf)
    {
      cerr << "can't open file: " << fname << endl;
      return -1;
    }

  outf.put('B');
  outf.put('M');  // file type  "BM"
  putLong(0x36L + nRows * nCols * 3);	//file size
  putLong(0L);				//reserved
  putLong(0x36L);			// offset to image
  putLong(40L);				//
  putLong((unsigned long) nCols);	//number of columns
  putLong((unsigned long) nRows);	//number of rows
  putShort(1);				//
  putShort(24);				// bits per pixel
  putLong(0L);				// not compressed
  putLong((unsigned long)nRows * nCols * 3);	// total bytes in image
  putLong(0L);				//
  putLong(0L);				//
  putLong(0L);				// numbers of colors : Not used
  putLong(0L);				//

  int count = 0;
  int numPadBytes = ((3 * nCols + 3)/4) * 4 - 3 * nCols;
  for (int y = 0; y < nRows; y ++)	// for each row
  {
    for (int x = 0; x < nCols ; x++)	// for each pixel
    {
      mRGB &p = pixel[count++];
      outf.put(p.b);
      outf.put(p.g);
      outf.put(p.r);
    }
    for (int k = 0; k < numPadBytes; k ++)	// padding
      outf.put(' ');
  }
  outf.close(); return 1;
}//writeBMPFile


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

class Noise{
public:
	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(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 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 marble(float strength,Point3& p)
	{
		float turbul = turbulence(10, p);
		float val = sin(6 * p.z + strength * turbul);
		return mySpline(val);
	}
private:
	float* noiseTable; // array of noise values
	unsigned char * index; //Pseudo random indices
	float 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 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)];
	}
}; // end of Noise class

// 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 @@@@@@@@@@@@@@@@
class Affine4 {// manages homogeneous affine transformations
// including inverse transformations
// and a stack to put them on
// used by Scene class to read SDL files
public:
	float m[16]; // hold a 4 by 4 matrix
	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 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 set(Affine4 a)// set this matrix to a
	{
		for(int i = 0; i < 16; i++)
			m[i]=a.m[i];
	}
	//<<<<<<<<<<<< postMult >>>>>>>>>>>
	void 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
	}
}; // end of Affine4 class

//@@@@@@@@@@ AffineNode class @@@@@@@@@@@
class AffineNode{
// used by Scene class to read SDL files
public:
	Affine4 * affn;
	Affine4 * invAffn;
	AffineNode * next;
	AffineNode()
	{
		next = NULL;
		affn = new Affine4; // new affine with identity in it
		invAffn = new Affine4; // and for the inverse
	}
	~AffineNode() //destructor
	{
		delete affn;
		delete invAffn;
	}
};
//@@@@@@@@@@@@@@@@ AffineStack class @@@@@@@@@@@@
class AffineStack{
// used by Scene class to read SDL files
public:
	AffineNode * tos;
	AffineStack()//default constructor;puts identity on top
	{
		tos = new AffineNode; // node with identity in it
		tos->next = NULL;
	}
	void dup()
	{
		AffineNode* tmp = new AffineNode;			
		tmp->affn = new Affine4(*(tos->affn));
		tmp->invAffn = new Affine4(*(tos->invAffn));
		tmp->next = tos;
		tos = tmp;
	}
	void setIdentity() // make top item the identity matrix
	{ 
		assert(tos != NULL);
		tos->affn->setIdentityMatrix();
		tos->invAffn->setIdentityMatrix();
	}
	void popAndDrop()
	{
		if(tos == NULL) return; // do nothing
		AffineNode *tmp = tos;
		tos = tos->next;
		delete tmp; // should call destructor, which deletes trices
	}	
	void releaseAffines()
	{ // pop and drop all remaining items 
	while(tos) popAndDrop();
	}
	void 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 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 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);
	}
}; // end of AffineStack class



//@@@@@@@@@@@@@@@@@@@ PointCluster class @@@@@@@@@@@@@@@@@@
class PointCluster{
public: // holds array of points for the bounding hull of a shape
	int num;
	Point3* pt;
	PointCluster() {num = 0; pt = NULL;}
	PointCluster(int n)// make a cluster of n points
	{
		pt = new Point3[n]; assert(pt);
		num = n;
  }
};
//@@@@@@@@@@@@@@@@@@@@@@@@@@ SphereInfo @@@@@@@@@@@@@@@@@@@@@
class SphereInfo{// holds the center and radius of a sphere
public:
	Point3 center;
	float radSq;
	void set(float x, float y, float z, float rsq)
	{
		center.set(x,y,z);
		radSq = rsq;
	}
};
//@@@@@@@@@@@@@@@@@@@@@@@@ Cuboid @@@@@@@@@@@@@@@@@@@@@@@@
class Cuboid{ // holds six border values of a cuboid
public:
	float left, top, right, bott, front, back;
	void 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 set(Cuboid& c)
	{
		left = c.left;top = c.top; right = c.right; bott = c.bott;
		front = c.front; back = c.back;
	}
	void print(){
		cout << "\n(l,t,r,b,f,bk) = (" << left << "," << top << "," << right << "," << bott << "," << front << "," << back << ")";
	}
};

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