//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: 
	mRGB* pixel; // array of pixels

public:
	int nRows, nCols; // dimensions of the pixmap
	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);
	}
	//<<<<<<<<<<<<<<<<< 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;
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
}


// 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 @@@@@@@@@@@@@@@@@@@@@@@@@