// raytrace.h
// Ray Tracer Classes Module
// Copyright (C) 1997,1998,1999 Y.Nagatani

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <iomanip.h>
#include <math.h>

typedef	double	Real;
typedef	char	*String;
typedef	int	Bool;
#define	False	0
#define	True	1

void Error (String szMessage)
{
	cerr << szMessage;
	exit (1);
}

inline Real rRand (void) {
	// rand is 0 to (256*128)	  for System V or Solaris
	// rand is 0 to (256*256*256*128) for SunOs, BSD or Linux system
	return	((Real) rand()) / (RAND_MAX);
}

void Check_rRand (void)
{
	Real	A	= 0;
	int	nMax	= 100;
	for (int n = 0 ; n < nMax ; n++)
		A += rRand ();
	A /= nMax;
	if (0.3 < A && A < 0.7)
		return;		// randomness OK!
	
	cerr << "rRand: Bad Randomness. Average100 = " << A << endl;
	exit (1);
}


const Real	rSmallestValue = 0.00001;	// smallest real value
						// for zero Crossing check

// Classes defined in this file
class Vector;
class Line;
class Matrix;
class RGB;
class Object;
class Triangle;
class Spher;
class Node;
class BackGround;
class List;
class PostScript;

class Vector {
public:
	Real	x,	y,	z;
	
	Vector (void)			// default constructor
	{
		x = 0;
		y = 0;
		z = 0;
	}
	Vector (Real X, Real Y, Real Z)	// constructor
	{
		x = X;
		y = Y;
		z = Z;
	}
	Vector (Real s)			// constructor
	{
		x = s;
		y = s;
		z = s;
	}
	Vector (const Vector &v)	// copy constructor
	{
		x = v.x;
		y = v.y;
		z = v.z;
	}
	Vector & operator = (const Vector &v)	// substitution
	{
		if (this == &v)
			return	*this;
		
		x = v.x;
		y = v.y;
		z = v.z;
		
		return	*this;
	}
	~Vector (void)			// destructor
	{
		// No Function
	}
	Vector & set (Real X = 0, Real Y = 0, Real Z = 0)
	{
		x = X;
		y = Y;
		z = Z;
		return	*this;
	}
	Vector & operator += (const Vector &v)
	{
		x += v.x;
		y += v.y;
		z += v.z;
		return	*this;
	}
	Vector & operator -= (const Vector &v)
	{
		x -= v.x;
		y -= v.y;
		z -= v.z;
		return	*this;
	}
	Vector & operator *= (const Real s)
	{
		x *= s;
		y *= s;
		z *= s;
		return	*this;
	}
	Vector & operator /= (const Real s)
	{
		x /= s;
		y /= s;
		z /= s;
		return	*this;
	}
	Vector operator + (const Vector &v) const
	{
		return	Vector (x+v.x, y+v.y, z+v.z);
	}
	Vector operator - (const Vector &v) const
	{
		return	Vector (x-v.x, y-v.y, z-v.z);
	}
	Vector operator * (const Real s) const
	{
		return	Vector (x*s, y*s, z*s);
	}
	Vector operator / (const Real s) const
	{
		return	Vector (x/s, y/s, z/s);
	}
	friend Vector operator * (const Real s, const Vector &v)
	{
		return	Vector (s*v.x, s*v.y, s*v.z);
	}
	Vector out (const Vector &v) const	// Outer product
	{
		Vector	w (y * v.z - z * v.y,
			   z * v.x - x * v.z,
			   x * v.y - y * v.x);
		return	w;
	}
	Real in (const Vector &v) const		// Inner product
	{
		return	x*v.x + y*v.y + z*v.z;
	}
	Real nin (const Vector &v) const	// Normalized Inner product
	{
		return	in (v) / ( norm () * v.norm () );
	}
	Real norm (void) const
	{
		return	sqrt (x*x + y*y + z*z);
	}
	Vector unit (void) const		// unit direction
	{
		return	Vector (*this / norm ());
	}
	Vector & unitize (void)			// unitize
	{
		*this /= norm ();
		return	*this;
	}
	Vector GetOrthogonalOne (void) const	// Get One Orthogonal Unit Vec
	{
		Vector	vA (unit ());			// unit this Vector
		Vector	v1 (1,0,0);			// some any unit vecor
		Real	A1 = vA .in (v1);		// vA element of v1 vec
		if (fabs(A1) >= 1.0-rSmallestValue) {	// if vA // v1
			v1 = Vector (0,1,0);		// another unit vecor
			A1 = vA .in (v1);
		}
		v1 = v1 - A1 * vA;			// v1 otrho to vA
		return	v1.unit ();
	}
	Vector GetRotateAxesLine (const Line &L, Real theta) const;
	Vector RotateAxesLine (const Line &L, Real theta)
	{
		*this = (*this).GetRotateAxesLine (L, theta);
		return	*this;
	}
	
	friend ostream & operator << (ostream& s, const Vector &v)
	{
		s << '['
		  << v.x << ','
		  << v.y << ','
		  << v.z << ']';
		return	s;
	}
};


class Line {
public:
	Vector	v,	w;
	Line (void)				// default constructor
	{
		v = 0;
		w = 0;
	}
	Line (const Vector v1, const Vector v2)	// constructor
	{
		v = v1;
		w = v2;
	}
	Line (Real x1, Real y1, Real z1,
	      Real x2, Real y2, Real z2)	// constructor
	{
		v.set (x1,y1,z1);
		w.set (x2,y2,z2);
	}
	Line (const Line &l)			// copy constructor
	{
		v = l.v;
		w = l.w;
	}
	Line & operator = (const Line &l)	// substitution
	{
		if (this == &l)
			return	*this;
		v = l.v;
		w = l.w;
		return	*this;
	}
	~Line (void)				// destructor
	{
		// No Function
	}
	Vector dir (void) const 	// line direction
	{
		return	Vector (w - v);
	}
	Vector ndir (void) const 	// line normalized direction
	{
		return	dir() / length ();
	}
	Real length (void) const 	// line length
	{
		return	dir () . norm ();
	}
	Real rCross (const Line &L) const
	{
		return	(dir () .out (L.v - v)) .in (L.dir ());
	}
	Bool operator > (const Line &L) const
	{
		return	(rCross (L) > 0.0);
	}
	Bool operator < (const Line &L) const
	{
		return	(rCross (L) < 0.0);
	}
	Bool operator == (const Line &L) const
	{
		return	(rCross (L) == 0.0);
	}
	friend ostream & operator << (ostream& s, const Line &L)
	{
		s << '['
		  << L.v << ','
		  << L.w << ','
		  << L.dir () << ']';
		return	s;
	}
	
};


Vector Vector::GetRotateAxesLine (const Line &L, Real theta) const
	{
		Vector	RotCenter = L.v;
		Vector	RotAxes   = L.ndir (),			// unit Axes
			RotN1 (RotAxes .GetOrthogonalOne ()),	// unit ortho
			RotN2 (RotAxes .out (RotN1));
		// (RotAxes, RotN1, RotN2) Complete Ortho Normal System
		
		Vector	V =  *this - RotCenter;
		
		Real	Vr = V .in (RotAxes),	// elements of CONS
			V1 = V .in (RotN1),
			V2 = V .in (RotN2);
		
		Real	W1 = V1*cos(theta) - V2*sin(theta),
			W2 = V1*sin(theta) + V2*cos(theta);
		
		Vector	W (Vr*RotAxes + W1*RotN1 + W2*RotN2);
		return	W + RotCenter;
	}


class Matrix {
public:
	Real	m[3][3];
	Matrix (void)				// default constructor
	{
		for (int i = 0 ; i < 3 ; i++)
			for (int j = 0 ; j < 3 ; j++)
				m[i][j] = 0;
	}
	Matrix (const Vector v0,
		const Vector v1,
		const Vector v2)	// constructor
	{
		m[0][0] = v0.x;	m[0][1] = v1.x;	m[0][2] = v2.x;
		m[1][0] = v0.y;	m[1][1] = v1.y;	m[1][2] = v2.y;
		m[2][0] = v0.z;	m[2][1] = v1.z;	m[2][2] = v2.z;
	}
	Matrix (Real A, Real B, Real C,
		Real D, Real E, Real F,
		Real G, Real H, Real I)	// constructor
	{
		m[0][0] = A;	m[0][1] = B;	m[0][2] = C;
		m[1][0] = D;	m[1][1] = E;	m[1][2] = F;
		m[2][0] = G;	m[2][1] = H;	m[2][2] = I;
	}
	Matrix (const Matrix &M)		// copy constructor
	{
		for (int i = 0 ; i < 3 ; i++)
			for (int j = 0 ; j < 3 ; j++)
				m[i][j] = M.m[i][j];
	}
	Matrix & operator = (const Matrix &M)	// substitution
	{
		if (this == &M)
			return	*this;
		for (int i = 0 ; i < 3 ; i++)
			for (int j = 0 ; j < 3 ; j++)
				m[i][j] = M.m[i][j];
		return	*this;
	}
	~Matrix (void)				// destructor
	{
		// No Function
	}
	Matrix & operator *= (const Matrix &M)
	{
		*this = *this * M;
		return	*this;
	}
	Matrix operator * (const Matrix &M) const	// Matrix*Matrix 
	{
		Matrix	N;
		for (int i = 0 ; i < 3 ; i++)
			for (int j = 0 ; j < 3 ; j++)
				for (int k = 0 ; k < 3 ; k++)
					N.m[i][j] += m[i][k] * M.m[k][j];
		return	N;
	}
	Vector operator * (const Vector &v) const	// Matrix*Vector
	{
		Vector	V;
		V.x = m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z;
		V.y = m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z;
		V.z = m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z;
		return	V;
	}
	Real det (void) const				// determinant
	{
		return	+ m[0][0]*m[1][1]*m[2][2]
			- m[0][0]*m[1][2]*m[2][1]
			+ m[0][1]*m[1][2]*m[2][0]
			- m[0][1]*m[1][0]*m[2][2]
			+ m[0][2]*m[1][0]*m[2][1]
			- m[0][2]*m[1][1]*m[2][0];
	}
	Real co (int i, int j) const		// co-matrix elements
	{
		Real	d[2][2];
		
		for (int k = 0 ; k < 3 ; k++) {
			if (k == i)
				continue;
			for (int l = 0 ; l < 3 ; l++) {
				if (l == j)
					continue;
				int	n0,	m0;
				n0 = (k < i ? k : k-1);
				m0 = (l < j ? l : l-1);
				d[n0][m0] = m[k][l];
			}
		}
		return	((i+j) % 2 == 0 ? 1.0 : -1.0)
			* (d[0][0]*d[1][1] - d[0][1]*d[1][0]);
	}
	Matrix tco (void) const			// transposed co-matrix
	{
		Matrix	N;
		for (int i = 0 ; i < 3 ; i++)
			for (int j = 0 ; j < 3 ; j++)
				N.m[i][j] = co (j,i);
		return	N;
	}
	
	friend ostream & operator << (ostream& s, const Matrix &M)
	{
		s << '[';
		for (int i = 0 ; i < 3 ; i++) {
			s << '[';
			for (int j = 0 ; j < 3 ; j++) {
				s << M.m[i][j];
				if (j < 2)
					s << ',';
			}
			s << ']';
			if (i < 2)
				s << ',';
		}
		s << ']';
		return	s;
	}
};


class RGB {
public:
	Real	rRed,
		rGreen,
		rBlue;
	RGB (void)
	{
		rRed	= 0;
		rGreen	= 0;
		rBlue	= 0;
	}
	RGB (Real R, Real G, Real B)
	{
		rRed	= R;
		rGreen	= G;
		rBlue	= B;
	}
	RGB (Real a)
	{
		rRed	= a;	// Gray Color Set
		rGreen	= a;
		rBlue	= a;
	}
	RGB (const RGB &C)	// copy constructor
	{
		rRed	= C.rRed;
		rGreen	= C.rGreen;
		rBlue	= C.rBlue;		
	}
	RGB & operator = (const RGB &C)		// substitution
	{
		if (this == &C)
			return	*this;
		rRed	= C.rRed;
		rGreen	= C.rGreen;
		rBlue	= C.rBlue;		
		return	*this;
	}
	~RGB (void)				// destructor
	{
		// No Function
	}
	RGB & operator += (const RGB &C)	// RGB += RGB
	{
		rRed	+= C.rRed;
		rGreen	+= C.rGreen;
		rBlue	+= C.rBlue;
		return	*this;
	}
	RGB & operator *= (const RGB &C)	// RGB *= RGB
	{					// filtering
		rRed	*= C.rRed;
		rGreen	*= C.rGreen;
		rBlue	*= C.rBlue;
		return	*this;
	}
	friend RGB operator + (const RGB &C1, const RGB &C2)	// RGB + RGB
	{							// add
		RGB	D;
		
		D.rRed	 = C1.rRed   + C2.rRed;
		D.rGreen = C1.rGreen + C2.rGreen;
		D.rBlue	 = C1.rBlue  + C2.rBlue;
		
		return	D;
	}
	friend RGB operator * (const RGB &C1, const RGB &C2)	// RGB * RGB
	{							// filtering
		RGB	D;
		
		D.rRed	 = C1.rRed   * C2.rRed;
		D.rGreen = C1.rGreen * C2.rGreen;
		D.rBlue	 = C1.rBlue  * C2.rBlue;
		
		return	D;
	}
	friend ostream & operator << (ostream& s, const RGB C)
	{
		s << '['
		  << "R:" << C.rRed   << ','
		  << "G:" << C.rGreen << ','
		  << "B:" << C.rBlue  << ']';
		return	s;
	}
	RGB & cutin (void)	// cut in 0.0 to 1.0
	{
		rRed	= cutin (rRed);
		rGreen	= cutin (rGreen);
		rBlue	= cutin (rBlue);
		return	*this;
	}
private:
	Real cutin (Real a)
	{
		a = (a > 1.0 ? 1.0 : a);
		a = (a < 0.0 ? 0.0 : a);
		return	a;
	}
};


class Object {
private:
	void ErrorUndefVirtual (String name) const
	{
		cerr << "Call undefined virtual function ("
		     << name << ") in class (Object)" << endl;
		exit (1);
	}
public:
	// New : memory allocation
	virtual Object * copyNew (void) const // = 0;
	{
		ErrorUndefVirtual ("copyNew");
		return	new Object (*this);
	}
	
	// Normal direction on vPoint 
	virtual Vector normal (const Vector &vPoint) const // = 0;
	{
		ErrorUndefVirtual ("normal");
		return	0;
	}
	
	// Normalized normal direction on on vPoint
	virtual Vector nnormal (const Vector &vPoint) const // = 0;
	{
		ErrorUndefVirtual ("nnormal");
		return	0;
	}
	
	// Get Cross Point on Line L as vCross,
	//     Distacne from L.v as rDistance
	// and return Cross or not as Bool.
	virtual Bool CrossPoint (const Line &L,
				 Vector &vCross,
				 Real &rDistance) const // = 0;
	{
		ErrorUndefVirtual ("CrossPoint");
		return	False;
	}
	
	virtual Object & operator += (const Vector &v)
	{
		ErrorUndefVirtual ("+=");
		return	*this;
	}
	virtual Object & operator -= (const Vector &v)
	{
		ErrorUndefVirtual ("-=");
		return	*this;
	}
	
	virtual Object & RotateAxesLine (const Line &L, Real theta)
	{
		ErrorUndefVirtual ("RotateAxesLine");
		return	*this;
	}
	
	
	// Color Propaties
	RGB	rgbSelf;	// Object Self Shining Color
	
				// Object Reflection Properties:
	RGB	rgbRefl;	//  Color Filtering Ratio
	Bool	bReflBlur;	//  Reflection Bluring (Bool)
	int	nReflBlur;	//  Nubmer of sum directons for Bluring
	Real	rReflBlurAngle;	//  Integral Angle (radian)
				//	0	: point
				//	M_PI	: half angle integral
				//	2*M_PI	: total angle integral
	
				// Total Direction Reflection Properties:
	RGB	rgbTotalBlur;	//  Color Filtering Ratio
	Bool	bTotalBlur;	//  Bluring (Bool)
	int	nTotalBlur;	//  Nubmer of sum directons for Bluring
	
	RGB	rgbTrans;	// Object Transparent Color Filtering Ratio
	
	Real	rRefraction;	// Refraction Ratio	(* 屈折率 *)
	
	
	Object & clear (void)
	{
		rgbSelf		= RGB (0,0,0);
		
		rgbRefl		= RGB (0,0,0);
		bReflBlur	= False;
		nReflBlur	= 0;
		rReflBlurAngle	= 0;
		
		rgbTotalBlur	= RGB (0,0,0);
		bTotalBlur	= False;
		nTotalBlur	= 0;
		
		rgbTrans	= RGB (0,0,0);
		
		rRefraction	= 1.0;
		
		return	*this;
	}
	
	Object (RGB Self=0 , RGB Refl=0, RGB Trans=0)
	{
		clear ();
		
		rgbSelf		= Self;
		rgbRefl		= Refl;
		rgbTrans	= Trans;
	}
	
	virtual ~Object (void)			// destructor
	{
		// No Function
	}
	
	Object & setRGB (RGB Self=0 , RGB Refl=0, RGB Trans=0)
	{
		rgbSelf		= Self;
		rgbRefl		= Refl;
		rgbTrans	= Trans;
		
		return	*this;
	}
	virtual RGB getRGBSelf (Vector vPoint)
	{
		return	rgbSelf;
	}
	virtual RGB getRGBRefl (Vector vPoint)
	{
		return	rgbRefl;
	}
	virtual RGB getRGBTrans (Vector vPoint)
	{
		return	rgbTrans;
	}
	Object & setRefraction (Real rRR = 1.0)
	{
		rRefraction	= rRR;
		
		return	*this;
	}
	Object & setReflBlur (int nBlur = 10, Real rAngle = 0.05)
	{
		bReflBlur	= True;
		// rgbRefl	= rgb;
		nReflBlur	= nBlur;
		rReflBlurAngle	= rAngle;
		
		return	*this;
	}
	Object & setTotalBlur (RGB rgb = 0.5, int nBlur = 10)
	{
		bTotalBlur	= True;
		rgbTotalBlur	= rgb;
		nTotalBlur	= nBlur;
		
		return	*this;
	}
	virtual ostream& print (ostream& s) const
	{
		s << '['
		  << "Self:"		<< rgbSelf		<< ','
		  << "Refl:"		<< rgbRefl		<< ','
		  << "Trans:"		<< rgbTrans		<< ",\n";
		
		if (bReflBlur)
			s << "nReflBlur:"	<< nReflBlur		<< ','
			  << "rReflBlurAngle:"	<< rReflBlurAngle	<< ',';
		else
			s << "(NoBlurRef),";
		
		if (bTotalBlur)
			s << "rgbTotalBlur:"	<< rgbTotalBlur		<< ','
			  << "nTotalBlur:"	<< nTotalBlur		<< ',';
		else
			s << "(NoBlur),";
		
		s << "Refraction:"	<< rRefraction;
		
		s << ']';
		return	s;
	}
};

ostream & operator << (ostream& s, const Object &O)
{
	O.print (s);
	return	s;
}


class Triangle : public Object {
public:
	Vector	a,	b,	c;
	Triangle (void)					// default constructor
		: Object ()
	{
		a = 0;
		b = 0;
		c = 0;
	}
	Triangle (const Vector v1,
		  const Vector v2,
		  const Vector v3)			// constructor
		: Object ()
	{
		a = v1;
		b = v2;
		c = v3;
	}
	Triangle (Real x1, Real y1, Real z1,
		  Real x2, Real y2, Real z2,
		  Real x3, Real y3, Real z3)		// constructor
		: Object ()
	{
		a.set (x1,y1,z1);
		b.set (x2,y2,z2);
		c.set (x3,y3,z3);
	}
	Triangle (const Triangle &T)			// copy constructor
		: Object (T)
	{
		a = T.a;
		b = T.b;
		c = T.c;
	}
	Triangle & operator = (const Triangle &T)	// substitution
	{
		if (this == &T)
			return	*this;
		
		Object::operator = (T);
		
		a = T.a;
		b = T.b;
		c = T.c;
		return	*this;
	}
	virtual ~Triangle (void)			// destructor
	{
		// No Function
	}
	virtual Object & operator += (const Vector &v)
	{
		a += v;
		b += v;
		c += v;
		return	*this;
	}
	virtual Object & operator -= (const Vector &v)
	{
		a -= v;
		b -= v;
		c -= v;
		return	*this;
	}
	Triangle operator + (const Vector &v) const
	{
		Triangle	T (*this);
		T.a += v;
		T.b += v;
		T.c += v;
		return		T;
	}
	Triangle operator - (const Vector &v) const
	{
		Triangle	T (*this);
		T.a -= v;
		T.b -= v;
		T.c -= v;
		return		T;
	}
	
	// Object Normal Direction at vPoint.
	// But in Triangle, we can neglect vPoint
	virtual Vector normal (void) const
	{
		return	Vector ((b-a). out (c-a));
	}
	virtual Vector nnormal (void) const			// Normalized
	{
		return	normal () .unit ();
	}
	virtual Vector normal (const Vector &vPoint) const
	{
		return	normal();
	}
	virtual Vector nnormal (const Vector &vPoint) const	// Normalized
	{
		return	normal () .unit ();
	}
	
	
	virtual Bool CrossPoint (const Line &L,
				 Vector &vCross,
				 Real &rDistance) const
	{
		Matrix	M (b-a, c-a, L.dir () * (-1));
		Real	det = M.det ();
		if (det == 0.0)
			return	False;	// Triangle degenerate or
					// Triangle and Line are Pararell
		
		Vector	vParam = M.tco () * (L.v - a) / det;
		if (vParam.x < 0 || vParam.y < 0 ||
		    vParam.x + vParam.y > 1)
			return	False;	// Line Cross Plane
					// but not Cross Triangle
		
		if (vParam.z <= rSmallestValue)
			return	False;	// Line Cross Triangle
					// but negative-Direction of Line
					// or 0 point of Line
		
		vCross = L.v + vParam.z * L.dir ();
		rDistance = vParam.z * L.length ();
		return	True;
	}
	
	virtual ostream& print (ostream& s) const
	{
		s << '[';
		Object::print(s) << '\n';
		s << '['
		  << a << ','
		  << b << ','
		  << c << ';'
		  << normal () << ']';
		s << ']';
		return	s;
	}
	
	virtual Object * copyNew (void) const
	{
		return	new Triangle (*this);
	}
	
	virtual Object & RotateAxesLine (const Line &L, Real theta)
	{
		a.RotateAxesLine (L,theta);
		b.RotateAxesLine (L,theta);
		c.RotateAxesLine (L,theta);
		return	*this;
	}
	Triangle GetRotateAxesLine (const Line &L, Real theta) const
	{
		Triangle	T (*this);
		T.RotateAxesLine (L, theta);
		return	T;
	}
	void getBaseCoefs (Vector vPoint, Real * Pe1, Real * Pe2)
	{
		// vPoint must be on Triangle Plane
		
		Vector	vPos = vPoint - a;	// Position Vector
		Vector	vE1  = b      - a,
			vE2  = c      - a;	// Base Vector
		
		Real	E11 = vE1.in (vE1),
			E12 = vE1.in (vE2),
			E21 = E12,
			E22 = vE2.in(vE2);		// Metric
		
		Real	detE = E11*E22 - E21*E12;	// Det Metric
		
		Real	e11 =   E22 / detE,
			e12 = - E12 / detE,
			e21 = - E21 / detE,
			e22 =   E11 / detE;		// Inv Metric
		
		Real	PE1 = vPos.in (vE1),
			PE2 = vPos.in (vE2);
		
		*Pe1 = e11*PE1 + e12*PE2;
		*Pe2 = e21*PE1 + e22*PE2;
	}
};


class Quadrilateral : public Triangle {
public:
	Quadrilateral (void)				// default constructor
		: Triangle ()
	{
	}
	Quadrilateral (const Vector v1,
		  const Vector v2,
		  const Vector v3)			// constructor
		: Triangle (v1, v2, v3)
	{
	}
	Quadrilateral (Real x1, Real y1, Real z1,
		  Real x2, Real y2, Real z2,
		  Real x3, Real y3, Real z3)		// constructor
		: Triangle (x1,y1,z1, x2,y2,z2, x3,y3,z3)
	{
	}
	Quadrilateral (const Quadrilateral &T)		// copy constructor
		: Triangle (T)
	{
	}
	virtual Object * copyNew (void) const
	{
		return	new Quadrilateral (*this);
	}
	virtual Bool CrossPoint (const Line &L,
				 Vector &vCross,
				 Real &rDistance) const
	{
		Matrix	M (b-a, c-a, L.dir () * (-1));
		Real	det = M.det ();
		if (det == 0.0)
			return	False;	// Quadrilateral degenerate or
					// Quadrilateral and Line are Pararell
		
		Vector	vParam = M.tco () * (L.v - a) / det;
		if (vParam.x < 0 || vParam.y < 0 ||
		    vParam.x > 1 || vParam.y > 1)
			return	False;	// Line Cross Plane
					// but not Cross Area
		
		if (vParam.z <= rSmallestValue)
			return	False;	// Line Cross Area
					// but negative-Direction of Line
					// or 0 point of Line
		
		vCross = L.v + vParam.z * L.dir ();
		rDistance = vParam.z * L.length ();
		return	True;
	}
};


class Spher : public Object {
public:
	Vector	c;	// center position
	Real	r;	// radius
	
	Spher (void)					// default constructor
		: Object ()
	{
		c = 0;
		r = 0;
	}
	Spher (const Vector vCenter,
		  const Real rRadius)			// constructor
		: Object ()
	{
		c = vCenter;
		r = rRadius;
	}
	Spher (Real x, Real y, Real z, Real rRadius)	// constructor
		: Object ()
	{
		c.set (x,y,z);
		r = rRadius;
	}
	Spher (const Spher &S)				// copy constructor
		: Object (S)
	{
		c = S.c;
		r = S.r;
	}
	Spher & operator = (const Spher &S)		// substitution
	{
		if (this == &S)
			return	*this;
		
		Object::operator = (S);
		
		c = S.c;
		r = S.r;
		return	*this;
	}
	virtual ~Spher (void)			// destructor
	{
		// No Function
	}
	virtual Object & operator += (const Vector &v)
	{
		c += v;
		return	*this;
	}
	virtual Object & operator -= (const Vector &v)
	{
		c -= v;
		return	*this;
	}
	Spher operator + (const Vector &v) const
	{
		Spher	S (*this);
		S.c += v;
		return		S;
	}
	Spher operator - (const Vector &v) const
	{
		Spher	S (*this);
		S.c -= v;
		return		S;
	}
	
	// Object Normal Direction at vPoint.
	virtual Vector normal (const Vector &vPoint) const
	{
		return	vPoint - c;
	}
	virtual Vector nnormal (const Vector &vPoint) const	// Normalized
	{
		return	normal (vPoint) .unit ();
	}
	
	virtual Bool CrossPoint (const Line &L,
				 Vector &vCross,
				 Real &rDistance) const
	{
		Vector	vDir = L.ndir ();	// normalized ray direction;
		Vector	vS   = c - L.v;		// Spher Center direction
		
		Real	DS = vDir .in (vS);	// inner product
		Real	rS  = vS.norm ();	// distance to Center
		
		Real	rCheck = DS*DS + r*r - rS*rS;
		
		if (rCheck < 0)
			return	False;	// Line does not cross to Spher
		
		Real	tMax = DS + sqrt (rCheck),	// 2 solutions
			tMin = DS - sqrt (rCheck);
		
		if (tMax <= rSmallestValue)
			return	False;	// Line Cross Spher
					// but negative-Direction of Line
					// or 0 point of Line
		Real	t = (tMin <= rSmallestValue ? tMax : tMin);
		
		vCross = L.v + t * vDir;
		rDistance = t;
		return	True;
	}
	
	virtual ostream& print (ostream& s) const
	{
		s << '[';
		Object::print(s) << '\n';
		s << '['
		  << "Cent:" << c << ','
		  << "Radi:" << r << ']';
		s << ']';
		return	s;
	}
	
	virtual Spher * copyNew (void) const
	{
		return	new Spher (*this);
	}
	
	virtual Object & RotateAxesLine (const Line &L, Real theta)
	{
		c.RotateAxesLine (L,theta);
		return	*this;
	}
	Spher GetRotateAxesLine (const Line &L, Real theta) const
	{
		Spher	S (*this);
		S.RotateAxesLine (L, theta);
		return	S;
	}
};


class BackGround {
private:
	RGB	rgbBG;
	
	Vector	vLightDir;
	RGB	rgbDir;
	Real	nSpotShape;
public:	
	BackGround (void)
	{
		rgbBG	 = RGB (0,0,0);
		
		vLightDir = Vector (0,0,1);	// must be non-zero vec.
		rgbDir	  = RGB (0,0,0);
		nSpotShape = 2;
	}
	RGB getRGB (const Vector &vDir) const
	{
		// normalized inner product
		Real	a = - (vLightDir .nin (vDir));
		if (a < 0)
			a = 0;
		
		// b = a ^ nSpotShape;
		Real	b = (a > 0 ? exp (nSpotShape * log(a)) : 0);
		
		return	rgbBG + b*rgbDir;
	}
	BackGround & setRGB (const RGB &rgb)
	{
		rgbBG = rgb;
		return	*this;
	}
	BackGround & setLightDir (const Vector &V,
				  const RGB &rgb,
				  int nShape = 2)
	{
		vLightDir = V;
		rgbDir	  = rgb;
		nSpotShape = nShape;
		return	*this;
	}
	BackGround & RotateAxesLine (const Line &L, Real theta)
	{
		vLightDir.RotateAxesLine (Line (0, L.dir ()), theta);
		return	*this;
	}
};


class Node {
	friend class List;
	friend ostream & operator << (ostream& s, const List &L);
private:
	Node	*next;
	Object	*O;
public:
	Node (const Object &toNew, Node *node = NULL)
	{
		next  = node;
		O = toNew.copyNew ();
//		cerr << "ORG:" << toNew	<< endl;
//		cerr << "NEW:" << *O	<< endl;
	}
};


class List {
private:
	Node		*root;
	BackGround	BG;
public:
	List (void)
	{
		root = NULL;
	}
	List (const List &L)
	{
		root	= NULL;
		BG	= L.BG;
		
		for (Node *node = L.root ; node ; node = node->next)
			add (*(node->O));
	}
	List & add (const Object &toNew)
	{
		Node	*NewNode = new Node (toNew, root);
		if (NewNode == NULL)
			Error ("Memery Shotage in Add List.\n");
		
		root = NewNode;
		return	*this;
	}
	List & add (const List &L)
	{
		for (Node *node = L.root ; node ; node = node->next)
			add (*(node->O));
		
		return	*this;
	}

	List & addTail (const Object &toNew)
	{
		Node	*NewNode = new Node (toNew, NULL);
		if (NewNode == NULL)
			Error ("Memery Shotage in Add List to Tail.\n");
		
		if (root != NULL) {
			Node	*node = root;
			for (; node->next ; node = node->next);
			node->next = NewNode;
		} else {
			root = NewNode;
		}
		return	*this;
	}
	List & addTail (const List &L)
	{
		for (Node *node = L.root ; node ; node = node->next)
			addTail (*(node->O));
		
		return	*this;
	}
	List & operator += (const Vector &V)
	{
		for (Node *node = root ; node ; node = node->next)
			*(node->O) += V;
		
		return	*this;
	}
	List & operator -= (const Vector &V)
	{
		for (Node *node = root ; node ; node = node->next)
			*(node->O) -= V;
		
		return	*this;
	}
	List & RotateAxesLine (const Line &L, Real theta)
	{
		BG.RotateAxesLine (Line (0, L.dir ()), theta);
		
		for (Node *node = root ; node ; node = node->next)
			node->O->RotateAxesLine (L, theta);
		
		return	*this;
	}
	
	int nSize (void) const
	{
		int	n = 0;
		for (Node *node = root ; node ; node = node->next)
			n++;
		return	n;
	}
	friend ostream & operator << (ostream& s, const List &L)
	{
		s << '[';
		for (const Node *node = L.root ; node ; node = node->next)
			s << *(node->O) << ',' << '\n';
		s << ']';
		return	s;
	}
	List & setBackGround (const RGB &bg)
	{
		BG.setRGB (bg);
		return	*this;
	}
	List & setBackGroundDir (const Vector &vDir,
				 const RGB &rgbLight,
				 int nSpotShape = 2)
	{
		BG.setLightDir (vDir, rgbLight, nSpotShape);
		return	*this;
	}
	// get Most Neighbourfood Cross Point, Distance and TriObj.
	Bool CrossPoint (const Line &L,
			 Vector	&vCross,	// cross point
			 Real	&rDistance,	// distance
			 Object*&objCrossed) const
	{
		Bool	bCrossed = False;
		Real	rMinimumDist;
		Vector	vMinimumCross;
		const Node	*node,	*nodeCross;
		for (node = root ; node ; node = node->next) {
			
			Vector	vNewCross;
			Real	rNewDist;
			
			if (!node->O->CrossPoint (L, vNewCross, rNewDist))
				continue;
			
			if (!bCrossed) {
				vMinimumCross = vNewCross;
				rMinimumDist  = rNewDist;
				nodeCross     = node;
				bCrossed = True;
			}
			if (rMinimumDist > rNewDist) {
				vMinimumCross = vNewCross;
				rMinimumDist  = rNewDist;
				nodeCross     = node;
			}
		}
		if (bCrossed) {
			vCross = vMinimumCross;
			rDistance = rMinimumDist;
			objCrossed = nodeCross->O;
			return	True;
		}
		return	False;
	}
	// get Most Neighbourfood Cross Point, Distance and TriObj.
	RGB GetRGB (int nLevel, const Line &L) const
	{
		// Cross to No Objects then return BackGround
		Vector	vCross;
		Real	rDistance;
		Object	*pobjCrossed;
		RGB	rgbBG	= BG.getRGB (L.dir ());
		if (!CrossPoint (L, vCross, rDistance, pobjCrossed))
			return	rgbBG;
		
		Vector	ViewRay		= L.ndir (),
			NormalVec	= pobjCrossed->nnormal (vCross);
		Real	RN		= ViewRay. in (NormalVec);
		
		// Common Color Buffer
		RGB	rgbNew (0,0,0);
		
		// Self Shining Color
		RGB	rgbSelf = (RN > 0 ?
			    pobjCrossed->getRGBSelf (vCross) :	// reverse side
			    pobjCrossed->getRGBSelf (vCross));	// face side
		rgbNew += rgbSelf;
		
		// return here for 0-th Level Color
		if (nLevel <= 0)
			return	rgbSelf;
		
		// Reflection Color from (n-1)-th Level
		Vector	ReflDir	= ViewRay - 2.0*RN * NormalVec;
		RGB	rgbRefl	 = GetRGB (nLevel-1,
					   Line (vCross, vCross + ReflDir));
		rgbNew += rgbRefl  * pobjCrossed->getRGBRefl (vCross);
		
		// Transparent / Refraction Color from (n-1)-th Level
		Vector	vTrans = ViewRay;
		Real	rRefraction = pobjCrossed->rRefraction;	// Refraction Ratio
		if (bRefractionDir (ViewRay, NormalVec, rRefraction, vTrans)) {
			RGB	rgbTrans = GetRGB (nLevel-1,
						   Line (vCross, vCross + vTrans));
			rgbNew += rgbTrans * pobjCrossed->getRGBTrans (vCross);
		}
		
		// restrict Blur by Level
		if (nLevel <= 2)
			return	rgbNew;
		
		// for Reflection Bluring
		if (pobjCrossed->bReflBlur) {
			
			RGB	rgbBlurR;
			rgbBlurR = GetBlurRGB (vCross, ReflDir,
						pobjCrossed->nReflBlur,
						pobjCrossed->rReflBlurAngle);
			rgbNew += rgbBlurR * pobjCrossed->getRGBRefl (vCross);
		}
		
		// for Total Direction Bluring
		if (pobjCrossed->bTotalBlur) {
			
			Vector	NormalVecHere		// Normal Vec Dir. to here
				= (RN > 0 ?
				   NormalVec * (-1) :	// reverse side
				   NormalVec);		// face side
			RGB	rgbBlurT;
			rgbBlurT = GetBlurRGB (vCross, NormalVecHere,
						pobjCrossed->nTotalBlur);
			rgbNew += rgbBlurT * pobjCrossed->rgbTotalBlur;
		}
		
		return	rgbNew;
	}
	RGB GetBlurRGB (const Vector &vCrossPoint,	// eye point
			const Vector &vCenterDir,	// Bluring center dir
			int   nBluring,			// # of integral dir
			Real  rIntegralAngle = M_PI/2	// integral angle
		) const
	{
		Vector	vCenter (vCenterDir.unit ()),	// unit Center Dir
			vDir1 (vCenter .GetOrthogonalOne ()),	// unit ortho
			vDir2 (vCenter .out (vDir1));
		// (vCenter, vDir1, vDir2) Complete Orthogonal System
		
		RGB	rgbSum (0,0,0);
		for (int n = 0 ; n < nBluring ; n++) {
 			Real	s = rRand (),		// [0,0]-[1,1]
 				t = rRand ();		// random value
			
			t = 1.0 - t* (1.0 - cos (rIntegralAngle));
							// restrict integral
			
			Real	rPhi	= 2.0*M_PI * s,
				rRadius	= 1.0 * sqrt (1/(t*t) - 1.0);
			Vector	vViewDir = vCenter
					 + vDir1 * rRadius * cos (rPhi)
					 + vDir2 * rRadius * sin (rPhi);
			rgbSum += GetRGB (1, Line (vCrossPoint,
						   vCrossPoint + vViewDir));
		}
		rgbSum *= (1.0 / nBluring);
		return	rgbSum;
	}
	Bool bRefractionDir (const Vector &vViewDir,
			     const Vector &vNormalVec,
			     Real rRefractionRatio,
			     Vector &vRefractionDir) const
	{
		Vector	vD = vViewDir. unit (),		// unit Dir
			vN = vNormalVec. unit ();	// unit Normal
		Real	DN = vD .in (vN);		// inner product
		Real	rRR;
		
		if (DN < 0) {
			rRR = rRefractionRatio;		// outside to inside
		} else {
			rRR = 1.0 / rRefractionRatio;	// inside to outside
			vN *= -1;
			DN *= -1;
		}
		
		Real	rCheck = DN*DN + rRR - 1.0;
		if (rCheck < 0)
			return	False;			// over critical angle
		Real	t = DN + sqrt (rCheck);
		
		vRefractionDir = vD - t*vN;
		return	True;
	}
};


class PostScript {
private:
	Line	CenterViewLine;
	Vector	xDir,	yDir;
	
	int	nRecursiveLevel;	// Reflection Recursive Level
					// for most simple case it be 0.
	
	int	xSize,		ySize;
	Bool	bEps;			// output as EPS file
	Bool	bPPM;			// output as PPM raw file
	Bool	bShowStatus;		// show calculating status
	
public:
	PostScript (Line CenterViewLine1, Vector xD, Vector yD,
		    int xS = 100, int yS = 100)
	{
		CenterViewLine = CenterViewLine1;
		xDir	= xD;
		yDir	= yD;
		
		nRecursiveLevel = 0;
		
		xSize	= xS;
		ySize	= yS;
		bEps	= False;
		bPPM	= False;
		bShowStatus = False;
	}
	int nCutIn (int n)
	{
		if (n < 0)
			return	0;
		if (n > 255)
			return	255;
		return	n;
	}
	RGB GetRGB (const List &LS, int x, int y)
	{
		Vector	EyePoint (CenterViewLine.v),
			CenterPoint (CenterViewLine.w);
		
		Real	rX = 2.0 * x / xSize - 1.0,
			rY = 2.0 * y / ySize - 1.0;
		
		Line	ViewLine (EyePoint, CenterPoint
					    + xDir*rX + yDir*rY);
		
		return	LS.GetRGB (nRecursiveLevel, ViewLine);
	}
	PostScript & putImagFile (const List &LS,
				  String szFileNameHead) // const
	{
		char	szFileName [128],	*szExt;
		
		if (bPPM) {
			szExt = ".ppm";
		} else if (bEps) {
			szExt = ".eps";
		} else {
			szExt = ".ps";
		}
		sprintf (szFileName, "%s%s", szFileNameHead, szExt);
		ofstream	fout (szFileName);
		
		if (bShowStatus)
			cerr << szFileName << endl;
		
		if (bPPM) {
			putAsPPMRaw (LS, fout);
		} else {
			put (LS, fout);
		}
	}
	PostScript & putPS (const List &LS, ostream &s = cout) // const
	{
		putPrologue (s);
		put (LS,s);
		putEpilogue (s);
	}
	PostScript & put (const List &LS, ostream &s = cout) // const
	{
		s << hex;
		
		for (int y = ySize -1 ; y >= 0 ; y--) {
			for (int x = 0 ; x < xSize ; x++) {
				
				RGB rgb = GetRGB (LS, x, y);
				//rgb.cutin ();
				
				int	nR = nCutIn ((int)(256.0*rgb.rRed)),
					nG = nCutIn ((int)(256.0*rgb.rGreen)),
					nB = nCutIn ((int)(256.0*rgb.rBlue));
				
 				s << setw (2) << setfill ('0') << nR
 				  << setw (2) << setfill ('0') << nG
 				  << setw (2) << setfill ('0') << nB;
				
				if (x % 8 == 7)
					s << '\n';
			}
			if ((xSize-1) % 8 != 7)
				s << '\n';
			
			if (bShowStatus && y % 10 == 0)
				cerr << setfill (' ')  << setw (4) << y;
			if (bShowStatus && y % 100 == 0) 
				cerr << endl;
		}
		if (bShowStatus) {
			// cerr << endl;
		}
		
		s << dec << endl;
		return	*this;
	}
	PostScript & putAsPPMRaw (const List &LS, ostream &s = cout) // const
	{
		s << "P6" << endl;
		s << "# Image Generated by Y.Nagatani RayTracer" << endl;
		s << xSize << " " << ySize << endl;
		s << 255 << endl;
		
		for (int y = ySize -1 ; y >= 0 ; y--) {
			for (int x = 0 ; x < xSize ; x++) {
				
				RGB rgb = GetRGB (LS, x, y);
				
				int	nR = nCutIn ((int)(256.0*rgb.rRed)),
					nG = nCutIn ((int)(256.0*rgb.rGreen)),
					nB = nCutIn ((int)(256.0*rgb.rBlue));
				
				unsigned char	szBuf [4];
				
				s << ((unsigned char) nR)
				  << ((unsigned char) nG)
				  << ((unsigned char) nB);
			}
			if (bShowStatus && y % 10 == 0)
				cerr << setfill (' ')  << setw (4) << y;
			if (bShowStatus && y % 100 == 0) 
				cerr << endl;
		}
		if (bShowStatus) {
			// cerr << endl;
		}
		// s << endl;
		return	*this;
	}
	PostScript & putPrologue (ostream &s = cout)
	{
		Real	xPaperSize = 210.0 * 72.0 / 25.4,	// (point)
			yPaperSize = 297.0 * 72.0 / 25.4,
			xOutSize   = 500,			// (point)
			yOutSize   = 500,
			xSt  = (xPaperSize - xOutSize)/2,	// (point)
			ySt  = (yPaperSize - yOutSize)/2,
			xEnd = xSt + xOutSize,			// (point)
			yEnd = ySt + yOutSize;
		
		if (!bEps) {	
			s <<	"%!PS-Adobe-2.0\n";
		} else {
			s <<	"%!PS-Adobe-2.0 EPSF-2.0\n";
		}
		s <<	"%%Creator: raytrace (C) 1992,1997 Y.Nagatani\n"
			"%%Title: Ray Tracing\n";
		
		if (!bEps) {
			s <<	"%%Pages: 1\n"
				"%%DocumentPaperSizes: A4\n"
				"%%Orientation: Portrait\n";
		} else {
			
			s <<	"%%BoundingBox: "
			  <<	(int) xSt  << ' ' << (int) ySt  << ' '
			  <<	(int) xEnd << ' ' << (int) yEnd << '\n';
		}
		s <<	"%%EndComments\n"
			"\n"
			"/mm {72 25.4 div mul} def\n"
			"/xPaperSize 210 mm def		% A3 size\n"
			"/yPaperSize 297 mm def\n"
			"\n"
			"%%EndProlog\n"
			"%%Page: 1\n"
			"\n"
			"/DefaultState save def\n"
			"20 dict begin\n"
			"\n"
			"/xNMax " << xSize    << " def	% Matrix Size\n"
			"/yNMax " << ySize    << " def\n"
			"\n"
			"/xSize " << xOutSize << " def	% Region Size (pt)\n"
			"/ySize " << yOutSize << " def\n"
			"\n"
			"/Init {			% initialize func.\n"
			"	xPaperSize xSize neg add 2 div\n"
			"	yPaperSize ySize neg add 2 div translate\n"
			"	\n"
			"	xSize xNMax div \n"
			"	ySize yNMax div scale\n"
			"	1 setlinewidth\n"
			"} bind def\n"
			"\n"
			"/Buf xNMax 3 mul string def	% line buffer\n"
			"\n"
			"%% Main\n"
			"gsave\n"
			"Init\n"
			"\n"
			"xNMax yNMax 8			% image size & depth\n"
			"[1 0 0 -1 0 yNMax]		% affine matrix\n"
			"{currentfile Buf readhexstring pop}\n"
			"false 3 colorimage\n"
			"\n";
		return	*this;
	}
	PostScript & putEpilogue (ostream &s = cout)
	{
		s <<	"\n"
			"grestore\n"
			"showpage\n"
			"\n"
			"end\n"
			"DefaultState restore\n"
			"%%Trailer\n"
			"%%EOF\n";
		return	*this;
	}
	PostScript & setEps (Bool b = True)
	{
		bEps = b;
		return	*this;
	}
	PostScript & setPPM (Bool b = True)
	{
		bPPM = b;
		return	*this;
	}
	PostScript & setShowStatus (Bool b = True)
	{
		bShowStatus = b;
		return	*this;
	}
	PostScript & setRecursiveLevel (int n = 0)
	{
		nRecursiveLevel = n;
		return	*this;
	}
};

////// Patterned Triangle and Patterned Quadrilateral //////

class TrianglePatt : public Triangle {
public:
	TrianglePatt (void)				// default constructor
		: Triangle ()
	{
	}
	TrianglePatt (const Vector v1,
		  const Vector v2,
		  const Vector v3)			// constructor
		: Triangle (v1, v2, v3)
	{
	}
	TrianglePatt (Real x1, Real y1, Real z1,
		  Real x2, Real y2, Real z2,
		  Real x3, Real y3, Real z3)		// constructor
		: Triangle (x1,y1,z1, x2,y2,z2, x3,y3,z3)
	{
	}
	TrianglePatt (const Triangle &Q)	// copy constructor
		: Triangle (Q)
	{
	}
	virtual Object * copyNew (void) const
	{
		return	new TrianglePatt (*this);
	}
	virtual RGB getRGBSelf (Vector vPoint)
	{
		// vPoint must be on Triangle Plane
		Real	x,	y;
		getBaseCoefs (vPoint, &x, &y);
		
		if (x < 0 || y < 0 ||
		    x > 1 || y > 1)
			return	RGB (0,0,0.0);
			// Line Cross Plane
			// but not Cross Area
		
		Real	A = 5.0*x,
			B = 5.0*y;
		Real	r = (A - floor (A)),
			g = (B - floor (B)),
			b = 1.0;
		
		if (r + g > 1)
			r = g = b = 0.0;
		
		return	RGB (r, g, b);
	}
	virtual RGB getRGBTrans (Vector vPoint)
	{
		// vPoint must be on Triangle Plane
		Real	x,	y;
		getBaseCoefs (vPoint, &x, &y);
		
		if (x < 0 || y < 0 ||
		    x > 1 || y > 1)
			return	RGB (0,0,0.0);
			// Line Cross Plane but not Cross Area
		
		Real	A = 5.0*x,
			B = 5.0*y;
		Real	r = (A - floor (A)),
			g = (B - floor (B)),
			b = 1.0;
		int	nA = ((int) A) % 2,
			nB = ((int) B) % 2;
		
		if (r + g > 1) {
			r = g = b = 1.0;
			if (!nA && !nB)
				r = 0.0;
			if ( nA && !nB)
				g = 0.0;
			if (!nA &&  nB)
				b = 1.0;
			if ( nA &&  nB)
				r = g = b = 0.2;
		} else {
			r = g = b = 0.0;
		}
		
		return	RGB (r, g, b);
	}
};


class QuadrilateralPatt : public Quadrilateral {
public:
	QuadrilateralPatt (void)			// default constructor
		: Quadrilateral ()
	{
	}
	QuadrilateralPatt (const Vector v1,
		  const Vector v2,
		  const Vector v3)			// constructor
		: Quadrilateral (v1, v2, v3)
	{
	}
	QuadrilateralPatt (Real x1, Real y1, Real z1,
		  Real x2, Real y2, Real z2,
		  Real x3, Real y3, Real z3)		// constructor
		: Quadrilateral (x1,y1,z1, x2,y2,z2, x3,y3,z3)
	{
	}
	QuadrilateralPatt (const Quadrilateral &Q)	// copy constructor
		: Quadrilateral (Q)
	{
	}
	virtual Object * copyNew (void) const
	{
		return	new QuadrilateralPatt (*this);
	}
	virtual RGB getRGBSelf (Vector vPoint)
	{
		// vPoint must be on Triangle Plane
		Real	x,	y;
		getBaseCoefs (vPoint, &x, &y);
		
		if (x < 0 || y < 0 ||
		    x > 1 || y > 1)
			return	RGB (0,0,0.0);
			// Line Cross Plane but not Cross Area
		
		Real	A = 15.0*x,
			B = 15.0*y;
		Real	r = (A - floor (A)),
			g = (B - floor (B));
		Real	b = 0.0;
		int	a1 = (((int)A) % 2),
			b1 = (((int)B) % 2);
		if (a1 && b1) {
			r = g = b = 0.0;
		}
		if (!a1 && !b1) {
			b = 1.0;
		}
		return	RGB (r, g, b);
	}
};

////// PPM Raw Ascii Image Loader //////  QuadrilateralImag

class PPMRawImage {
public:
	int		xSize,	ySize;
	unsigned char	*im;
	PPMRawImage (void) {
		xSize = 0;
		ySize = 0;
		im = NULL;
	}
	PPMRawImage (char *szFileName) {
		char	szBuffer [256];
		FILE	*fp = fopen (szFileName, "r");
		Bool	bAscii = True;
		if (fp == NULL) {
			cerr << "Cannot Open File : " << szFileName << endl;
			exit (1);
		}
		char	* ret = fgets (szBuffer, 255, fp);
		if (ret == NULL ||
		     szBuffer[0] != 'P' ||
		    (szBuffer[1] != '3' &&
		     szBuffer[1] != '6')) {
			cerr << "Non PPM (raw/ascii) File : "
			     << szFileName << endl;
			exit (1);
		}
		if (szBuffer[1] == '3')
			bAscii = True;
		else
			bAscii = False;
		
		cerr << "Reading Image File <" << szFileName
		     << "> as PPM(" << (bAscii ? "Ascii" : "Raw")
		     <<") File" << endl;
		
		for (;;) {
			if (0 == fgets (szBuffer, 255, fp)) {
				cerr << "Non PPM File (Line Shortage): "
				     << szFileName << endl;
				exit (1);
			}
			if (szBuffer[0] == '#')
				continue;
			if (2 == sscanf (szBuffer, "%d%d", &xSize, &ySize)) {
				cerr << "Image Size is "
				     << xSize << "x" << ySize << endl;
				break;
			} else {
				cerr << "Cannot get Image Size." << endl;
				exit (1);
			}
		}
		if (NULL == fgets (szBuffer, 255, fp)) {
			cerr << "Non PPM File (Line Shortage): "
			     << szFileName << endl;
			exit (1);
		}
		
		im = new unsigned char [xSize * ySize * 3 + 15];
		if (im == NULL) {
			cerr << "Cannot Allocate Memory for Image."
			     << endl;
			exit (0);
		}
		if (bAscii)
			ReadAsPPMAscii (fp);
		else
			ReadAsPPMRaw (fp);
	}
	void ReadAsPPMAscii (FILE *fp) {
		int	pos = 0;
		for (;;) {
			char	szBuffer [256];
			if (NULL == fgets (szBuffer, 255, fp)) {
				cerr << "Non PPM File (Line Shortage)."
				     << endl;
				exit (1);
			}
			int	nGet;
			int	f [5 * 3];
			nGet = sscanf (szBuffer,
				       "%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d",
				       f+ 0, f+ 1, f+ 2, f+3, f+ 4, f+ 5,
				       f+ 6, f+ 7, f+ 8, f+9, f+10, f+11,
				       f+12, f+13, f+14);
			for (int n = 0 ; n < nGet ; n++) {
				im [pos++] = (unsigned char) f [n];
			}
			if (pos >= xSize * ySize * 3)
				break;
		}		
	}	
	void ReadAsPPMRaw (FILE *fp) {
		int	pos = 0;
		for (;;) {
			int	nCh = fgetc (fp);
			if (0 == EOF) {
				cerr << "Non PPM File (Line Shortage)."
				     << endl;
				exit (1);
			}
			
			im [pos++] = (unsigned char) nCh;
			
			if (pos >= xSize * ySize * 3)
				break;
		}		
	}	
	Real getValue (int x, int y, int c) const {
		if (0 <= x && x < xSize &&
		    0 <= y && y < ySize) {
			y = (ySize - 1) - y;
			return	(Real) im[3*(x + xSize*y) + c] / 256.0;
		} else {
			return	0.0;
		}
	}
	Real getValue (Real rX, Real rY, int c) const {
		return	getValue ((int)((Real) xSize * rX),
				  (int)((Real) ySize * rY),
				  c);
	}
	void show (void) {
		for (int y = 0 ; y < ySize ; y++) {
			for (int x = 0 ; x < xSize ; x++) {
				cerr << "[" << (int) im[3*(x + xSize*y) + 0]
				     << "," << (int) im[3*(x + xSize*y) + 1]
				     << "," << (int) im[3*(x + xSize*y) + 2]
				     << "] ";
			}
		}
	}
};


class QuadrilateralImag : public Quadrilateral {
public:
	PPMRawImage const	*ppm;
	Bool	bNoTranspareceWithoutBlack;
	
	QuadrilateralImag (void)			// default constructor
		: Quadrilateral ()
	{
		bNoTranspareceWithoutBlack = False;
	}
	QuadrilateralImag (const Vector v1,
			   const Vector v2,
			   const Vector v3,
			   const PPMRawImage & ppm)	// constructor
		: Quadrilateral (v1, v2, v3)
	{
		this->ppm = & ppm;
		bNoTranspareceWithoutBlack = False;
	}
	QuadrilateralImag (Real x1, Real y1, Real z1,
			   Real x2, Real y2, Real z2,
			   Real x3, Real y3, Real z3,
			   const PPMRawImage & ppm)	// constructor
		: Quadrilateral (x1,y1,z1, x2,y2,z2, x3,y3,z3)
	{
		this->ppm = & ppm;
		bNoTranspareceWithoutBlack = False;
	}
	QuadrilateralImag (const Quadrilateral &Q,
			   const PPMRawImage & ppm)	// copy constructor
		: Quadrilateral (Q)
	{
		this->ppm = & ppm;
		bNoTranspareceWithoutBlack = False;
	}
	virtual Object * copyNew (void) const
	{
		return	new QuadrilateralImag (*this);
	}
	virtual RGB getRGBSelf (Vector vPoint)
	{
		// vPoint must be on Triangle Plane
		Real	x,	y;
		getBaseCoefs (vPoint, &x, &y);
		
		if (x < 0 || y < 0 ||
		    x > 1 || y > 1)
			return	RGB (0,0,0.0);
			// Line Cross Plane but not Cross Area
		
		Real	r = ppm->getValue (x,y, 0),
			g = ppm->getValue (x,y, 1),
			b = ppm->getValue (x,y, 2);
		return	RGB (r, g, b);
	}
	virtual RGB getRGBTrans (Vector vPoint)
	{
		if (!bNoTranspareceWithoutBlack)
			return	Quadrilateral::getRGBTrans (vPoint);
		
		// vPoint must be on Triangle Plane
		Real	x,	y;
		getBaseCoefs (vPoint, &x, &y);
		
		if (x < 0 || y < 0 ||
		    x > 1 || y > 1)
			return	RGB (1,1,1);
			// Line Cross Plane but not Cross Area
			// then perfect trancepearent
		
		Real	r = ppm->getValue (x,y, 0),
			g = ppm->getValue (x,y, 1),
			b = ppm->getValue (x,y, 2);
		Real	rSmall = 0.001;
		
		if (r > rSmall || g > rSmall || b > rSmall)
			return	RGB (0,0,0);	// No Trancepearent
		else
			return	RGB (1,1,1);	// Perfect Trancepearent
						// for Black Color
	}
	QuadrilateralImag & setNoTranspareceWithoutBlack (Bool bSet = True)
	{
		bNoTranspareceWithoutBlack = bSet;
		return	*this;
	}
};