Index: GTP/trunk/App/Demos/Vis/KdTreeDemo/OGRE/scripts/TestKdTree.vcproj
===================================================================
--- GTP/trunk/App/Demos/Vis/KdTreeDemo/OGRE/scripts/TestKdTree.vcproj	(revision 2361)
+++ GTP/trunk/App/Demos/Vis/KdTreeDemo/OGRE/scripts/TestKdTree.vcproj	(revision 2746)
@@ -86,5 +86,5 @@
 				OptimizeForWindowsApplication="TRUE"
 				AdditionalIncludeDirectories="..\include;&quot;$(OGRE_PATH)\Samples\Common\include&quot;;&quot;$(OGRE_PATH)\OgreMain\include&quot;;..\..\..\..\..\..\Lib\Vis\OnlineCullingCHC\IVReader\include;..\..\..\..\..\..\Lib\Vis\OnlineCullingCHC\OGRE\include;..\..\..\..\..\..\Lib\Vis\OnlineCullingCHC\include"
-				PreprocessorDefinitions="WIN32;NDEBUG;_WINDOWS;GTP_VISIBILITY_MODIFIED_OGRE;GAMETOOLS_ILLUMINATION_MODULE"
+				PreprocessorDefinitions="WIN32;NDEBUG;_WINDOWS;GTP_VISIBILITY_MODIFIED_OGRE;GAMETOOLS_ILLUMINATION_MODULE1"
 				StringPooling="TRUE"
 				MinimalRebuild="TRUE"
@@ -263,5 +263,5 @@
 				Name="VCAuxiliaryManagedWrapperGeneratorTool"/>
 		</Configuration>
-<Configuration
+		<Configuration
 			Name="Release_noqt|Win32"
 			OutputDirectory="Release"
Index: GTP/trunk/App/Demos/Vis/HillyTerrain/OGRE/TestCullingTerrain.vcproj
===================================================================
--- GTP/trunk/App/Demos/Vis/HillyTerrain/OGRE/TestCullingTerrain.vcproj	(revision 2501)
+++ GTP/trunk/App/Demos/Vis/HillyTerrain/OGRE/TestCullingTerrain.vcproj	(revision 2746)
@@ -83,5 +83,5 @@
 				OptimizeForWindowsApplication="TRUE"
 				AdditionalIncludeDirectories="&quot;$(OGRE_PATH)\Dependencies\include&quot;;&quot;$(OGRE_PATH)\OgreMain\include&quot;;&quot;$(OGRE_PATH)\Samples\Common\include&quot;;&quot;$(OGRE_PATH)\Dependencies\include\CEGUI&quot;;&quot;$(OGRE_PATH)\PlugIns\OctreeSceneManager\include&quot;;&quot;$(OGRE_PATH)\Samples\Common\CEGUIRenderer\include&quot;;..\..\Ogre\include;..\..\..\..\..\Lib\Vis\OnlineCullingCHC\include;..\include;..\..\..\..\..\Lib\Vis\OnlineCullingCHC\OGRE\include;..\..\..\..\..\Lib\Vis\OnlineCullingCHC\IVReader\include"
-				PreprocessorDefinitions="_WINDOWS;_STLP_USE_DYNAMIC_LIB;OGRE_LIBRARY_IMPORTS;_RELEASE;WIN32;_STLP_RELEASE;GTP_VISIBILITY_MODIFIED_OGRE;GAMETOOLS_ILLUMINATION_MODULE"
+				PreprocessorDefinitions="_WINDOWS;_STLP_USE_DYNAMIC_LIB;OGRE_LIBRARY_IMPORTS;_RELEASE;WIN32;_STLP_RELEASE;GTP_VISIBILITY_MODIFIED_OGRE;GAMETOOLS_ILLUMINATION_MODULE1"
 				RuntimeLibrary="2"
 				RuntimeTypeInfo="FALSE"
@@ -142,5 +142,5 @@
 				OptimizeForWindowsApplication="TRUE"
 				AdditionalIncludeDirectories="&quot;$(OGRE_PATH)\Dependencies\include&quot;;&quot;$(OGRE_PATH)\OgreMain\include&quot;;&quot;$(OGRE_PATH)\Samples\Common\include&quot;;&quot;$(OGRE_PATH)\Dependencies\include\CEGUI&quot;;&quot;$(OGRE_PATH)\PlugIns\OctreeSceneManager\include&quot;;&quot;$(OGRE_PATH)\Samples\Common\CEGUIRenderer\include&quot;;..\..\Ogre\include;..\..\..\..\..\Lib\Vis\OnlineCullingCHC\include;..\include;..\..\..\..\..\Lib\Vis\OnlineCullingCHC\OGRE\include;..\..\..\..\..\Lib\Vis\OnlineCullingCHC\IVReader\include;&quot;$(LIBXML_HOME)\include&quot;;d:\work\Collada\include;..\..\..\..\..\Lib\Vis\Preprocessing\src\Timer"
-				PreprocessorDefinitions="_WINDOWS;_STLP_USE_DYNAMIC_LIB;OGRE_LIBRARY_IMPORTS;_RELEASE;WIN32;_STLP_RELEASE;GTP_VISIBILITY_MODIFIED_OGRE;GAMETOOLS_ILLUMINATION_MODULE"
+				PreprocessorDefinitions="_WINDOWS;_STLP_USE_DYNAMIC_LIB;OGRE_LIBRARY_IMPORTS;_RELEASE;WIN32;_STLP_RELEASE;GTP_VISIBILITY_MODIFIED_OGRE;GAMETOOLS_ILLUMINATION_MODUL"
 				StringPooling="TRUE"
 				RuntimeLibrary="2"
Index: GTP/trunk/App/Demos/Vis/HillyTerrain/OGRE/TestCullingTerrainApplication.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/HillyTerrain/OGRE/TestCullingTerrainApplication.cpp	(revision 2555)
+++ GTP/trunk/App/Demos/Vis/HillyTerrain/OGRE/TestCullingTerrainApplication.cpp	(revision 2746)
@@ -234,4 +234,6 @@
 	mVizCamera = mSceneMgr->createCamera("VizCam");
 	mVizCamera->setPosition(mCamNode->getPosition());
+
+	///////////////
 
 	mVizCamera->setNearClipDistance(1);
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/DataTypes.c
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/DataTypes.c	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/DataTypes.c	(revision 2746)
@@ -18,19 +18,27 @@
 }
 
-void addVector3(Vector3 result, const Vector3 a, const Vector3 b) {
+void addVector3(Vector3 result, const Vector3 a, const Vector3 b) 
+{
 	int i;
-	for(i = 0; i < 3; i++) {
-		result[i] = a[i]+b[i];
-	}
-}
-
-void diffVector3(Vector3 result, const Vector3 a, const Vector3 b) {
+	
+	for (i = 0; i < 3; ++ i) 
+	{
+		result[i] = a[i] + b[i];
+	}
+}
+
+
+void diffVector3(Vector3 result, const Vector3 a, const Vector3 b) 
+{
 	int i;
-	for(i = 0; i < 3; i++) {
-		result[i] = a[i]-b[i];
-	}
-}
-
-void copyVector3Values(Vector3 result, const double x, const double y, const double z) {
+	for (i = 0; i < 3; ++ i) 
+	{
+		result[i] = a[i] - b[i];
+	}
+}
+
+
+void copyVector3Values(Vector3 result, const double x, const double y, const double z) 
+{
 	result[0] = x; 
 	result[1] = y;
@@ -38,15 +46,23 @@
 }
 
-void copyMatrix(Matrix4x4 result, const Matrix4x4 input) {
-	if(input != result) {
+
+void copyMatrix(Matrix4x4 result, const Matrix4x4 input) 
+{
+	if (input != result) 
+	{
 		memcpy(result,input,4*4*sizeof(double));
 	}
 }
 
-void planesSetSize(struct VecPlane* p, const int size) {
-	if(0 != p) {
-		if(size == p->size) {
+
+void planesSetSize(struct VecPlane* p, const int size) 
+{
+	if (0 != p) 
+	{
+		if(size == p->size) 
+		{
 			return;
 		}
+
 		p->plane = (struct Plane*)realloc(p->plane,size*sizeof(struct Plane));
 		p->size = size;
@@ -54,13 +70,20 @@
 }
 
-void emptyVecPlane(struct VecPlane* p) {
+
+void emptyVecPlane(struct VecPlane* p) 
+{
 	planesSetSize(p,0);
 }
 
-void vecPointSetSize(struct VecPoint* poly, const int size) {
-	if(0 != poly) {
-		if(size == poly->size) {
+
+void vecPointSetSize(struct VecPoint* poly, const int size) 
+{
+	if (0 != poly) 
+	{
+		if (size == poly->size) 
+		{
 			return;
 		}
+		
 		poly->points = (Vector3*)realloc(poly->points,size*sizeof(Vector3));
 		poly->size = size;
@@ -68,70 +91,96 @@
 }
 
-void emptyVecPoint(struct VecPoint* poly) {
+
+void emptyVecPoint(struct VecPoint* poly) 
+{
 	vecPointSetSize(poly,0);
 }
 
-void append2VecPoint(struct VecPoint* poly, const Vector3 p) {
-	if(0 != poly) {
+
+void append2VecPoint(struct VecPoint* poly, const Vector3 p) 
+{
+	if (0 != poly) 
+	{
 		int size = poly->size;
 		vecPointSetSize(poly,size+1);
-		copyVector3(poly->points[size],p);
-	}
-}
-
-void copyVecPoint(struct VecPoint* poly, const struct VecPoint poly2) {
-	if(0 != poly) {
+		copyVector3(poly->points[size], p);
+	}
+}
+
+
+void copyVecPoint(struct VecPoint* poly, const struct VecPoint poly2) 
+{
+	if (0 != poly) 
+	{
 		int i;
 		vecPointSetSize(poly,poly2.size);
-		for(i= 0; i < poly2.size; i++) {
-			copyVector3(poly->points[i],poly2.points[i]);
-		}
-	}
-}
-
-void swapVecPoint(struct VecPoint* poly1, struct VecPoint* poly2) {
-	if(0 != poly1 && 0 != poly2 && poly1 != poly2) {
-		{
-			Vector3* points = poly1->points;
-			poly1->points = poly2->points;
-			poly2->points = points;
-		}
-		{
-			int size = poly1->size;
-			poly1->size = poly2->size;
-			poly2->size = size;
-		}
-	}
-}
-
-void objectSetSize(struct Object* obj, const int size) {
-	if(0 != obj) {
+		for (i = 0; i < poly2.size; ++ i) 
+		{
+			copyVector3(poly->points[i], poly2.points[i]);
+		}
+	}
+}
+
+
+void swapVecPoint(struct VecPoint* poly1, struct VecPoint* poly2) 
+{
+	if(0 != poly1 && 0 != poly2 && poly1 != poly2) 
+	{
+		Vector3* points = poly1->points;
+		poly1->points = poly2->points;
+		poly2->points = points;
+		
+		int size = poly1->size;
+		poly1->size = poly2->size;
+		poly2->size = size;
+	}
+}
+
+
+void objectSetSize(struct Object* obj, const int size) 
+{
+	if (0 != obj) 
+	{
 		int i;
-		if(size == obj->size) {
+		if(size == obj->size) 
+		{
 			return;
 		}
+		
 		//dispose if shrinking
-		for(i = size; i < obj->size; i++) {
-			emptyVecPoint( &(obj->poly[i]) ); 
-		}
+		for(i = size; i < obj->size; i++) 
+		{
+			emptyVecPoint(&(obj->poly[i])); 
+		}
+
 		//allocate new place
 		obj->poly = (struct VecPoint*)realloc(obj->poly,size*sizeof(struct VecPoint));
+		
 		//initialize new place
-		for(i = obj->size; i < size; i++) {
+		for(i = obj->size; i < size; i++) 
+		{
 			obj->poly[i] = VECPOINT_NULL;
 		}
+		
 		obj->size = size;
 	}
 }
 
-void emptyObject(struct Object* obj) {
-	objectSetSize(obj,0);
-}
-
-void copyObject(struct Object* obj, const struct Object objIn) {
-	if(0 != obj) {
+
+void emptyObject(struct Object* obj) 
+{
+	objectSetSize(obj, 0);
+}
+
+
+void copyObject(struct Object* obj, const struct Object objIn) 
+{
+	if(0 != obj) 
+	{
 		int i;
 		objectSetSize(obj,objIn.size);
-		for(i = 0; i < objIn.size; i++) {
+
+		for(i = 0; i < objIn.size; i++)
+		{
 			copyVecPoint( &(obj->poly[i]), objIn.poly[i]);
 		}
@@ -139,14 +188,18 @@
 }
 
-void append2Object(struct Object* obj, const struct VecPoint poly) {
-	if(0 != obj) {
+
+void append2Object(struct Object* obj, const struct VecPoint poly) 
+{	
+	if (0 != obj) 
+	{
 		int size = obj->size;
-		objectSetSize(obj,size+1);
-		copyVecPoint( &(obj->poly[size]) ,poly);
-	}
-}
-
-
-void convObject2VecPoint(struct VecPoint* points,const struct Object obj) {
+		objectSetSize(obj, size + 1);
+		copyVecPoint(&(obj->poly[size]), poly);
+	}
+}
+
+
+void convObject2VecPoint(struct VecPoint* points,const struct Object obj) 
+{
 	if(0 != points) {
 		int i, j;
@@ -160,4 +213,5 @@
 	}
 }
+
 
 void calcAABoxPoints(Vector3x8 points, const struct AABox b) {
@@ -173,4 +227,5 @@
 }
 
+
 void calcAABoxPlanes(struct VecPlane* planes, const struct AABox b) {
 	if(0 != planes) {
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Geometry.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Geometry.h	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Geometry.h	(revision 2746)
@@ -16,7 +16,6 @@
 	Geometry();
 	Geometry(Vector3 translation, float xRot, float yRot, float zRot, float scale, int objectType);
-	//! renders this geometry
+	//! Renders this geometry
 	void Render();
-	
 	//! sets rotations in degree : executed in the order x, y, and z rotation
 	void SetRotations(float xRot, float yRot, float zRot);
@@ -33,20 +32,16 @@
 	void GetRotations(float &xRot, float &yRot, float &zRot);
 
-	// --- material settings
-	void SetAmbientColor(float ambientR,   float ambientG,  float ambientB);
-	void SetDiffuseColor(float diffuseR,   float diffuseG,  float diffuseB);
-	void SetSpecularColor(float specularR, float specularG, float specularB);
-
+	//! Set pointer to the material
+	void SetMaterial(Material *mat);
+	//! Returns the bounding box
 	const AABox& GetBoundingVolume();
-		
 	//! last visited flag, important for rendering each geometry only once
 	void SetLastVisited(int lastVisited);
 	//! returns lastvisited flag
 	int GetLastVisited();
-
 	// --- sets the object type to one of teapot, torus, sphere
 	void SetObjectType(int type);
 
-	enum{TEAPOT, TORUS, SPHERE, NUM_OBJECTS};
+	
 	
 	static int sDisplayList[NUM_OBJECTS];
@@ -60,4 +55,5 @@
 
 protected:
+
 	//! generates the display list this a object
 	void GenerateList();
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Matrix4x4.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Matrix4x4.cpp	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Matrix4x4.cpp	(revision 2746)
@@ -0,0 +1,637 @@
+#include "Matrix4x4.h"
+#include "Vector3.h"
+
+// standard headers
+#include <iomanip>
+using namespace std;
+
+//namespace GtpVisibilityPreprocessor {
+
+
+Matrix4x4::Matrix4x4(const Vector3 &a, const Vector3 &b, const Vector3 &c)
+{
+	// first index is column [x], the second is row [y]
+	x[0][0] = a.x;
+	x[1][0] = b.x;
+	x[2][0] = c.x;
+	x[3][0] = 0.0f;
+
+	x[0][1] = a.y;
+	x[1][1] = b.y;
+	x[2][1] = c.y;
+	x[3][1] = 0.0f;
+
+	x[0][2] = a.z;
+	x[1][2] = b.z;
+	x[2][2] = c.z;
+	x[3][2] = 0.0f;
+
+	x[0][3] = 0.0f;
+	x[1][3] = 0.0f;
+	x[2][3] = 0.0f;
+	x[3][3] = 1.0f;
+}
+
+
+void Matrix4x4::SetColumns(const Vector3 &a, const Vector3 &b, const Vector3 &c)
+{
+	// first index is column [x], the second is row [y]
+	x[0][0] = a.x;
+	x[1][0] = a.y;
+	x[2][0] = a.z;
+	x[3][0] = 0.0f;
+
+	x[0][1] = b.x;
+	x[1][1] = b.y;
+	x[2][1] = b.z;
+	x[3][1] = 0.0f;
+
+	x[0][2] = c.x;
+	x[1][2] = c.y;
+	x[2][2] = c.z;
+	x[3][2] = 0.0f;
+
+	x[0][3] = 0.0f;
+	x[1][3] = 0.0f;
+	x[2][3] = 0.0f;
+	x[3][3] = 1.0f;
+}
+
+// full constructor
+Matrix4x4::Matrix4x4(float x11, float x12, float x13, float x14,
+					 float x21, float x22, float x23, float x24,
+					 float x31, float x32, float x33, float x34,
+					 float x41, float x42, float x43, float x44)
+{
+	// first index is column [x], the second is row [y]
+	x[0][0] = x11;
+	x[1][0] = x12;
+	x[2][0] = x13;
+	x[3][0] = x14;
+
+	x[0][1] = x21;
+	x[1][1] = x22;
+	x[2][1] = x23;
+	x[3][1] = x24;
+
+	x[0][2] = x31;
+	x[1][2] = x32;
+	x[2][2] = x33;
+	x[3][2] = x34;
+
+	x[0][3] = x41;
+	x[1][3] = x42;
+	x[2][3] = x43;
+	x[3][3] = x44;
+}
+
+
+// inverse matrix computation gauss_jacobiho method .. from N.R. in C
+// if matrix is regular = computatation successfull = returns 0
+// in case of singular matrix returns 1
+int Matrix4x4::Invert()
+{
+	int indxc[4],indxr[4],ipiv[4];
+	int i,icol,irow,j,k,l,ll,n;
+	float big,dum,pivinv,temp;
+	// satisfy the compiler
+	icol = irow = 0;
+
+	// the size of the matrix
+	n = 4;
+
+	for ( j = 0 ; j < n ; j++) /* zero pivots */
+		ipiv[j] = 0;
+
+	for ( i = 0; i < n; i++)
+	{
+		big = 0.0;
+		for (j = 0 ; j < n ; j++)
+			if (ipiv[j] != 1)
+				for ( k = 0 ; k<n ; k++)
+				{
+					if (ipiv[k] == 0)
+					{
+						if (fabs(x[k][j]) >= big)
+						{
+							big = fabs(x[k][j]);
+							irow = j;
+							icol = k;
+						}
+					}
+					else
+						if (ipiv[k] > 1)
+							return 1; /* singular matrix */
+				}
+				++(ipiv[icol]);
+				if (irow != icol)
+				{
+					for ( l = 0 ; l<n ; l++)
+					{
+						temp = x[l][icol];
+						x[l][icol] = x[l][irow];
+						x[l][irow] = temp;
+					}
+				}
+				indxr[i] = irow;
+				indxc[i] = icol;
+				if (x[icol][icol] == 0.0)
+					return 1; /* singular matrix */
+
+				pivinv = 1.0f / x[icol][icol];
+				x[icol][icol] = 1.0f ;
+				for ( l = 0 ; l<n ; l++)
+					x[l][icol] = x[l][icol] * pivinv ;
+
+				for (ll = 0 ; ll < n ; ll++)
+					if (ll != icol)
+					{
+						dum = x[icol][ll];
+						x[icol][ll] = 0.0;
+						for ( l = 0 ; l<n ; l++)
+							x[l][ll] = x[l][ll] - x[l][icol] * dum ;
+					}
+	}
+	for ( l = n; l--; )
+	{
+		if (indxr[l] != indxc[l])
+			for ( k = 0; k<n ; k++)
+			{
+				temp = x[indxr[l]][k];
+				x[indxr[l]][k] = x[indxc[l]][k];
+				x[indxc[l]][k] = temp;
+			}
+	}
+
+	return 0 ; // matrix is regular .. inversion has been succesfull
+}
+
+// Invert the given matrix using the above inversion routine.
+Matrix4x4
+Invert(const Matrix4x4& M)
+{
+	Matrix4x4 InvertMe = M;
+	InvertMe.Invert();
+	return InvertMe;
+}
+
+
+// Transpose the matrix.
+void
+Matrix4x4::Transpose()
+{
+	for (int i = 0; i < 4; i++)
+		for (int j = i; j < 4; j++)
+			if (i != j) {
+				float temp = x[i][j];
+				x[i][j] = x[j][i];
+				x[j][i] = temp;
+			}
+}
+
+// Transpose the given matrix using the transpose routine above.
+Matrix4x4
+Transpose(const Matrix4x4& M)
+{
+	Matrix4x4 TransposeMe = M;
+	TransposeMe.Transpose();
+	return TransposeMe;
+}
+
+// Construct an identity matrix.
+Matrix4x4
+IdentityMatrix()
+{
+	Matrix4x4 M;
+
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			M.x[i][j] = (i == j) ? 1.0f : 0.0f;
+	return M;
+}
+
+// Construct a zero matrix.
+Matrix4x4
+ZeroMatrix()
+{
+	Matrix4x4 M;
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			M.x[i][j] = 0;
+	return M;
+}
+
+// Construct a translation matrix given the location to translate to.
+Matrix4x4
+TranslationMatrix(const Vector3& Location)
+{
+	Matrix4x4 M = IdentityMatrix();
+	M.x[3][0] = Location.x;
+	M.x[3][1] = Location.y;
+	M.x[3][2] = Location.z;
+	return M;
+}
+
+// Construct a rotation matrix.  Rotates Angle radians about the
+// X axis.
+Matrix4x4
+RotationXMatrix(float Angle)
+{
+	Matrix4x4 M = IdentityMatrix();
+	float Cosine = cos(Angle);
+	float Sine = sin(Angle);
+	M.x[1][1] = Cosine;
+	M.x[2][1] = -Sine;
+	M.x[1][2] = Sine;
+	M.x[2][2] = Cosine;
+	return M;
+}
+
+// Construct a rotation matrix.  Rotates Angle radians about the
+// Y axis.
+Matrix4x4 RotationYMatrix(float angle)
+{
+	Matrix4x4 m = IdentityMatrix();
+
+	float cosine = cos(angle);
+	float sine = sin(angle);
+
+	m.x[0][0] = cosine;
+	m.x[2][0] = -sine;
+	m.x[0][2] = sine;
+	m.x[2][2] = cosine;
+
+	return m;
+}
+
+
+// Construct a rotation matrix.  Rotates Angle radians about the
+// Z axis.
+Matrix4x4 RotationZMatrix(float angle)
+{
+	Matrix4x4 m = IdentityMatrix();
+
+	float cosine = cos(angle);
+	float sine = sin(angle);
+	
+	m.x[0][0] = cosine;
+	m.x[1][0] = -sine;
+	m.x[0][1] = sine;
+	m.x[1][1] = cosine;
+
+	return m;
+}
+
+
+// Construct a yaw-pitch-roll rotation matrix.	Rotate Yaw
+// radians about the XY axis, rotate Pitch radians in the
+// plane defined by the Yaw rotation, and rotate Roll radians
+// about the axis defined by the previous two angles.
+Matrix4x4 RotationYPRMatrix(float Yaw, float Pitch, float Roll)
+{
+	Matrix4x4 M;
+	float ch = cos(Yaw);
+	float sh = sin(Yaw);
+	float cp = cos(Pitch);
+	float sp = sin(Pitch);
+	float cr = cos(Roll);
+	float sr = sin(Roll);
+
+	M.x[0][0] = ch * cr + sh * sp * sr;
+	M.x[1][0] = -ch * sr + sh * sp * cr;
+	M.x[2][0] = sh * cp;
+	M.x[0][1] = sr * cp;
+	M.x[1][1] = cr * cp;
+	M.x[2][1] = -sp;
+	M.x[0][2] = -sh * cr - ch * sp * sr;
+	M.x[1][2] = sr * sh + ch * sp * cr;
+	M.x[2][2] = ch * cp;
+	for (int i = 0; i < 4; i++)
+		M.x[3][i] = M.x[i][3] = 0;
+	M.x[3][3] = 1;
+
+	return M;
+}
+
+// Construct a rotation of a given angle about a given axis.
+// Derived from Eric Haines's SPD (Standard Procedural 
+// Database).
+Matrix4x4 RotationAxisMatrix(const Vector3& axis, float angle)
+{
+	Matrix4x4 M;
+
+	float cosine = cos(angle);
+	float sine = sin(angle);
+	float one_minus_cosine = 1 - cosine;
+
+	M.x[0][0] = axis.x * axis.x + (1.0f - axis.x * axis.x) * cosine;
+	M.x[0][1] = axis.x * axis.y * one_minus_cosine + axis.z * sine;
+	M.x[0][2] = axis.x * axis.z * one_minus_cosine - axis.y * sine;
+	M.x[0][3] = 0;
+
+	M.x[1][0] = axis.x * axis.y * one_minus_cosine - axis.z * sine;
+	M.x[1][1] = axis.y * axis.y + (1.0f - axis.y * axis.y) * cosine;
+	M.x[1][2] = axis.y * axis.z * one_minus_cosine + axis.x * sine;
+	M.x[1][3] = 0;
+
+	M.x[2][0] = axis.x * axis.z * one_minus_cosine + axis.y * sine;
+	M.x[2][1] = axis.y * axis.z * one_minus_cosine - axis.x * sine;
+	M.x[2][2] = axis.z * axis.z + (1.0f - axis.z * axis.z) * cosine;
+	M.x[2][3] = 0;
+
+	M.x[3][0] = 0;
+	M.x[3][1] = 0;
+	M.x[3][2] = 0;
+	M.x[3][3] = 1;
+
+	return M;
+}
+
+
+// Constructs the rotation matrix that rotates 'vec1' to 'vec2'
+Matrix4x4 RotationVectorsMatrix(const Vector3 &vecStart, const Vector3 &vecTo)
+{
+	Vector3 vec = CrossProd(vecStart, vecTo);
+
+	if (Magnitude(vec) > Limits::Small) {
+		// vector exist, compute angle
+		float angle = acos(DotProd(vecStart, vecTo));
+		// normalize for sure
+		vec.Normalize();
+		return RotationAxisMatrix(vec, angle);
+	}
+
+	// opposite or colinear vectors
+	Matrix4x4 ret = IdentityMatrix();
+	if (DotProd(vecStart, vecTo) < 0.0f)
+		ret *= -1.0f; // opposite vectors
+
+	return ret;
+}
+
+
+// Construct a scale matrix given the X, Y, and Z parameters
+// to scale by.  To scale uniformly, let X==Y==Z.
+Matrix4x4 ScaleMatrix(float X, float Y, float Z)
+{
+	Matrix4x4 M = IdentityMatrix();
+
+	M.x[0][0] = X;
+	M.x[1][1] = Y;
+	M.x[2][2] = Z;
+
+	return M;
+}
+
+
+// Construct a rotation matrix that makes the x, y, z axes
+// correspond to the vectors given.
+Matrix4x4 GenRotation(const Vector3 &x, const Vector3 &y, const Vector3 &z)
+{
+	Matrix4x4 M = IdentityMatrix();
+
+	//  x  y
+	M.x[0][0] = x.x;
+	M.x[1][0] = x.y;
+	M.x[2][0] = x.z;
+
+	M.x[0][1] = y.x;
+	M.x[1][1] = y.y;
+	M.x[2][1] = y.z;
+
+	M.x[0][2] = z.x;
+	M.x[1][2] = z.y;
+	M.x[2][2] = z.z;
+
+	return M;
+}
+
+// Construct a quadric matrix.	After Foley et al. pp. 528-529.
+Matrix4x4
+QuadricMatrix(float a, float b, float c, float d, float e,
+			  float f, float g, float h, float j, float k)
+{
+	Matrix4x4 M;
+
+	M.x[0][0] = a;  M.x[0][1] = d;  M.x[0][2] = f;  M.x[0][3] = g;
+	M.x[1][0] = d;  M.x[1][1] = b;  M.x[1][2] = e;  M.x[1][3] = h;
+	M.x[2][0] = f;  M.x[2][1] = e;  M.x[2][2] = c;  M.x[2][3] = j;
+	M.x[3][0] = g;  M.x[3][1] = h;  M.x[3][2] = j;  M.x[3][3] = k;
+
+	return M;
+}
+
+// Construct various "mirror" matrices, which flip coordinate
+// signs in the various axes specified.
+Matrix4x4
+MirrorX()
+{
+	Matrix4x4 M = IdentityMatrix();
+	M.x[0][0] = -1;
+	return M;
+}
+
+Matrix4x4
+MirrorY()
+{
+	Matrix4x4 M = IdentityMatrix();
+	M.x[1][1] = -1;
+	return M;
+}
+
+Matrix4x4
+MirrorZ()
+{
+	Matrix4x4 M = IdentityMatrix();
+	M.x[2][2] = -1;
+	return M;
+}
+
+Matrix4x4
+RotationOnly(const Matrix4x4& x)
+{
+	Matrix4x4 M = x;
+	M.x[3][0] = M.x[3][1] = M.x[3][2] = 0;
+	return M;
+}
+
+// Add corresponding elements of the two matrices.
+Matrix4x4&
+Matrix4x4::operator+= (const Matrix4x4& A)
+{
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			x[i][j] += A.x[i][j];
+	return *this;
+}
+
+// Subtract corresponding elements of the matrices.
+Matrix4x4&
+Matrix4x4::operator-= (const Matrix4x4& A)
+{
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			x[i][j] -= A.x[i][j];
+	return *this;
+}
+
+// Scale each element of the matrix by A.
+Matrix4x4&
+Matrix4x4::operator*= (float A)
+{
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			x[i][j] *= A;
+	return *this;
+}
+
+// Multiply two matrices.
+Matrix4x4&
+Matrix4x4::operator*= (const Matrix4x4& A)
+{
+	Matrix4x4 ret = *this;
+
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++) {
+			float subt = 0;
+			for (int k = 0; k < 4; k++)
+				subt += ret.x[i][k] * A.x[k][j];
+			x[i][j] = subt;
+		}
+		return *this;
+}
+
+// Add corresponding elements of the matrices.
+Matrix4x4
+operator+ (const Matrix4x4& A, const Matrix4x4& B)
+{
+	Matrix4x4 ret;
+
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			ret.x[i][j] = A.x[i][j] + B.x[i][j];
+	return ret;
+}
+
+// Subtract corresponding elements of the matrices.
+Matrix4x4
+operator- (const Matrix4x4& A, const Matrix4x4& B)
+{
+	Matrix4x4 ret;
+
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			ret.x[i][j] = A.x[i][j] - B.x[i][j];
+	return ret;
+}
+
+// Multiply matrices.
+Matrix4x4
+operator* (const Matrix4x4& A, const Matrix4x4& B)
+{
+	Matrix4x4 ret;
+
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++) {
+			float subt = 0;
+			for (int k = 0; k < 4; k++)
+				subt += A.x[i][k] * B.x[k][j];
+			ret.x[i][j] = subt;
+		}
+		return ret;
+}
+
+// Transform a vector by a matrix.
+Vector3
+operator* (const Matrix4x4& M, const Vector3& v)
+{
+	Vector3 ret;
+	float denom;
+
+	ret.x = v.x * M.x[0][0] + v.y * M.x[1][0] + v.z * M.x[2][0] + M.x[3][0];
+	ret.y = v.x * M.x[0][1] + v.y * M.x[1][1] + v.z * M.x[2][1] + M.x[3][1];
+	ret.z = v.x * M.x[0][2] + v.y * M.x[1][2] + v.z * M.x[2][2] + M.x[3][2];
+	denom = M.x[0][3] + M.x[1][3] + M.x[2][3] + M.x[3][3];
+	if (denom != 1.0)
+		ret /= denom;
+	return ret;
+}
+
+// Apply the rotation portion of a matrix to a vector.
+Vector3
+RotateOnly(const Matrix4x4& M, const Vector3& v)
+{
+	Vector3 ret;
+	float denom;
+
+	ret.x = v.x * M.x[0][0] + v.y * M.x[1][0] + v.z * M.x[2][0];
+	ret.y = v.x * M.x[0][1] + v.y * M.x[1][1] + v.z * M.x[2][1];
+	ret.z = v.x * M.x[0][2] + v.y * M.x[1][2] + v.z * M.x[2][2];
+	denom = M.x[0][3] + M.x[1][3] + M.x[2][3] + M.x[3][3];
+	if (denom != 1.0)
+		ret /= denom;
+	return ret;
+}
+
+// Scale each element of the matrix by B.
+Matrix4x4
+operator* (const Matrix4x4& A, float B)
+{
+	Matrix4x4 ret;
+
+	for (int i = 0; i < 4; i++)
+		for (int j = 0; j < 4; j++)
+			ret.x[i][j] = A.x[i][j] * B;
+	return ret;
+}
+
+// Overloaded << for C++-style output.
+ostream&
+operator<< (ostream& s, const Matrix4x4& M)
+{
+	for (int i = 0; i < 4; i++) { // y
+		for (int j = 0; j < 4; j++) { // x
+			//  x  y 
+			s << setprecision(4) << setw(10) << M.x[j][i];
+		}
+		s << '\n';
+	}
+	return s;
+}
+
+// Rotate a direction vector...
+Vector3
+PlaneRotate(const Matrix4x4& tform, const Vector3& p)
+{
+	// I sure hope that matrix is invertible...
+	Matrix4x4 use = Transpose(Invert(tform));
+
+	return RotateOnly(use, p);
+}
+
+// Transform a normal
+Vector3
+TransformNormal(const Matrix4x4& tform, const Vector3& n)
+{
+	Matrix4x4 use = NormalTransformMatrix(tform);
+
+	return RotateOnly(use, n);
+}
+
+Matrix4x4 
+NormalTransformMatrix(const Matrix4x4 &tform)
+{
+	Matrix4x4 m = tform;
+	// for normal translation vector must be zero!
+	m.x[3][0] = m.x[3][1] = m.x[3][2] = 0.0; 
+	// I sure hope that matrix is invertible...
+	return Transpose(Invert(m));
+}
+
+
+Vector3 GetTranslation(const Matrix4x4 &M)
+{
+	return Vector3(M.x[3][0], M.x[3][1], M.x[3][2]);
+}
+
+//}
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Timers.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Timers.h	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Timers.h	(revision 2746)
@@ -1,11 +1,4 @@
-// ================================================================
-// $Id: Timers.h,v, 1.1.1.1 2004/10/12 12:49:13 matt Exp $
-// ****************************************************************
-// 
-/** \file time.h
-
-Various time measuring routines
-
-@author Jiri Bittner
+/** Various time measuring routines
+    @author Jiri Bittner
 */
 
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Matrix4x4.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Matrix4x4.h	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Matrix4x4.h	(revision 2746)
@@ -0,0 +1,154 @@
+#ifndef __MATRIX4x4_H
+#define __MATRIX4x4_H
+
+#include <iostream>
+
+
+//namespace GtpVisibilityPreprocessor {
+
+class Vector3;
+
+
+
+class Matrix4x4
+{
+public:
+	float x[4][4]; // first index is column [x], the second is row [y]
+
+	Matrix4x4() { }
+
+	// here xXY - 'X' is row, and 'Y' is column - classical mathematical notation
+	Matrix4x4(float x11, float x12, float x13, float x14,
+		float x21, float x22, float x23, float x24,
+		float x31, float x32, float x33, float x34,
+		float x41, float x42, float x43, float x44);
+
+	Matrix4x4(const Vector3 &a, const Vector3 &b, const Vector3 &c);
+
+
+	// Assignment operators
+	Matrix4x4& operator+= (const Matrix4x4 &A); 	// add-to
+	Matrix4x4& operator-= (const Matrix4x4 &A); 	// subtract-from
+	Matrix4x4& operator*= (const Matrix4x4 &A); 	// multiply by matrix
+	Matrix4x4& operator*= (float A); 		// scale by scalar
+
+	// Fundamental operations
+	/// Invert the matrix .. returns 0 = regular
+	int Invert();		
+	/// Transpose the matrix
+	void Transpose();	
+
+	void
+		SetColumns(const Vector3 &a, const Vector3 &b, const Vector3 &c);
+
+	float Det3x3() const {
+		return (x[0][0]*x[1][1]*x[2][2] + \
+			x[1][0]*x[2][1]*x[0][2] + \
+			x[2][0]*x[0][1]*x[1][2] - \
+			x[2][0]*x[1][1]*x[0][2] - \
+			x[0][0]*x[2][1]*x[1][2] - \
+			x[1][0]*x[0][1]*x[2][2]);
+	}
+
+
+	friend Matrix4x4 Invert(const Matrix4x4 &M);	// Invert a given matrix
+	friend Matrix4x4 Transpose(const Matrix4x4 &M);	// Transpose a given matrix
+
+	// Create various types of matrix.
+	friend Matrix4x4 IdentityMatrix();
+	friend Matrix4x4 ZeroMatrix();
+	friend Matrix4x4 TranslationMatrix(const Vector3 &Location);
+	friend Matrix4x4 RotationXMatrix(float Angle);
+	friend Matrix4x4 RotationYMatrix(float Angle);
+	friend Matrix4x4 RotationZMatrix(float Angle);
+	friend Matrix4x4 RotationYPRMatrix(float Yaw, float Pitch, float Roll);
+	// about axis 'axis' by angle 'Angle'
+	friend Matrix4x4 RotationAxisMatrix(const Vector3 &axis, float Angle);
+	// create the rotation matrix that rotates 'vecFrom' to 'vecTo'
+	friend Matrix4x4 RotationVectorsMatrix(const Vector3 &vecFrom,
+		const Vector3 &vecTo);
+
+	friend Matrix4x4 ScaleMatrix(float X, float Y, float Z);
+	friend Matrix4x4 GenRotation(const Vector3 &x, const Vector3 &y,
+		const Vector3 &z);
+	friend Matrix4x4 QuadricMatrix(float a, float b, float c, float d, float e,
+		float f, float g, float h, float j, float k);
+	// returns matrix for transforming normal
+	friend Matrix4x4 NormalTransformMatrix(const Matrix4x4 &M);
+
+	friend Matrix4x4 MirrorX();
+	friend Matrix4x4 MirrorY();
+	friend Matrix4x4 MirrorZ();
+	friend Matrix4x4 RotationOnly(const Matrix4x4 &x);
+
+	// Binary operators
+	friend Matrix4x4 operator+ (const Matrix4x4 &A, const Matrix4x4 &B);
+	friend Matrix4x4 operator- (const Matrix4x4 &A, const Matrix4x4 &B);
+	friend Matrix4x4 operator* (const Matrix4x4 &A, float B);
+	friend Matrix4x4 operator* (const Matrix4x4 &A, const Matrix4x4 &B);
+
+	// friends returning Vector3
+	friend Vector3 operator*(const Matrix4x4 &M, const Vector3 &v);
+	friend Vector3 RotateOnly(const Matrix4x4 &M, const Vector3 &v);
+	friend Vector3 TransformNormal(const Matrix4x4 &M, const Vector3 &v);
+	friend Vector3 GetTranslation(const Matrix4x4 &M);
+
+	// Construct rotation description according VRML'97 specification
+	//  const CVector4D SFRotation(void) const;
+
+
+
+	// Overloaded output operator.
+	friend std::ostream& operator<< (std::ostream &s, const Matrix4x4 &M);
+};
+
+// forward declaration
+Matrix4x4 IdentityMatrix();
+Matrix4x4 Invert(const Matrix4x4 &M);	// Invert a given matrix
+Matrix4x4 Transpose(const Matrix4x4 &M);	// Transpose a given matrix
+
+// Create various types of matrix.
+Matrix4x4 IdentityMatrix();
+Matrix4x4 ZeroMatrix();
+Matrix4x4 TranslationMatrix(const Vector3 &Location);
+Matrix4x4 RotationXMatrix(float Angle);
+Matrix4x4 RotationYMatrix(float Angle);
+Matrix4x4 RotationZMatrix(float Angle);
+Matrix4x4 RotationYPRMatrix(float Yaw, float Pitch, float Roll);
+// about axis 'axis' by angle 'Angle'
+Matrix4x4 RotationAxisMatrix(const Vector3 &axis, float Angle);
+// create the rotation matrix that rotates 'vecFrom' to 'vecTo'
+Matrix4x4 RotationVectorsMatrix(const Vector3 &vecFrom,
+								const Vector3 &vecTo);
+
+Matrix4x4 ScaleMatrix(float X, float Y, float Z);
+Matrix4x4 GenRotation(const Vector3 &x, const Vector3 &y,
+					  const Vector3 &z);
+Matrix4x4 QuadricMatrix(float a, float b, float c, float d, float e,
+						float f, float g, float h, float j, float k);
+// returns matrix for transforming normal
+Matrix4x4 NormalTransformMatrix(const Matrix4x4 &M);
+
+Matrix4x4 MirrorX();
+Matrix4x4 MirrorY();
+Matrix4x4 MirrorZ();
+Matrix4x4 RotationOnly(const Matrix4x4 &x);
+
+// Binary operators
+Matrix4x4 operator+ (const Matrix4x4 &A, const Matrix4x4 &B);
+Matrix4x4 operator- (const Matrix4x4 &A, const Matrix4x4 &B);
+Matrix4x4 operator* (const Matrix4x4 &A, float B);
+Matrix4x4 operator* (const Matrix4x4 &A, const Matrix4x4 &B);
+
+// friends returning Vector3
+Vector3 operator*(const Matrix4x4 &M, const Vector3 &v);
+Vector3 RotateOnly(const Matrix4x4 &M, const Vector3 &v);
+Vector3 TransformNormal(const Matrix4x4 &M, const Vector3 &v);
+Vector3 GetTranslation(const Matrix4x4 &M);
+
+
+//}
+
+
+
+#endif
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/RenderTraverser.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/RenderTraverser.cpp	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/RenderTraverser.cpp	(revision 2746)
@@ -3,9 +3,6 @@
 #include "Timers.h"
 
-/*extern "C"
-{
-	#include "MathStuff.h"
-}*/
-
+
+#if TOIMPLEMENT
 
 RenderTraverser::RenderTraverser(): 
@@ -546,2 +543,5 @@
 	printf("using opt %d\n", mUseOptimization);
 }
+
+
+#endif
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.cpp	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.cpp	(revision 2746)
@@ -0,0 +1,2592 @@
+#if TOIMPLEMENT
+
+#include <queue>
+#include <stack>
+#include "Bvh.h"
+#include "yare.h"
+#include "PerfTimer.h"
+#include "Camera.h"
+#include "Settings.h"
+#include "Context.h"
+#include "TriangleBvh.h"
+#include "NodeGeometry.h"
+#include "Viewer.h"
+
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+#define INVALID_TEST ((unsigned int)-1)
+#define MAX_FLOAT 1e20f
+
+#define TYPE_INTERIOR -2
+#define TYPE_LEAF -3
+#define USE_TIGHTER_BOUNDS_FOR_ALL 0
+
+const static bool useVbos = true;
+unsigned int Bvh::sCurrentVboId = -1;
+
+using namespace std;
+
+
+
+/*
+3 x---------x 2
+  |\         \
+  | \         \
+  |  \         \
+  |	4 x---------x 5
+  |   |         | 
+0 x   |     x 1 |
+   \  |         |
+    \ |         |
+     \|         |
+  	7 x---------x 6             
+*/
+
+static unsigned int sIndices[] =	
+{7, // degenerated
+ 7, 6, 4, 5, 3, 2, 0, 1,
+ 1, 4, // degenerated
+ 4, 3, 7, 0, 6, 1, 5, 2,
+ 2 // degenerated
+};
+
+
+const static int sNumIndicesPerBox = 20;
+
+/*	Order of vertices
+	0 = (0, 0, 0)
+	1 = (1, 0, 0)
+	2 = (1, 1, 0)
+	3 = (0, 1, 0)
+	4 = (0, 1, 1)
+	5 = (1, 1, 1)
+	6 = (1, 0, 1)
+	7 = (0, 0, 1)
+*/
+
+
+BvhNode::BvhNode(BvhNode *parent): 
+mParent(parent),
+mAxis(-1), 
+mDepth(0), 
+mPlaneMask(0),
+mPreferredPlane(0),
+mVizBox(NULL),
+mLastRenderedFrame(-999),
+mFirst(-1), 
+mLast(-1),
+mNumTestNodes(1),
+mTestNodesIdx(-1),
+mIndicesPtr(-1)
+{
+}
+
+
+void BvhNode::ResetVisibility()
+{
+	mVisibility.Reset();
+	mMeasurements.Reset();
+
+	mLastRenderedFrame = -999;
+}
+
+
+BvhNode::~BvhNode() 
+{
+	if (mVizBox) mVizBox->unref();
+}
+
+
+Box *BvhNode::GetOrCreateVizBox() 
+{
+	if (!mVizBox) 
+	{
+		mVizBox = new Box(mBox); 
+		mVizBox->ref();
+	}
+
+	return mVizBox;
+}
+
+
+void BvhNode::VisibilityInfo::Reset()
+{
+	mIsVisible = false;
+	mIsFullyVisible = false;
+
+	mVisiblePixels = 0;
+	mAssumedVisibleFrames = 0;
+	mRemainingVisibleFrames = 0;
+
+	mLastVisitedFrame = -1;
+	mLastTestedFrame = -1;
+	mTurnedVisibleFrame = 0;
+
+	mTimesInvisible = 0;
+	mTimesTested = 0;
+	mTimesChangedClassification = 0;
+	mAvgChangedClassification = 0;
+	mIsFrustumCulled = false;
+
+	mLastTestedVisibleFrame = 0;
+	mIsNew = true;
+}
+
+
+void BvhNode::SetFullyVisible(const bool visible)
+{
+	mVisibility.mIsFullyVisible = visible;
+}
+
+
+BvhLeaf::~BvhLeaf()
+{
+	if (mTriangleBvh) delete mTriangleBvh;
+}
+
+
+/***********************************************************/
+/*                class Bvh implementation                 */
+/***********************************************************/
+
+
+
+Bvh::Bvh(const GeometryVector &geometry, 
+		 DistanceSortedRenderAction *const renderer):
+mCamera(NULL), 
+mFrameId(-1), 
+mRoot(NULL), 
+mVertices(NULL),
+mIndices(NULL),
+mUseTighterBoundsOnlyForLeafTests(false),
+mRenderer(renderer),
+mTestIndices(NULL),
+mCurrentIndicesPtr(0),
+mCollectTighterBoundsWithMaxLevel(true)
+{
+	mGeometry = new NodeGeometry*[geometry.size()];
+	mGeometrySize = geometry.size();
+
+	BoundingBox sceneBox;
+
+	// compute scene extent
+	for (size_t i = 0; i < mGeometrySize; ++ i)
+	{
+		mGeometry[i] = geometry[i];
+		sceneBox.combine(&mGeometry[i]->mBox);
+
+		mBvhStats.mTriangles += mGeometry[i]->mNumTriangles;
+	}
+
+
+	////////
+	//-- create root
+
+	BvhLeaf *leaf = new BvhLeaf(NULL);
+	leaf->mDepth = 0;
+	leaf->mFirst = 0;
+	leaf->mLast = (int)geometry.size() - 1;
+	leaf->mBox = sceneBox;
+
+	OUT1("geometry in root: " << leaf->mLast);
+
+	mRoot = leaf;
+	mNumNodes = 1;
+
+	// parameters
+	mMaxGeometry = Settings::Global()->get_nvocc_bvh_max_objects();
+	mMaxTriangles = Settings::Global()->get_nvocc_bvh_max_triangles();
+	mMaxDepth = Settings::Global()->get_nvocc_bvh_max_depth();
+	mSplitType = Settings::Global()->get_nvocc_bvh_split_type();
+
+	float sceneArea = sceneBox.getSurface();
+	float minAreaRatio = Settings::Global()->get_nvocc_bvh_min_area();
+
+	mMinArea = minAreaRatio * sceneArea;
+
+	mMaxDepthForTestingChildren = Settings::Global()->get_nvocc_bvh_max_depth_for_testing_children();
+	mUseTighterBoundsOnlyForLeafTests = (Settings::Global()->get_nvocc_use_tighter_bounds() >= 1);
+	
+	mAreaRatioThreshold = Settings::Global()->get_nvocc_area_ratio_threshold();
+	mVolRatioThreshold = Settings::Global()->get_nvocc_vol_ratio_threshold();
+}
+
+
+Bvh::~Bvh() 
+{
+	if (mVertices) delete []mVertices; 
+	if (mIndices) delete [] mIndices;
+	if (mTestIndices) delete [] mTestIndices;
+	if (mGeometry) delete [] mGeometry;
+
+	if (mRoot) delete mRoot;
+}
+
+
+void Bvh::Construct()
+{
+	PerfTimer constructTimer;
+	constructTimer.Entry();
+
+	OUT1("Info: Constructing BVH");
+
+	mRoot = SubdivideLeaf((BvhLeaf *)mRoot, 0);
+
+	constructTimer.Exit();
+
+	OUT1("Info: Bvh done in " << constructTimer.TotalCount() << " seconds.");
+	OUT1("Info: Number Of BVH Nodes = " << mNumNodes);
+
+	// update stats on leaf nodes
+	UpdateNumLeaves(mRoot);
+	
+	mAvgDepth = 1.0f + log((float)GetNumLeaves()) / log(2.0f);
+
+	BvhLeafContainer leaves;
+	leaves.reserve(GetNumLeaves());
+	CollectLeaves(mRoot, leaves);
+
+	// compute tighter boundaries for the leaves
+	PostProcessLeaves(leaves);
+	// compute unique ids
+	ComputeIds();
+	/// pull up some data from the leaves
+	UpdateInteriors(mRoot);
+	// recompute the boundaries once
+	RecomputeBounds();
+	// do stats
+	ComputeBvhStats();
+
+	PrintBvhStats();
+}
+
+
+float Bvh::EvaluateSahCost(BvhLeaf *leaf, const int axis, const float position)
+{
+	// count triangles on the left / right
+	int left = 0;
+	int right = 0;
+	BoundingBox leftBox, rightBox;
+
+	const int candidates = std::max(50, (int)(leaf->CountPrimitives() / 20));
+
+	const float finc = leaf->CountPrimitives() / (float)candidates;
+
+	int i = leaf->mFirst;
+	float fi = leaf->mFirst + 0.5f;
+
+	BoundingBox box;
+
+	for (; i <= leaf->mLast;) 
+	{
+		box = mGeometry[i]->GetBoundingBox();
+
+		if (box.getCenter()[axis] < position) 
+		{
+			++ left;
+			// update the box
+			leftBox.combine(&box);
+		} 
+		else 
+		{
+			++ right;
+			rightBox.combine(&box);
+		}
+		 
+		fi += finc;
+		i = (int)fi;
+	}
+
+	float bW = 1.0f;
+	float leftRatio = left / (float)leaf->CountPrimitives();
+	float rightRatio = right / (float)leaf->CountPrimitives();
+
+	float saLeft = 0.0f;
+	float saRight = 0.0f;
+
+	// not a valid split
+	if (!left || !right)
+		return 1e25f;
+
+	saLeft = leftBox.getSurface();
+	saRight = rightBox.getSurface();
+
+
+	return
+		saLeft  * ((1.0f - bW) + bW * leftRatio) + 
+		saRight * ((1.0f - bW) + bW * rightRatio);
+}
+
+
+float Bvh::SelectPlaneSah(BvhLeaf *leaf, int &axis, float &minCost)
+{
+	minCost = MAX_FLOAT;
+	float bestPos = minCost;
+	int bestAxis = 0;
+
+	//  cout<<"Evaluating axis..."<<endl<<flush;
+
+	const int initialPlanes = 3;
+
+	// initiate the costs
+	for (axis = 0; axis < 3; ++ axis) 
+	{
+		const float size = leaf->mBox.getUpper()[axis] - leaf->mBox.getLower()[axis];
+
+		for (int i = 0; i < initialPlanes; ++ i) 
+		{
+			const float pos = leaf->mBox.getLower()[axis] + (i + 1) * size / (initialPlanes + 1);
+			const float cost = EvaluateSahCost(leaf, axis, pos);
+
+			if (cost < minCost) 
+			{
+				minCost = cost;
+				bestPos = pos;
+				bestAxis = axis;
+			}
+		}
+	}
+
+	axis = bestAxis;
+
+	//  cout<<axis<<endl<<flush;
+	const float shrink = 0.5f;
+
+	// now hierarchically refine the initial interval
+	float size = shrink * 
+		(leaf->mBox.getUpper()[axis] - leaf->mBox.getLower()[axis]) / (initialPlanes + 1);
+
+	const int steps = 4;
+	float cost;
+
+	for (int i = 0; i < steps; ++ i) 
+	{
+		const float left = bestPos - size;
+		const float right = bestPos + size;
+
+		cost = EvaluateSahCost(leaf, axis, left);
+
+		if (cost < minCost) 
+		{
+			minCost = cost;
+			bestPos = left;
+		}
+
+		cost = EvaluateSahCost(leaf, axis, right);
+
+		if (cost < minCost) 
+		{
+			minCost = cost;
+			bestPos = right;
+		}
+
+		size = shrink * size;
+	}
+
+	//OUT1("best axis: " << axis << " " << bestPos << " " << minCost);
+
+	return bestPos;
+}
+
+
+void Bvh::ResizeVisibilityBuffers(const int maxSize)
+{
+	std::stack<BvhNode *> nodeStack;
+
+	nodeStack.push(mRoot);
+
+	while (!nodeStack.empty()) 
+	{
+		BvhNode *node = nodeStack.top();
+		nodeStack.pop();
+
+		static_cast<BvhLeaf *>(node)->mMeasurements.Reset();
+
+		if (!node->IsLeaf()) 
+		{
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+			
+			nodeStack.push(interior->mBack);
+			nodeStack.push(interior->mFront);
+		}
+	}
+}
+
+
+int Bvh::SortTriangles(BvhLeaf *leaf, 
+					   const int axis, 
+					   const float position)
+{
+	int i = leaf->mFirst;
+	int j = leaf->mLast;
+
+    while (1) 
+	{
+		while (mGeometry[i]->GetBoundingBox().getCenter()[axis] < position) 
+			++ i;
+	
+		while (position < mGeometry[j]->GetBoundingBox().getCenter()[axis])
+			-- j;
+
+		if (i < j) 
+		{
+			std::swap(mGeometry[i], mGeometry[j]);
+			
+			++ i;
+			-- j;
+		}
+		else
+		{
+			break;
+		}
+	}
+
+	return j;
+}
+
+
+int Bvh::SortTrianglesSpatialMedian(BvhLeaf *leaf, const int axis)
+{
+	// spatial median
+	float x = leaf->mBox.getCenter()[axis];
+
+	return SortTriangles(leaf, axis, x);
+}
+
+
+int Bvh::SortTrianglesObjectMedian(BvhLeaf *leaf, const int axis, float &pos)
+{
+	// Now distribute the objects into the subnodes
+	// Use QuickMedian sort
+	int l = leaf->mFirst;
+	int r = leaf->mLast;
+	int k = (l + r) / 2;
+
+	while (l < r) 
+	{
+		int i = l;
+		int j = r;
+
+		// get some estimation of the pivot
+		pos = mGeometry[(l + r) / 2]->GetBoundingBox().getCenter()[axis];
+
+		while (1) 
+		{
+			while ((i  <= leaf->mLast) && (mGeometry[i]->GetBoundingBox().getCenter()[axis] < pos))
+				++ i;
+
+			while((j >= leaf->mFirst) && (pos < mGeometry[j]->GetBoundingBox().getCenter()[axis]))
+				-- j;
+
+			if (i <= j) 
+			{
+				std::swap(mGeometry[i], mGeometry[j]);
+				++ i;
+				-- j;
+			}
+			else
+			{
+				break;
+			}
+		}
+
+		// now check the extents
+		if (i >= k)
+			r = j;
+		else
+			l = i;
+	}
+
+	return k;
+}
+
+
+int Bvh::SortTrianglesSurfaceArea(BvhLeaf *leaf, const float sa)
+{
+	int i = leaf->mFirst;
+	int j = leaf->mLast;
+
+	while(1) 
+	{
+		while ((i <= j) && (mGeometry[i]->GetBoundingBox().getSurface() < sa)) 
+			++ i;
+
+		while ((i <= j) && (sa < mGeometry[j]->GetBoundingBox().getSurface()))
+			-- j;
+
+		if (i < j) 
+		{
+			swap(mGeometry[i], mGeometry[j]);
+			++ i;
+			-- j;
+		}
+		else
+			break;
+	}
+
+	return j;
+}
+
+
+bool Bvh::TerminationCriteriaMet(BvhLeaf *leaf) const
+{
+	const bool terminationCriteriaMet = 
+		(leaf->mBox.getSurface() < mMinArea) ||
+		(leaf->CountPrimitives() <= mMaxGeometry) ||
+		(CountTriangles(leaf) <= mMaxTriangles) ||
+		(leaf->mDepth > mMaxDepth);
+
+	return terminationCriteriaMet;
+}
+
+
+BvhNode *Bvh::SubdivideLeaf(BvhLeaf *leaf, const int parentAxis)
+{
+	if (TerminationCriteriaMet(leaf))
+		return leaf;
+	
+	//int axis = leaf->mBox.MajorAxis();
+	int axis = (parentAxis + 1) % 3;
+
+	// position of the split in the partailly sorted array of triangles
+	// corresponding to this leaf
+	int split = -1;
+	const int scale = 20;
+
+	float pos = -1.0f;
+
+	// Spatial median subdivision
+	switch (mSplitType) 
+	{
+	case SPATIAL_MEDIAN:
+		split = SortTrianglesSpatialMedian(leaf, axis);
+
+		if ( 
+			((split - leaf->mFirst) < leaf->CountPrimitives() / scale) ||
+			((leaf->mLast - split) < leaf->CountPrimitives() / scale) ) 
+		{
+			split = SortTrianglesObjectMedian(leaf, axis, pos);
+		}
+
+		pos = leaf->mBox.getCenter()[axis];
+		break;
+	case OBJECT_MEDIAN:
+		// Object median subdivision: approximately the same number
+		// of objects on the left of the the splitting point.
+		split = SortTrianglesObjectMedian(leaf, axis, pos);
+		break;
+	case SAH:
+		{
+			float cost;
+			pos = SelectPlaneSah(leaf, axis, cost);
+
+			if (pos != MAX_FLOAT)
+			{
+				split = SortTriangles(leaf, axis, pos);
+			}
+
+			if ((pos == MAX_FLOAT) || (split == leaf->mLast))
+			{
+				split = -1;
+				split = SortTrianglesObjectMedian(leaf, axis, pos);
+			}
+		}
+		break;
+	case SAH_OR_SIZE: 
+		{
+			// split by size instead
+			const float saThreshold = 0.2f * leaf->GetBox().getSurface();
+			split = SortTrianglesSurfaceArea(leaf, saThreshold);
+
+			if ((split == leaf->mLast) || (split == leaf->mFirst - 1)) 
+			{
+				// use SAH
+				float cost;
+				pos = SelectPlaneSah(leaf, axis, cost);
+
+				if (pos != MAX_FLOAT)
+					split = SortTriangles(leaf, axis, pos);
+				else
+					split = SortTrianglesObjectMedian(leaf, axis, pos);
+			}
+			else
+			{
+				// note: no position is computed!!
+				//OUT1("sorted by size");
+			}
+		}
+		break;
+	default:
+		OUT1("should not come here");
+		break;
+	}
+
+	if (1 && ((split == leaf->mLast)))
+	{
+		// no reduction: we should never come here
+		OUT1("error: no reduction " << leaf->CountPrimitives() << " " << leaf->mFirst << " " << split << " " << leaf->mLast);
+		return leaf;
+	}
+
+	// create two more nodes
+	mNumNodes += 2;
+	BvhInterior *parent = new BvhInterior(leaf->GetParent());
+	BvhLeaf *front = new BvhLeaf(parent);
+	
+	parent->mAxis = axis;
+	parent->mBox = leaf->mBox;
+	parent->mDepth = leaf->mDepth;
+
+	parent->mBack = leaf;
+	parent->mFront = front;
+	parent->mPosition = pos;
+	
+	// now assign the triangles to the subnodes
+	front->mFirst = split + 1;
+	front->mLast = leaf->mLast;
+	front->mDepth = leaf->mDepth + 1;
+	leaf->mLast = split;
+	leaf->mDepth = front->mDepth;
+	leaf->mParent = parent;
+	
+	// reset box
+	leaf->mBox = BoundingBox();
+
+	UpdateLeafBox(static_cast<BvhLeaf *>(parent->mBack));
+	UpdateLeafBox(static_cast<BvhLeaf *>(parent->mFront));
+
+	// some output
+	const int n = 500;
+	if ((GetNumLeaves() % n) == n - 1)
+		OUT1("created " << GetNumLeaves() << " leaves");
+	
+#if VERBOSE_OUTPUT
+	OUT1("box: " << parent->mBox.getUpper() << " " << parent->mBox.getLower());
+	OUT1("depth: " << (int)parent->mDepth);
+	OUT1("axis: " << axis);
+	OUT1("bc="<<((BvhLeaf *)parent->mBack)->Count());
+	OUT1("fc="<<((BvhLeaf *)parent->mFront)->Count());
+	OUT1("back: " << parent->mBack->mBox.getSurface());
+	OUT1("front: " << parent->mFront->mBox.getSurface());
+#endif
+
+	// recursively continue subdivision
+	parent->mBack = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mBack), axis);
+	parent->mFront = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mFront), axis);
+
+	return parent;
+}
+
+
+void Bvh::UpdateLeafBox(BvhLeaf *leaf)
+{
+	for (int i = leaf->mFirst; i <= leaf->mLast; ++ i) 
+	{
+		leaf->mBox.combine(&mGeometry[i]->mBox);
+	} // for
+}
+
+
+void Bvh::UpdateNumLeaves(BvhNode *node) const
+{
+	if (node->IsLeaf())
+	{
+		node->mNumLeaves = 1;
+	}
+	else
+	{
+		BvhNode *f = static_cast<BvhInterior *>(node)->GetFront();
+		BvhNode *b = static_cast<BvhInterior *>(node)->GetBack();
+
+		UpdateNumLeaves(f);
+		UpdateNumLeaves(b);
+
+		node->mNumLeaves = f->mNumLeaves + b->mNumLeaves;
+	}
+}
+
+
+void Bvh::UpdateBoxes(BvhNode *node)
+{
+	if (node->IsLeaf()) 
+	{
+		UpdateLeafBox(static_cast<BvhLeaf *>(node));
+	}
+	else 
+	{
+		BvhInterior *interior = static_cast<BvhInterior *>(node);
+	
+		UpdateBoxes(interior->mBack);
+		UpdateBoxes(interior->mFront);
+
+		node->mBox = interior->mBack->mBox;
+		node->mBox.combine(&interior->mFront->mBox);
+	}
+}
+
+
+void Bvh::UpdateFullVisibility(BvhNode *node) const
+{		
+	// node not visited in this frame => no change
+	if (node->GetLastVisitedFrame() != mFrameId) 
+		return;
+
+	// leaf node: terminate recursion
+	if (node->IsLeaf())
+	{
+		node->SetFullyVisible(node->IsVisible());	
+		return;
+	}
+
+	BvhInterior *interior = static_cast<BvhInterior *>(node);
+	
+	// recursive traversal
+	UpdateFullVisibility(interior->mBack);
+	UpdateFullVisibility(interior->mFront);
+	
+	interior->SetFullyVisible(interior->mBack->IsFullyVisible() && 
+							  interior->mFront->IsFullyVisible());
+}
+
+
+void Bvh::PullUpLastVisited(BvhNode *node, const int frameId) const
+{		
+	BvhNode *parent = node->GetParent();
+
+	while (parent && (parent->GetLastVisitedFrame() != frameId))
+	{
+		parent->SetLastVisitedFrame(frameId);
+		parent = parent->GetParent();
+	}
+}
+
+
+void Bvh::MakeParentsVisible(BvhNode *node)
+{
+	BvhNode *parent = node->GetParent();
+
+	while (parent && (!parent->IsVisible()))
+	{
+		parent->SetVisible(true);
+		parent = parent->GetParent();
+	}
+}
+
+
+void Bvh::CollectLeaves(BvhNode *node, BvhLeafContainer &leaves)
+{
+	stack<BvhNode *> tStack;
+	tStack.push(node);
+
+	while (!tStack.empty())
+	{
+		BvhNode *node = tStack.top();
+
+		tStack.pop();
+
+		if (!node->IsLeaf())
+		{
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+
+			tStack.push(interior->mFront);
+			tStack.push(interior->mBack);
+		}
+		else
+		{
+			leaves.push_back(static_cast<BvhLeaf *>(node));
+		}
+	}
+}
+
+
+void Bvh::CollectNodes(BvhNode *node, BvhNodeContainer &nodes)
+{
+	stack<BvhNode *> tStack;
+
+	tStack.push(node);
+
+	while (!tStack.empty())
+	{
+		BvhNode *node = tStack.top();
+		tStack.pop();
+
+		nodes.push_back(node);
+
+		if (!node->IsLeaf())
+		{
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+
+			tStack.push(interior->mFront);
+			tStack.push(interior->mBack);
+		}
+	}
+}
+
+
+typedef pair<BvhNode *, int> tPair;
+
+void Bvh::CollectNodes(BvhNode *root, 
+					   const int depth, 
+					   HierarchyNodeContainer &nodes)
+{
+	stack<tPair> tStack;
+	tStack.push(tPair(root, 0));
+
+	while (!tStack.empty())
+	{
+		BvhNode *node = tStack.top().first;
+		const int d = tStack.top().second;
+
+		tStack.pop();
+
+		// found depth => take this node
+		if (d == depth)
+		{
+			nodes.push_back(node);
+		}
+		else if (node->IsLeaf())
+		{
+			// found leaf
+			BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
+
+			// there is no further subdivision on triangle level
+			if (leaf->mTriangleBvh->GetNumNodes() == 1)
+			{
+				nodes.push_back(leaf);
+			}
+			else // more than a root => search in local bvh
+			{
+				leaf->mTriangleBvh->
+					CollectNodes(leaf->mTriangleBvh->GetRoot(), depth - d, nodes);
+			}
+		}
+		else // interior
+		{
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+
+			tStack.push(tPair(interior->mFront, d + 1));
+			tStack.push(tPair(interior->mBack, d + 1));
+		}
+	}
+}
+
+
+void Bvh::CollectAllNodes(BvhNode *node, HierarchyNodeContainer &nodes)
+{
+	BvhNodeContainer bvhNodes;
+
+	// collect nodes of geometry bvh
+	CollectNodes(mRoot, bvhNodes);
+
+	
+	// first collect nodes of object bvh
+	BvhNodeContainer::const_iterator lit, lit_end = bvhNodes.end();
+
+	for (lit = bvhNodes.begin(); lit != lit_end; ++ lit)
+	{
+		nodes.push_back(*lit);
+	}
+
+	// collect nodes of local bvh: all except roots, 
+	// because they are equivalent to geometry bvh leaves
+	for (lit = bvhNodes.begin(); lit != lit_end; ++ lit)
+	{
+		if ((*lit)->IsLeaf())
+		{
+			BvhLeaf *leaf = static_cast<BvhLeaf *>(*lit);
+
+			TriangleBvhNodeContainer hnodes;
+			leaf->mTriangleBvh->CollectNodes(leaf->mTriangleBvh->GetRoot(), hnodes);
+
+			TriangleBvhNodeContainer::const_iterator hit, hit_end = hnodes.end();
+
+			for (hit = hnodes.begin(); hit != hit_end; ++ hit)
+			{
+				TriangleBvhNode *n = *hit;
+
+				// don't include root of local bvh - it has the same bounds as the leaves
+				if (n != leaf->mTriangleBvh->GetRoot())
+					nodes.push_back(n);
+			}
+		}
+	}
+}
+
+
+void Bvh::CollectVisibleLeaves(BvhNode *node, BvhLeafContainer &leaves)
+{
+	stack<BvhNode *> tStack;
+	tStack.push(node);
+
+	while (!tStack.empty())
+	{
+		BvhNode *node = tStack.top();
+		tStack.pop();
+
+		if (IsWithinViewFrustum(node))
+		{
+			if (!node->IsLeaf())
+			{
+				BvhInterior *interior = static_cast<BvhInterior *>(node);
+
+				tStack.push(interior->mFront);
+				tStack.push(interior->mBack);
+			}
+			else
+			{
+				leaves.push_back(static_cast<BvhLeaf *>(node));
+			}
+		}
+	}
+}
+
+
+BvhLeaf *Bvh::GetRandomLeaf(BvhNode *node)
+{
+	stack<BvhNode *> nodeStack;
+  
+	nodeStack.push(node);
+
+	int mask = rand();
+
+    while (!nodeStack.empty()) 
+	{
+		BvhNode *node = nodeStack.top();
+		nodeStack.pop();
+		
+		if (node->IsLeaf()) 
+		{
+			// random leaf
+			return static_cast<BvhLeaf *>(node);
+		} 
+		else 
+		{
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+			
+			BvhNode *next;
+			
+			// random decision
+			if (mask & 1)
+				next = interior->mBack;
+			else
+				next = interior->mFront;
+
+			mask = mask >> 1;
+
+			nodeStack.push(next);
+		}
+	}
+  
+	// should never come here
+	return NULL;
+}
+
+
+BvhLeaf *Bvh::GetRandomVisibleLeaf(BvhNode *node)
+{
+	BvhLeafContainer leaves;
+	leaves.reserve(node->GetNumLeaves());
+
+	CollectVisibleLeaves(node, leaves);
+
+	if (leaves.empty())
+		return NULL;
+
+	const int r = (int)(((float)rand() / RAND_MAX) * ((float)leaves.size() - 0.5f));
+
+	return leaves[r];
+}
+
+
+void Bvh::CollectGeometry(BvhNode *node, GeometryVector &geometry)
+{
+	//geometry.reserve(node->CountPrimitives());
+
+	for (int i = node->mFirst; i <= node->mLast; ++ i)
+	{
+		geometry.push_back(mGeometry[i]);
+	}
+}
+
+
+NodeGeometry **Bvh::GetGeometry(BvhNode *node, int &geometrySize)
+{
+	geometrySize = node->CountPrimitives();
+	return mGeometry + node->mFirst;
+}
+
+
+void Bvh::UpdateInteriors(BvhNode *node)
+{
+	if (!node->IsLeaf())
+	{
+		BvhInterior *interior = static_cast<BvhInterior *> (node);
+		
+		// update the indices of the geometry so we can render interiors as well
+		UpdateInteriors(interior->GetBack());
+		UpdateInteriors(interior->GetFront());
+
+		// update area
+		interior->mFirst = min(interior->GetBack()->mFirst, interior->GetFront()->mFirst);
+		interior->mLast = max(interior->GetBack()->mLast, interior->GetFront()->mLast);
+	
+		interior->mArea = interior->GetBack()->mArea + interior->GetFront()->mArea;
+	}
+}
+
+
+int	Bvh::IsWithinViewFrustumLocal(BvhNode *node)
+{
+	bool bIntersect	= false;
+
+	if (node->GetParent())
+		node->mPlaneMask = node->GetParent()->mPlaneMask;
+
+
+	////////
+	//-- do the view frustum culling for the planes [mPreferredPlane - 5]
+
+	for (int i = node->mPreferredPlane; i < 6; ++ i)
+	{
+		//-- do the test only if necessary
+		if (node->mPlaneMask & (1 << i))
+		{
+			//-- test the n-vertex
+			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
+			{
+				//-- outside
+				node->mPreferredPlane = i;
+				return 0;
+			}
+
+			//-- test the p-vertex
+			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f)
+			{
+				//-- completely inside: children don't need to check against this plane no more
+				node->mPlaneMask^=	1 << i;
+			}
+			else 
+			{
+				bIntersect = true;
+			}
+		}
+	}
+
+	//////////
+	//-- do the view frustum culling for the planes [0 - m_iPreferredPlane)
+
+	for (int i = 0; i < node->mPreferredPlane; ++ i)
+	{
+		// do the test only if necessary
+		if (node->mPlaneMask & (1 << i))
+		{
+			//-- test the n-vertex
+			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
+			{
+				// outside
+				node->mPreferredPlane = i;
+				return 0;
+			}
+
+			//-- test the p-vertex
+			if (node->mBox.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f)
+			{
+				// completely inside: children don't need to check against this plane no more
+				node->mPlaneMask^= 1 << i;
+			}
+			else 
+			{
+				bIntersect = true;
+			}
+		}
+	}
+
+	return bIntersect ? -1 : 1;
+}
+
+
+int	Bvh::IsWithinViewFrustum(const BoundingBox &box, const int planeMask, const int preferredPlane)
+{
+	bool bIntersect	= false;
+
+	////////
+	//-- do the view frustum culling for the planes [mPreferredPlane - 5]
+
+	for (int i = preferredPlane; i < 6; ++ i)
+	{
+		//-- do the test only if necessary
+		if (planeMask & (1 << i))
+		{
+			//-- test the n-vertex
+			if (box.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
+				return 0;
+
+			//-- test the p-vertex
+			if (!(box.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f))
+				bIntersect = true;
+		}
+	}
+
+	//////////
+	//-- do the view frustum culling for the planes [0 - m_iPreferredPlane)
+
+	for (int i = 0; i < preferredPlane; ++ i)
+	{
+		// do the test only if necessary
+		if (planeMask & (1 << i))
+		{
+			//-- test the n-vertex
+			if (box.getDistance(mClipPlaneAABBVertexIndices[i][0], mFrustum.mClipPlane[i]) > 0.0f)
+				// outside
+				return 0;
+
+			//-- test the p-vertex
+			if (!(box.getDistance(mClipPlaneAABBVertexIndices[i][1], mFrustum.mClipPlane[i]) <= 0.0f))
+				bIntersect = true;
+		}
+	}
+
+	return bIntersect ? -1 : 1;
+}
+
+
+int	Bvh::IsWithinViewFrustum(BvhNode *node)
+{
+	if (node->mCache.mLastFrustumTestedFrameId == mFrameId)
+	{
+		return node->mCache.mFrustumCulled;
+	}
+
+	node->mCache.mLastFrustumTestedFrameId = mFrameId;
+
+	timeViewFrustumCulling.Entry();
+	
+	if (Settings::Global()->get_nvocc_use_tighter_bounds_for_frustum_culling())
+	{
+		const int intersect = IsWithinViewFrustumLocal(node);
+
+		if ((node->mNumTestNodes == 1) || (intersect != -1))
+		{
+			//OUT1("x");
+			node->mCache.mFrustumCulled = intersect;
+		}
+		else 
+		{
+			// maybe intersecting one of the tighter boxes
+			node->mCache.mFrustumCulled = 0;
+
+			for (int i = 0; i < node->mNumTestNodes; ++ i)
+			{
+				RenderableHierarchyNode *n = mTestNodes[node->mTestNodesIdx + i];
+
+				if (IsWithinViewFrustum(n->GetBox(), 
+										node->mPlaneMask, 
+										node->mPreferredPlane))
+				{
+					node->mCache.mFrustumCulled = -1;
+					break;
+				}
+			}
+		}
+	}
+	else
+	{
+		//OUT1("y");
+		node->mCache.mFrustumCulled = IsWithinViewFrustumLocal(node);
+	}
+
+	timeViewFrustumCulling.Exit();
+
+	return node->mCache.mFrustumCulled;
+}
+
+
+void Bvh::InitFrame(Camera *camera, const int currentFrameId, Viewer *viewer)
+{
+	// = 0011 1111 which means that at the beginning, all six planes have to frustum culled
+	mRoot->mPlaneMask = 0x3f;
+	mCamera = camera;
+
+	mFrameId = currentFrameId;
+
+	// to begin with, we must grab the plane equations of the six clipplanes of the viewfrustum
+	Matrix4f matViewing, matProjectionView;
+	Transform3D	transform;
+
+	camera->GetTransform(transform);
+	transform.get(matViewing);
+
+	camera->GetProjectionMatrix(matProjectionView);
+	matProjectionView *= matViewing;
+
+	Vec3f vec;
+	float fInvLength;
+
+
+	//////////
+	//-- extract the plane equations
+
+	for (int i = 0; i < 4; ++ i)
+	{
+		mFrustum.mClipPlane[0][i] = matProjectionView[i][3] - matProjectionView[i][0];	// right plane
+		mFrustum.mClipPlane[1][i] = matProjectionView[i][3] + matProjectionView[i][0];	// left plane
+		mFrustum.mClipPlane[2][i] = matProjectionView[i][3] + matProjectionView[i][1];	// bottom plane
+		mFrustum.mClipPlane[3][i] = matProjectionView[i][3] - matProjectionView[i][1];	// top plane
+		mFrustum.mClipPlane[4][i] = matProjectionView[i][3] - matProjectionView[i][2];	// far plane
+		mFrustum.mClipPlane[5][i] = matProjectionView[i][3] + matProjectionView[i][2];	// near plane
+	}
+
+	////////////
+	//-- normalize the coefficients and find the indices of the n- and p-vertices
+
+	for (int i = 0; i < 6; ++ i)
+	{
+		// the clipping planes look outward the frustum,
+		// so distances > 0 mean that a point is outside
+		fInvLength	= -1.0f / sqrt(	mFrustum.mClipPlane[i][0] * mFrustum.mClipPlane[i][0] +	
+									mFrustum.mClipPlane[i][1] * mFrustum.mClipPlane[i][1] +	
+									mFrustum.mClipPlane[i][2] * mFrustum.mClipPlane[i][2]);
+
+		mFrustum.mClipPlane[i][0] *= fInvLength;
+		mFrustum.mClipPlane[i][1] *= fInvLength;
+		mFrustum.mClipPlane[i][2] *= fInvLength;
+		mFrustum.mClipPlane[i][3] *= fInvLength;
+
+		vec.x = mFrustum.mClipPlane[i][0];
+		vec.y = mFrustum.mClipPlane[i][1];
+		vec.z = mFrustum.mClipPlane[i][2];
+
+		// n-vertex
+		mClipPlaneAABBVertexIndices[i][0] = BoundingBox::getIndexNearestVertex(vec);
+		// p-vertex
+		mClipPlaneAABBVertexIndices[i][1] = BoundingBox::getIndexFurthestVertex(vec);	
+	}
+
+
+	const float aspect = mCamera->GetAspect();
+	const float fov = mCamera->GetFov();
+	
+	mScale = sqrt(aspect) * 2.0f * tan(fov * mypi / 180.0f);
+
+	mNearPlane = mCamera->GetRealViewDirection();
+	mNearPlaneD = - mNearPlane * (mCamera->GetPosition() - Origin3f);
+}
+
+
+float Bvh::CalcDistance(BvhNode *node) const
+{
+	timeDistance.Entry();
+
+	float distance;
+
+#if USE_TIGHTER_BOUNDS_FOR_ALL
+	float minDist = 1e25f;
+
+	//OUT1("testing " << mNumTestNodes << " nodes");
+	for (int i = 0; i < node->mNumTestNodes; ++ i)
+	{
+		RenderableHierarchyNode *hnode = mTestNodes[node->mTestNodesIdx + i];
+
+		if(hnode->mCache.mLastDistanceTestedFrameId < mFrameId)
+		{
+			hnode->mCache.mDistance = 
+				hnode->GetBox().getMinVisibleDistance(mNearPlane, mNearPlaneD);
+			hnode->mCache.mLastDistanceTestedFrameId = mFrameId;
+		}
+
+		if (hnode->mCache.mDistance < minDist)
+			minDist = hnode->mCache.mDistance;
+	}
+
+	distance = minDist;
+
+#else
+
+	if(node->mCache.mLastDistanceTestedFrameId != mFrameId)
+	{
+		//OUT1("z " << node->mLastDistanceTestedFrameId  << " " << mFrameId);
+		node->mCache.mDistance = node->GetBox().getMinVisibleDistance(mNearPlane, mNearPlaneD);
+		node->mCache.mLastDistanceTestedFrameId = mFrameId;
+	}
+
+	distance = node->mCache.mDistance;
+#endif
+
+	timeDistance.Exit();
+
+	return distance;
+}
+
+
+float Bvh::GetSquareDistance(BvhNode *node) const
+{
+	timeDistance.Entry();
+
+	float distance;
+
+#if USE_TIGHTER_BOUNDS_FOR_ALL
+	const Vector3f nearplane = mCamera->GetRealViewDirection();
+	const float nearplaneD = - nearplane * (mCamera->GetPosition() - Origin3f);
+
+	float minDist = 1e25f;
+	
+	//OUT1("testing " << mNumTestNodes << " nodes");
+	for (int i = 0; i < node->mNumTestNodes; ++ i)
+	{
+		RenderableHierarchyNode *hnode = mTestNodes[node->mTestNodesIdx + i];
+		const float dist = GetMinSquareDistance(hnode->GetBox());
+
+		if (dist < minDist)
+			minDist = dist;
+	}
+
+	distance = minDist;
+#else
+	distance = GetMinSquareDistance(node->GetBox());
+#endif
+
+	timeDistance.Exit();
+
+	return distance;
+}
+
+
+float Bvh::GetMinSquareDistance(const BoundingBox &box) const
+{
+	float minDist = 1e25f;
+
+	const Point3f *pts = (const Point3f *)box.getVertexData();
+	const Point3f pos = mCamera->GetPosition();
+
+	for (int i = 0; i < 8; ++ i)
+	{
+		const float newDist = DistanceSquared(pts[i], pos);
+
+		if (newDist < minDist)
+			minDist = newDist;
+	}
+
+	return minDist;
+}
+
+
+bool Bvh::ExportBinTree(const string &filename)
+{
+	ofstream stream(filename.c_str(), ifstream::binary);
+	
+	if (!stream.is_open()) return false;
+
+	OUT1("exporting bvh");
+
+	std::queue<BvhNode *> tStack;
+	tStack.push(mRoot);
+
+	while(!tStack.empty())
+	{
+		BvhNode *node = tStack.front();
+		tStack.pop();
+		
+		if (node->IsLeaf())
+		{
+			//OUT1("l");
+			ExportBinLeaf(stream, static_cast<BvhLeaf *>(node));
+		}
+		else
+		{
+			//OUT1("i");
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+
+			ExportBinInterior(stream, interior);
+			
+			tStack.push(interior->mFront);
+			tStack.push(interior->mBack);
+		}
+	}
+
+	OUT1("... finished");
+
+	return true;
+}
+
+
+void Bvh::ExportBinLeaf(ofstream &stream, BvhLeaf *leaf)
+{
+	int type = TYPE_LEAF;
+	int first = leaf->mFirst;
+	int last = leaf->mLast;
+
+	stream.write(reinterpret_cast<char *>(&type), sizeof(int));
+	stream.write(reinterpret_cast<char *>(&first), sizeof(int));
+	stream.write(reinterpret_cast<char *>(&last), sizeof(int));
+}
+
+
+void Bvh::ExportBinInterior(ofstream &stream, BvhInterior *interior)
+{
+	int type = TYPE_INTERIOR;
+	stream.write(reinterpret_cast<char *>(&type), sizeof(int));
+	stream.write(reinterpret_cast<char *>(&interior->mAxis), sizeof(char));
+	stream.write(reinterpret_cast<char *>(&interior->mPosition), sizeof(float));
+}
+
+
+BvhLeaf *Bvh::ImportBinLeaf(ifstream &stream, BvhInterior *parent)
+{
+	BvhLeaf *leaf = new BvhLeaf(parent);
+
+	int first, last;
+
+	stream.read(reinterpret_cast<char *>(&first), sizeof(int));
+	stream.read(reinterpret_cast<char *>(&last), sizeof(int));
+
+	leaf->mFirst = first;
+	leaf->mLast = last;
+
+	return leaf;
+}
+
+
+BvhInterior *Bvh::ImportBinInterior(ifstream &stream, BvhInterior *parent)
+{
+	BvhInterior *interior = new BvhInterior(parent);
+
+	stream.read(reinterpret_cast<char *>(&interior->mAxis), sizeof(char));
+	stream.read(reinterpret_cast<char *>(&interior->mPosition), sizeof(float));
+
+	return interior;
+}
+
+
+BvhNode *Bvh::LoadNextNode(ifstream &stream, BvhInterior *parent)
+{
+	int nodeType;
+	stream.read(reinterpret_cast<char *>(&nodeType), sizeof(int));
+
+	if (nodeType == TYPE_LEAF)
+	{
+		//OUT1("l");
+		return ImportBinLeaf(stream, static_cast<BvhInterior *>(parent));
+	}
+
+	if (nodeType == TYPE_INTERIOR)
+	{
+		//OUT1("i");
+		return ImportBinInterior(stream, static_cast<BvhInterior *>(parent));
+	}
+
+	return NULL;
+}
+
+
+Bvh *Bvh::LoadFromFile(const string &filename, const GeometryVector &geom)
+{
+	// export binary version of mesh
+	queue<BvhNode *> tStack;
+	ifstream stream(filename.c_str(), ifstream::binary);
+
+	if (!stream.is_open()) return NULL;
+
+	OUT1("loading bvh");
+	Bvh *bvh = new Bvh();
+
+	bvh->mGeometrySize = geom.size();
+	bvh->mGeometry = new NodeGeometry*[bvh->mGeometrySize];
+
+	for (size_t i = 0; i < bvh->mGeometrySize; ++ i)
+		bvh->mGeometry[i] = geom[i];
+
+	bvh->mRoot = bvh->LoadNextNode(stream, NULL);
+
+	tStack.push(bvh->mRoot);
+	bvh->mNumNodes = 1;
+
+	while(!tStack.empty())
+	{
+		BvhNode *node = tStack.front();
+		tStack.pop();
+
+		if (!node->IsLeaf())
+		{
+			bvh->mNumNodes += 2;
+
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+
+            BvhNode *front = bvh->LoadNextNode(stream, interior);
+			BvhNode *back = bvh->LoadNextNode(stream, interior);
+	
+			interior->mFront = front;
+			interior->mBack = back;
+
+			front->mDepth = interior->mDepth + 1;
+			back->mDepth = interior->mDepth + 1;
+
+			tStack.push(front);			
+			tStack.push(back);
+		}
+	}
+
+	OUT1("... finished loading " << bvh->mNumNodes << " nodes, updating boxes");
+
+	//adjust bounding boxes
+	bvh->UpdateBoxes(bvh->mRoot);
+	bvh->UpdateNumLeaves(bvh->mRoot);
+
+	// compute unique ids
+	bvh->ComputeIds();
+	// update the indices of the geometry so we can render interiors as well
+	bvh->UpdateInteriors(bvh->mRoot);
+	// do this once so at least the current node bounding box is tested
+	bvh->RecomputeBounds();
+
+	return bvh;
+}
+
+
+float Bvh::ComputeScreenSpaceProjection(BvhNode *node) const
+{
+#if 0
+	int projection = 0;
+
+	for (int i = 0; i < node->mNumTestNodes; ++ i)
+	{
+		RenderableHierarchyNode *n = mTestNodes[node->mTestNodesIdx + i];
+
+		if(hnode->.mCache.mLastProjectionFrameId < mFrameId)
+		{
+			hnode->mCache.mProjection = 
+				ComputeScreenSpaceProjection(hnode->GetBox());
+			hnode->mCache.mLastProjectionFrameId = mCurrentFrameId;
+		}
+		
+		projection += hnode->mCache.mProjection;
+	}
+
+	return projection;
+
+#else
+
+	if(node->mCache.mLastProjectionFrameId < mFrameId)
+	{
+		node->mCache.mProjection = 
+			ComputeScreenSpaceProjection(node->GetBox());
+		node->mCache.mLastProjectionFrameId = mFrameId;
+	}
+
+	return node->mCache.mProjection;
+#endif
+}
+
+
+float Bvh::ComputeScreenSpaceProjection(const BoundingBox &box) const
+{
+#if 0
+	const float dist = CalcDistance(box);
+#else
+	//const float dist = sqrt(GetMinSquareDistance(box));
+	const float dist = Distance(mCamera->GetPosition(), box.getCenter());
+#endif
+
+	const float scale = dist ? dist * mScale : 1e-6;
+	const float f = 1.0f / scale;
+	const float coverage = f * f * box.getSurface() / 6.0f;
+
+#if 0
+	OUT1("scale: " << scale);
+	OUT1("box: " << box.getSurface());
+	OUT1("aspect: " << aspect);
+	OUT1("fov: " << fov);
+	OUT1("tan: " << tan(fov * mypi / 180.0f));
+	//OUT1("w: " << w << " h: " << h);
+	OUT1("cov: " << coverage << " pix: " << coverage * w * h << " f: " << f);
+#endif
+	return coverage;
+}
+
+
+void Bvh::RenderBoundingBox(BvhNode *node)
+{
+#if 1
+	static BvhNodeContainer dummy(1);
+	dummy[0] = node;
+	RenderBoundingBoxes(dummy);
+#else
+	RenderBoxIndividual(node);
+#endif
+}
+
+
+void Bvh::RenderBoundingBoxImmediate(const BoundingBox &box, const bool restartStrip)
+{
+	const Point3f *pU = box.getUpperPtr();
+	const Point3f *pL = box.getLowerPtr();
+
+	///////////
+	//-- render AABB as triangle strips
+
+	glVertex3f(pL->x, pL->y, pU->z);
+	glVertex3f(pU->x, pL->y, pU->z);
+	glVertex3f(pL->x, pU->y, pU->z);
+	glVertex3f(pU->x, pU->y, pU->z);
+	glVertex3f(pL->x, pU->y, pL->z); 
+	glVertex3f(pU->x, pU->y, pL->z);
+	glVertex3f(pL->x, pL->y, pL->z); 
+	glVertex3f(pU->x, pL->y, pL->z);
+
+	glPrimitiveRestartNV();
+
+	//-- render second half of AABB
+	glVertex3f(pL->x, pU->y, pU->z); 
+	glVertex3f(pL->x, pU->y, pL->z);
+	glVertex3f(pL->x, pL->y, pU->z);
+	glVertex3f(pL->x, pL->y, pL->z);
+	glVertex3f(pU->x, pL->y, pU->z); 
+	glVertex3f(pU->x, pL->y, pL->z);
+	glVertex3f(pU->x, pU->y, pU->z); 
+	glVertex3f(pU->x, pU->y, pL->z);
+
+	
+	// restart for all but last node
+	if (restartStrip)
+	{
+		glPrimitiveRestartNV();
+	}
+}
+
+
+int Bvh::RenderBoundingBoxesImmediate(const BvhNodeContainer &nodes)
+{
+	timeSetup.Entry();
+	int renderedBoxes = 0;
+	
+	glBegin(GL_TRIANGLE_STRIP);
+
+	BvhNodeContainer::const_iterator bit, bit_end = nodes.end();
+
+	for (bit = nodes.begin(); bit != bit_end; ++ bit)
+	{
+		BvhNode *node = *bit;
+
+		renderedBoxes += node->mNumTestNodes;
+#if 0
+		const bool restartStrip = false;
+		RenderBoundingBoxImmediate(node->GetBox(), restartStrip);
+#else
+		//OUT1("nodes: " << node->mTestNodes.size() << " l: " << node->IsLeaf());
+		for (int i = 0; i < node->mNumTestNodes; ++ i)
+		{
+			RenderableHierarchyNode *testNode = mTestNodes[node->mTestNodesIdx + i];
+			const BoundingBox &box = testNode->GetBox();
+
+			// restart for all but last node
+			const bool restartStrip = (bit != bit_end - 1) || (i != node->mNumTestNodes - 1);
+
+			RenderBoundingBoxImmediate(box, restartStrip);
+		}
+#endif
+	}
+	glEnd();
+	timeSetup.Exit();
+
+	return renderedBoxes;
+}
+
+
+int Bvh::RenderBoundingBoxes(const BvhNodeContainer &nodes)
+{
+	// always render in immediate mode if there is only one box to render
+	if (((nodes.size() == 1) && (nodes[0]->mNumTestNodes == 1)) || 
+		!Settings::Global()->get_nvocc_use_index_arrays())
+	{	
+		// render bounding boxes individually
+		return RenderBoundingBoxesImmediate(nodes);
+	}
+	else 
+	{
+		const int renderedBoxes = 
+			PrepareBoundingBoxesWithDrawArrays(nodes);
+		RenderBoundingBoxesWithDrawArrays(renderedBoxes);
+
+		return renderedBoxes;
+	}
+}
+
+
+int Bvh::PrepareBoundingBoxesWithDrawArrays(const BvhNodeContainer &nodes)
+{
+	//////
+	//-- for the first time we come here ...
+
+	if (!mIndices)
+	{	// create list of indices
+		CreateIndices();
+	}
+
+	if (sCurrentVboId == -1)
+	{	
+		// prepare the vbo 
+		PrepareVertices();
+	}
+
+	///////////////
+
+	timeSetup.Entry();
+
+	int numNodes = 0;
+
+	BvhNodeContainer::const_iterator nit, nit_end = nodes.end();
+
+	for (nit = nodes.begin(); nit != nit_end; ++ nit)
+	{
+		BvhNode *node = *nit;
+
+        const int numIndices = node->mNumTestNodes * sNumIndicesPerBox;
+		
+		// copy indices
+		memcpy(mIndices + numNodes * sNumIndicesPerBox, 
+			   mTestIndices + node->mIndicesPtr, 
+#if 0 //ALIGN_INDICES
+			   ((numIndices / 32) * 32 + 32) * sizeof(unsigned int));
+#else
+			   numIndices * sizeof(unsigned int));
+#endif
+		numNodes += node->mNumTestNodes;
+	}
+
+	timeSetup.Exit();
+
+	return numNodes;
+}
+
+
+void Bvh::RenderBoundingBoxesWithDrawArrays(const int numNodes)
+{
+	//////
+	//-- Rendering the vbo
+
+	timeIssueDrawElements.Entry();
+
+	if (useVbos)
+		// set the vertex pointer to the vertex buffer
+		glVertexPointer(3, GL_FLOAT, 0, (char *)NULL);		
+	else
+		glVertexPointer(3, GL_FLOAT, 0, mVertices);
+
+	// we do use the last or the first index (they are generate and only used to connect strips)
+	const int numElements = numNodes * sNumIndicesPerBox - 1;
+
+	// don't render first degenerate index
+	glDrawElements(GL_TRIANGLE_STRIP, numElements, GL_UNSIGNED_INT, mIndices + 1); 
+	
+	timeIssueDrawElements.Exit();
+}
+
+
+#define ALIGN_INDICES 1
+
+void Bvh::CreateIndices()
+{
+	// collect bvh nodes
+	BvhNodeContainer nodes;
+	CollectNodes(mRoot, nodes);
+
+	OUT1("creating new indices");
+
+	int numMaxIndices = 0;
+
+	BvhNodeContainer::const_iterator lit, lit_end = nodes.end();
+
+	for (lit = nodes.begin(); lit != lit_end; ++ lit)
+	{
+		int offset = (*lit)->mNumTestNodes * sNumIndicesPerBox;
+#if ALIGN_INDICES
+		// align with 32
+		offset = (offset / 32) * 32 + 32; 
+#endif
+		numMaxIndices += offset;
+	}
+
+
+	OUT1("creating global indices buffer");
+
+	if (mIndices) delete [] mIndices;
+	if (mTestIndices) delete [] mTestIndices;
+
+	// global buffer: create it once so we don't have
+	// to allocate memory from the chunks of the node
+	mIndices = new unsigned int[numMaxIndices];
+
+	// create new index buffer for the individual nodes
+	mTestIndices = new unsigned int[numMaxIndices];
+	
+	mCurrentIndicesPtr = 0;
+
+	for (lit = nodes.begin(); lit != lit_end; ++ lit)
+	{
+		BvhNode *node = *lit;
+		
+		// resize array
+		node->mIndicesPtr = mCurrentIndicesPtr;
+
+		int numIndices = 0;
+
+		// the bounding box of the test nodes are rendered instead of the root node
+		// => store their indices
+		for (int i = 0; i < node->mNumTestNodes; ++ i, numIndices += sNumIndicesPerBox)
+		{
+			RenderableHierarchyNode *testNode = mTestNodes[node->mTestNodesIdx + i];
+
+			// add indices to root node
+			for (int j = 0; j < sNumIndicesPerBox; ++ j)
+			{
+				mTestIndices[mCurrentIndicesPtr + numIndices + j] = sIndices[j] + testNode->GetId() * 8;
+			}
+		}
+
+		// align with 32
+#if ALIGN_INDICES
+		const int offset = (numIndices / 32) * 32 + 32;
+#else
+		const int offset = numIndices;
+#endif
+		mCurrentIndicesPtr += offset;
+	}
+}
+
+
+void Bvh::ComputeIds()
+{
+	// collect all nodes, also the nodes from local bvh
+	// warning: root nodes local bvh are not in there, as they 
+	// are equivalent geometry bvh leaves
+	HierarchyNodeContainer nodes;
+	CollectAllNodes(mRoot, nodes);
+
+	// assign ids to all nodes of the "regular" hierarchy
+	int i = 0;
+	HierarchyNodeContainer::const_iterator lit, lit_end = nodes.end();
+
+	for (lit = nodes.begin(); lit != lit_end; ++ lit, ++ i)
+	{
+		(*lit)->SetId(i);
+	}
+}
+
+
+void Bvh::PrepareVertices()
+{
+	// collect all nodes, also the nodes from local bvh
+	// warning: root nodes local bvh are not in there, as they 
+	// are equivalent geometry bvh leaves
+	HierarchyNodeContainer nodes;
+
+	nodes.reserve(GetNumNodes());
+	CollectAllNodes(mRoot, nodes);
+
+	const unsigned int bufferSize = 8 * (int)nodes.size();
+	mVertices = new Point3f[bufferSize];
+	
+	int i = 0;
+	
+	// store bounding box vertices
+	HierarchyNodeContainer::const_iterator lit, lit_end = nodes.end();
+
+	for (lit = nodes.begin(); lit != lit_end; ++ lit, i += 8)
+	{
+		RenderableHierarchyNode *node = *lit;
+		const Point3f *vertices = (const Point3f *)node->GetBox().getVertexData();
+
+		for (int j = 0; j < 8; ++ j)
+		{	
+			((Point3f *)mVertices)[node->GetId() * 8 + j] = vertices[j];
+		}
+	}
+
+	if (useVbos)
+	{
+		glGenBuffersARB(1, &sCurrentVboId);
+		glBindBufferARB(GL_ARRAY_BUFFER_ARB, sCurrentVboId);
+		glVertexPointer(3, GL_FLOAT, 0, (char *)NULL);
+	
+		glBufferDataARB(GL_ARRAY_BUFFER_ARB, 
+			            bufferSize * sizeof(Point3f), 
+			            mVertices, 
+						GL_STATIC_DRAW_ARB);
+
+		//glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);
+
+		// data handled by graphics driver from now on
+		DELAPTR(mVertices);
+
+		OUT1("***** created vbos *********");
+	}
+	else
+	{
+		sCurrentVboId = 0;
+        glVertexPointer(3, GL_FLOAT, 0, mVertices);
+
+		OUT1("created vertices");
+	}
+}
+
+
+void Bvh::SetMaxDepthForTestingChildren(const int maxDepth) 
+{
+	if (maxDepth != mMaxDepthForTestingChildren)
+	{
+		mMaxDepthForTestingChildren = maxDepth;
+		RecomputeBounds();
+	}
+}
+
+
+void Bvh::SetAreaRatioThresholdForTestingChildren(const float ratio) 
+{
+	if (ratio != mAreaRatioThreshold)
+	{
+		mAreaRatioThreshold = ratio;
+		RecomputeBounds();
+	}
+}
+
+
+void Bvh::SetVolRatioThresholdForTestingChildren(const float ratio) 
+{
+	if (ratio != mVolRatioThreshold)
+	{
+		mVolRatioThreshold = ratio;
+		RecomputeBounds();
+	}
+}
+
+
+void Bvh::SetUseTighterBoundsForTests(const bool tighterBoundsForTests)
+{
+	if (mUseTighterBoundsOnlyForLeafTests != tighterBoundsForTests)
+	{
+		mUseTighterBoundsOnlyForLeafTests = tighterBoundsForTests;
+		RecomputeBounds();
+	}
+}
+
+
+void Bvh::SetCollectTighterBoundsWithMaxLevel(bool t)
+{
+	if (mCollectTighterBoundsWithMaxLevel != t)
+	{
+		mCollectTighterBoundsWithMaxLevel = t;
+		RecomputeBounds();
+	}
+}
+
+
+void Bvh::RecomputeBounds()
+{
+	// clear old list
+	mTestNodes.clear();
+
+	// collect all nodes
+	BvhNodeContainer nodes;
+	CollectNodes(mRoot, nodes);
+
+	OUT1("recomputing bounds, children will be tested in depth " << mMaxDepthForTestingChildren);
+
+	if (mCollectTighterBoundsWithMaxLevel)
+		OUT1("creating tighter bounds using max level");
+	else
+		OUT1("creating tighter bounds using best set");
+
+	int success = 0;
+	BvhNodeContainer::const_iterator lit, lit_end = nodes.end();
+
+	for (lit = nodes.begin(); lit != lit_end; ++ lit)
+	{
+		BvhNode *node = *lit;
+
+		if (mCollectTighterBoundsWithMaxLevel)
+		{
+			//OUT1("creating tighter bounds using max level");
+
+			// recreate list of nodes that will be tested as a proxy ...
+			if (CreateNodeRenderList(node))
+				++ success;
+		}
+		else
+		{
+			//OUT1("creating tighter bounds using best set");
+
+			HierarchyNodeContainer hnodes;
+			CollectBestNodes(node, hnodes);	
+
+			// the new test nodes are added at the end of the vector
+			node->mTestNodesIdx = mTestNodes.size();
+
+			HierarchyNodeContainer::const_iterator cit;
+
+			// use the found nodes as nodes during the occlusion tests
+			for (cit = hnodes.begin(); cit != hnodes.end(); ++ cit)
+			{
+				RenderableHierarchyNode *child = *cit;
+				mTestNodes.push_back(child);
+			}
+
+			node->mNumTestNodes = (int)hnodes.size();
+		}
+		//OUT1("testnodes: " << node->mNumTestNodes);
+	}
+
+	const float p = 100.0f * (float)success / nodes.size();
+
+	OUT1("created tighter bounds for " << p << " percent of the nodes");
+
+	// recreate indices used for indirect mode rendering
+	if (mIndices)
+	{
+		CreateIndices();
+	}
+}
+
+	
+bool Bvh::CreateNodeRenderList(BvhNode *node)
+{
+	HierarchyNodeContainer children;
+
+	if (mUseTighterBoundsOnlyForLeafTests && !node->IsLeaf())
+	{
+		children.push_back(node);
+	}
+	else
+	{
+		// collect nodes that will be tested instead of the leaf node
+		// in order to get a tighter bounding box test
+		CollectNodes(node, mMaxDepthForTestingChildren, children);
+	}
+
+
+	// using the tighter bounds is not feasable in case
+	// that the tighter bounds represent nearly the same projected area
+	// as the old bounding box. Find this out using either volume or area
+	// heuristics
+
+	float vol = 0;
+	float area = 0;
+
+	HierarchyNodeContainer::const_iterator cit;
+
+	for (cit = children.begin(); cit != children.end(); ++ cit)
+	{
+		RenderableHierarchyNode *c = *cit;
+
+		area += c->GetBox().getSurface();
+		vol += c->GetBox().getVolume();
+	}
+
+	const float volRatio = vol / node->GetBox().getVolume();
+	const float areaRatio = area / node->GetBox().getSurface();
+	
+	bool success;
+
+	if ((areaRatio < mAreaRatioThreshold) &&
+		(volRatio < mVolRatioThreshold))
+	{
+		success = true;
+	}
+	else
+	{
+		// hack: only store node itself
+		children.clear();
+		children.push_back(node);
+
+		success = false;
+	}
+
+	// the new test nodes are added at the end of the vector
+	node->mTestNodesIdx = mTestNodes.size();
+
+	// use the found nodes as nodes during the occlusion tests
+	for (cit = children.begin(); cit != children.end(); ++ cit)
+	{
+		RenderableHierarchyNode *child = *cit;
+		mTestNodes.push_back(child);
+	}
+
+	node->mNumTestNodes = (int)children.size();
+
+	return success;
+}
+
+
+void Bvh::ResetNodeClassifications()
+{
+	BvhNodeContainer nodes;
+
+	nodes.reserve(GetNumNodes());
+	CollectNodes(mRoot, nodes);
+
+	BvhNodeContainer::const_iterator lit, lit_end = nodes.end();
+
+	for (lit = nodes.begin(); lit != lit_end; ++ lit)
+	{
+		(*lit)->ResetVisibility();
+	}
+}
+
+
+int Bvh::PostProcessLeaves(BvhLeafContainer &leaves)
+{
+	OUT1("creating tighter bounds for leaves");
+
+	BvhLeafContainer::const_iterator lit, lit_end = leaves.end();
+
+	int i = 0;
+	int tighter = 0;
+
+	for (lit = leaves.begin(); lit != lit_end; ++ lit, ++ i)
+	{
+		BvhLeaf *leaf = *lit;
+
+		int triangleCount;
+		// collect the triangles from the stored geometry
+		Point3f *triangles = CollectTriangles(leaf, triangleCount);
+
+		if (CreateTighterBoundsForLeaf(leaf, triangles, triangleCount))
+			++ tighter;
+
+		leaf->mArea = ComputeGeometryArea(leaf, triangles, triangleCount);
+
+		delete [] triangles;
+
+		if (i % 1000 == 999)
+			OUT1(i << " leaves processed ");
+	}
+
+	const float p = 100.0f * tighter / leaves.size();
+	OUT1("created tighter bounds for " << p << " percent of the leaves");
+
+	return tighter;
+}
+
+
+bool Bvh::CreateTighterBoundsForLeaf(BvhLeaf *leaf, Point3f *triangles, const int triangleCount)
+{
+	// create a local bvh over the triangles
+	if (leaf->mTriangleBvh) delete leaf->mTriangleBvh;
+
+	leaf->mTriangleBvh = new TriangleBvh((Triangle *)triangles, triangleCount);
+
+	const int maxDepth = Settings::Global()->get_nvocc_local_bvh_max_depth();
+	const int maxTriangles = Settings::Global()->get_nvocc_local_bvh_max_triangles();
+	const int minArea = Settings::Global()->get_nvocc_local_bvh_min_area();
+	const int splitType = Settings::Global()->get_nvocc_local_bvh_split_type();
+
+	leaf->mTriangleBvh->SetMaxDepth(maxDepth);
+	leaf->mTriangleBvh->SetMaxTriangles(maxTriangles);
+	leaf->mTriangleBvh->SetMinAreaRatio(minArea);
+	leaf->mTriangleBvh->SetSplitType(splitType);
+	
+	leaf->mTriangleBvh->SetVerboseOutput(false);
+
+	// construct
+	leaf->mTriangleBvh->Construct();
+
+	return true;
+}
+
+
+float Bvh::ComputeGeometryArea(BvhLeaf *leaf, Point3f *triangles, const int triangleCount) const
+{
+	float area = 0;
+
+	for (int i = 0; i < triangleCount; ++ i)
+	{
+		area += ((Triangle *)triangles)[i].GetArea();
+	}
+
+	return area;
+}
+
+
+Point3f *Bvh::CollectTriangles(BvhLeaf *leaf, int &triangleCount)
+{
+	// count triangles in the leaf
+	triangleCount = 0;
+
+	for (int geomIdx = leaf->mFirst; geomIdx <= leaf->mLast; ++ geomIdx)
+	{
+		Shape3D *sh = static_cast<Shape3D *>(mGeometry[geomIdx]->GetNode());
+		IndexedTriangleStripArray *strips = static_cast<IndexedTriangleStripArray *>(sh->getGeometry());
+
+		triangleCount += strips->getNumTriangles();
+	}
+
+	Triangle *triangles = new Triangle[triangleCount];
+	//OUT1("The leaf contains " << triangleCount << " triangles in " << (int)geometry.size() << " nodes");
+
+	int currentIdx = 0;
+	
+	for (int geomIdx = leaf->mFirst; geomIdx <= leaf->mLast; ++ geomIdx)
+	{
+		NodeGeometry *geom = mGeometry[geomIdx];
+
+		Shape3D *sh = static_cast<Shape3D *>(geom->GetNode());
+		
+		IndexedTriangleStripArray *strips = 
+			static_cast<IndexedTriangleStripArray *>(sh->getGeometry());
+
+		// get indices of containted triangles (strips are converted before)
+		int *indices;
+	
+		int tCount = strips->getNumTriangles();
+		int indexCount = tCount * 3;
+		
+		// compute triangles from strips
+		Point3f *coordinates = (Point3f *)strips->getCoordRef3f();
+		strips->getTriangleCoordinateIndices(0, &indices, indexCount);
+
+		for (int i = 0; i < tCount; ++ i)
+		{
+			triangles[i + currentIdx].mVertices[0] = coordinates[indices[i * 3 + 0]]; 
+			triangles[i + currentIdx].mVertices[1] = coordinates[indices[i * 3 + 1]]; 
+			triangles[i + currentIdx].mVertices[2] = coordinates[indices[i * 3 + 2]]; 
+
+			if (geom->mMatId != NV_RenderPredictorRenderAction::NO_MAT)
+			{
+				const Matrix4f &m = mRenderer->app_mats[geom->mMatId];
+				Transform3D tf(m);
+				tf.transform(triangles[i + currentIdx].mVertices[0]);
+				tf.transform(triangles[i + currentIdx].mVertices[1]);
+				tf.transform(triangles[i + currentIdx].mVertices[2]);
+			}
+		}
+
+		currentIdx += tCount;
+	}
+
+	return (Point3f *)triangles;
+}
+
+
+int Bvh::CountTriangles(BvhLeaf *leaf) const
+{
+	int triangleCount = 0;
+	
+	for (int i = leaf->mFirst; i <= leaf->mLast; ++ i)
+	{
+		triangleCount += mGeometry[i]->mNumTriangles;
+	}
+
+	return triangleCount;
+}
+
+
+void Bvh::ComputeBvhStats() 
+{
+	std::stack<BvhNode *> nodeStack;
+	nodeStack.push(mRoot);
+
+	while (!nodeStack.empty()) 
+	{
+		BvhNode *node = nodeStack.top();
+		nodeStack.pop();
+
+		if (node->IsLeaf()) 
+		{
+			BvhLeaf *leaf = static_cast<BvhLeaf *>(node);
+
+			mBvhStats.mLeafSA += leaf->mBox.getSurface();
+			mBvhStats.mLeafVol += leaf->mBox.getVolume();
+
+			TriangleBvh::TriangleBvhStats tStats;
+			leaf->mTriangleBvh->GetBvhStats(tStats);
+
+			mBvhStats.mBoundsLeafSA += tStats.mLeafSA;
+			mBvhStats.mBoundsInteriorSA += tStats.mInteriorSA;
+
+			mBvhStats.mBoundsLeafVol += tStats.mLeafVol;
+			mBvhStats.mBoundsInteriorVol += tStats.mInteriorVol;
+
+			mBvhStats.mBoundsLeavesCount += leaf->mTriangleBvh->GetNumLeaves();
+		} 
+		else 
+		{
+			mBvhStats.mInteriorSA += node->mBox.getSurface();
+			mBvhStats.mInteriorVol += node->mBox.getVolume();
+
+			BvhInterior *interior = static_cast<BvhInterior *>(node);
+			
+			nodeStack.push(interior->mBack);
+			nodeStack.push(interior->mFront);
+		}
+	}
+
+	mBvhStats.mGeometryRatio = mGeometrySize / (float)GetNumLeaves();
+	mBvhStats.mTriangleRatio = mBvhStats.mTriangles / (float)GetNumLeaves();
+}
+
+
+void Bvh::PrintBvhStats() const
+{
+	//OUT1("triangle bvh surface " << mBvhStats.mBoundsLeafSA);
+	//OUT1("triangle bvh volume " << mBvhStats.mBoundsLeafVol);
+
+	OUT1("bvh stats:");
+	OUT1("interiorNodesSA = " << mBvhStats.mInteriorSA / mRoot->mBox.getSurface());
+	OUT1("leafNodesSA = " << mBvhStats.mLeafSA / mRoot->mBox.getSurface());
+	OUT1("interiorNodesVolume = " << mBvhStats.mInteriorVol / mRoot->mBox.getVolume());
+	OUT1("leafNodesVolume = " << mBvhStats.mLeafVol / mRoot->mBox.getVolume() << "\n");
+
+	OUT1("boundsInteriorNodesSA = " << mBvhStats.mBoundsInteriorSA / mBvhStats.mLeafSA);
+	OUT1("boundsLeafNodesSA = " << mBvhStats.mBoundsLeafSA / mBvhStats.mLeafSA);
+	OUT1("boundsInteriorNodesVolume = " << mBvhStats.mBoundsInteriorVol / mBvhStats.mLeafVol);
+	OUT1("boundsLeafNodesVolume = " << mBvhStats.mBoundsLeafVol / mBvhStats.mLeafVol << "\n");
+
+	OUT1("boundsLeaves: " <<  (float)mBvhStats.mBoundsLeavesCount / (float)GetNumLeaves() << "\n");
+	OUT1("geometry per leaf: " <<  mBvhStats.mGeometryRatio);
+	OUT1("triangles per leaf: " <<  mBvhStats.mTriangleRatio);
+	OUT1("");
+}
+
+
+static void RenderBoxForViz(const BoundingBox &box)
+{
+	glBegin(GL_LINE_LOOP);
+	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
+	glEnd();
+
+	glBegin(GL_LINE_LOOP);
+	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
+	glEnd();
+
+	glBegin(GL_LINE_LOOP);
+	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
+	glEnd();
+
+	glBegin(GL_LINE_LOOP);
+	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
+	glEnd();
+
+	glBegin(GL_LINE_LOOP);
+	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getLowerPtr()->y, box.getUpperPtr()->z);
+	glEnd();
+
+	glBegin(GL_LINE_LOOP);
+	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getLowerPtr()->z);
+	glVertex3d(box.getUpperPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
+	glVertex3d(box.getLowerPtr()->x, box.getUpperPtr()->y, box.getUpperPtr()->z);
+
+	glEnd();
+}
+
+
+void Bvh::RenderBoundingBoxesForViz(const int mode)
+{
+	BvhLeafContainer leaves;
+	leaves.reserve(mRoot->GetNumLeaves());
+
+	CollectLeaves(mRoot, leaves);
+
+	BvhLeafContainer::const_iterator it, it_end = leaves.end();
+
+	for (it = leaves.begin(); it != it_end; ++ it)
+	{
+		BvhLeaf *leaf = *it;
+
+		if ((mode == 0) || (mode == 2))
+		{
+			glColor3f(1.0f, 1.0f, 1.0f);
+            RenderBoxForViz(leaf->GetBox());
+		}
+		
+		if ((mode == 1) || (mode == 2)) 
+		{
+			glColor3f(1.0f, 0, 0);
+			for (int i = 0; i < leaf->mNumTestNodes; ++ i)
+			{
+				RenderBoxForViz(mTestNodes[leaf->mTestNodesIdx + i]->GetBox());
+			}
+		}
+	}
+}
+
+
+int Bvh::CountTriangles(BvhNode *node) const
+{
+	int numTriangles = 0;
+
+	for (int i = node->mFirst; i <= node->mLast; ++ i)
+	{
+		numTriangles += mGeometry[i]->mNumTriangles;
+	}
+
+	return numTriangles;
+}
+
+
+float Bvh::GetGeometryArea(BvhNode *node) const
+{
+	return node->mArea;
+}
+
+
+float Bvh::GetAvgDepth() const
+{
+	return mAvgDepth;
+}
+
+
+static float GetNodePriority(RenderableHierarchyNode *node)
+{
+	if (node->IsLeaf())
+		return 0.0f;
+
+	float result;
+
+	if (node->GetType() == BVH_NODE)
+	{
+		BvhInterior *interior = static_cast<BvhInterior *>(node);
+		
+		// evaluate the priority of this node
+		const float correctionFactor = 0.0f;
+		//  1.0f*interior->box.FaceArea(node->axis);
+
+		result = 
+			interior->GetBox().getSurface() - 
+			(interior->GetBack()->GetBox().getSurface() + interior->GetFront()->GetBox().getSurface()) + 
+			correctionFactor;
+	}
+	else
+	{
+		TriangleBvhInterior *interior = static_cast<TriangleBvhInterior *>(node);
+		
+		// evaluate the priority of this node
+		const float correctionFactor = 0.0f;
+		//  1.0f*interior->box.FaceArea(node->axis);
+
+		result = 
+			interior->GetBox().getSurface() - 
+			(interior->GetBack()->GetBox().getSurface() + interior->GetFront()->GetBox().getSurface()) + 
+			correctionFactor;
+	}
+
+	return result;
+}
+
+
+struct NodeEntry
+{
+	NodeEntry(RenderableHierarchyNode *node, float p): mNode(node), mPriority(p)
+	{}
+
+	RenderableHierarchyNode *mNode;
+	float mPriority;
+};
+
+
+bool operator<(const NodeEntry &a, const NodeEntry &b)
+{
+	return a.mPriority < b.mPriority;
+}
+
+
+void Bvh::CollectBestNodes(RenderableHierarchyNode *node, HierarchyNodeContainer &nodes)
+{
+	int maxNodes = 32;
+	priority_queue<NodeEntry> nodeStack;
+
+	nodeStack.push(NodeEntry(node, GetNodePriority(node)));
+	nodes.clear();
+
+	float SA = node->GetBox().getSurface();
+	float saThreshold = 1.1f * SA;
+	int numNodes = 1;
+	OUT1(numNodes << " " << SA);
+
+	while (!nodeStack.empty() &&
+		int(nodes.size() + nodeStack.size()) < maxNodes)
+	{
+		NodeEntry entry = nodeStack.top();
+		nodeStack.pop();
+
+		RenderableHierarchyNode *current = entry.mNode;
+
+		if (current->IsLeaf())
+		{
+			// check if there further subdivision on triangle level
+			if ((current->GetType() == BVH_NODE) && (((BvhLeaf *)current)->mTriangleBvh->GetNumNodes() > 1))
+			{
+				// push back the root of the local bvh
+				nodeStack.push(NodeEntry(((BvhLeaf *)current)->mTriangleBvh->GetRoot(), entry.mPriority));
+			}
+			else // terminate traversal
+			{
+				nodes.push_back(current);
+			}
+		}
+		else // interior node
+		{
+			// surface area of child nodes
+			SA -= entry.mPriority;
+		
+			// finish the search if we have too much fillrate increase
+			if (SA > saThreshold)
+			{
+				OUT1("break at " << SA << " > " << saThreshold);
+				nodes.push_back(current);
+			}
+			else
+			{
+				++ numNodes;
+
+				OUT1(numNodes << " " << SA << " " << entry.mPriority);
+
+				if (current->GetType() == BVH_NODE)
+				{
+					BvhInterior *interior = static_cast<BvhInterior *>(current);
+
+					nodeStack.push(NodeEntry(interior->GetBack(), 
+						GetNodePriority(interior->GetBack())));
+					nodeStack.push(NodeEntry(interior->GetFront(),
+						GetNodePriority(interior->GetFront())));
+				}
+				else
+				{
+					TriangleBvhInterior *interior = static_cast<TriangleBvhInterior *>(current);
+
+					nodeStack.push(NodeEntry(interior->GetBack(), 
+						GetNodePriority(interior->GetBack())));
+					nodeStack.push(NodeEntry(interior->GetFront(),
+						GetNodePriority(interior->GetFront())));
+				}
+			}
+		}
+	}
+
+	while (!nodeStack.empty()) 
+	{
+		NodeEntry entry = nodeStack.top();
+	
+		nodeStack.pop();
+		RenderableHierarchyNode *current = entry.mNode;
+		nodes.push_back(current);
+	}
+} 
+
+#endif
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/RenderTraverser.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/RenderTraverser.h	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/RenderTraverser.h	(revision 2746)
@@ -1,8 +1,6 @@
-#ifndef NO_PRAGMA_ONCE
-#pragma once
-#endif
-
 #ifndef RENDERTRAVERSER_H
 #define RENDERTRAVERSER_H
+
+#if TOIMPLEMENT
 
 #include "glInterface.h"
@@ -12,5 +10,5 @@
 #include <stack>
 
-using namespace std;
+
 
 typedef stack<HierarchyNode *> TraversalStack;
@@ -133,4 +131,4 @@
 	bool mUseOptimization;
 };
-
+#endif
 #endif // RENDERTRAVERSER_H
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/HierarchyNode.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/HierarchyNode.cpp	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/HierarchyNode.cpp	(revision 2746)
@@ -1,2 +1,4 @@
+#if TOIMPLEMENT
+
 #include "HierarchyNode.h"
 #include "glInterface.h"
@@ -486,2 +488,4 @@
 	mDistance = distance;
 }
+
+#endif
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.h	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Bvh.h	(revision 2746)
@@ -0,0 +1,897 @@
+// This file has been written by Jiri Bittner, October 2006
+
+#ifndef __BVH_H
+#define __BVH_H
+
+#if TOIMPLEMENT
+
+#include "Bounds.h"
+#include "PerfTimer.h"
+#include "NV_RenderPredictorRenderAction.h"
+#include "VisibilityPredictor.h"
+#include "HierarchyNode.h"
+#include "FlexibleHeap.h"
+
+
+////////
+// Forward declarations
+
+class Scene;
+class Node;
+class Camera;
+class Box;
+class BvhLeaf;
+class BvhNode;
+class TriangleBvh;
+
+
+
+/** A node in the bv hierarchy.
+*/
+class BvhNode: public RenderableHierarchyNode, public Heapable
+{
+	friend class Bvh;
+
+public:
+
+	/** visibility related options
+	*/
+	struct VisibilityInfo
+	{
+		VisibilityInfo() { Reset(); }		
+		/** Reset the visibility info.
+		*/
+		void Reset();
+
+
+		/// if the node is visible
+		bool mIsVisible;
+		/// if all the leaves are visible
+		bool mIsFullyVisible;
+		
+		/// #pixels this node was visible for the last test
+		int mVisiblePixels;		
+		/// #frames this node is assumed to stay visible
+		int mAssumedVisibleFrames;
+		/// #frames this node is assumed to stay visible
+		int mRemainingVisibleFrames;
+		
+		/// the frame in which the node was last tested visible
+		int mLastTestedVisibleFrame;
+		/// the frame when this node was last touched during traversal
+		int mLastVisitedFrame;
+		/// frame id when the node turned visible
+		int mTurnedVisibleFrame;
+		/// when was the node last tested
+		int mLastTestedFrame;
+
+		/// #times this node was invisible (only valid if the node actually is invisible)
+		int mTimesInvisible;
+		/// #times this node was tested
+		int mTimesTested;
+		/// #times this node was tested
+		int mTimesChangedClassification;
+		/// The ratio of classification changes / tests
+		float mAvgChangedClassification;
+		
+		/// if the node is view frustum culled
+		bool mIsFrustumCulled;
+
+		bool mIsNew;
+	};
+
+	/** Constructor taking the parent node.
+	*/
+	BvhNode(BvhNode *parent);
+	
+	/** Returns true if this node is a leaf.
+	*/
+	//virtual bool IsPseudoLeaf() = 0;
+	/** Virtual destructor doing nothing.
+	*/
+	virtual ~BvhNode();
+	/** Returns unique id for this node.
+	*/
+	inline int GetId() {return mId;}
+	/** Depth of this node in the tree.
+	*/
+	inline int GetDepth() const {return (int)mDepth;}
+	/** Returns parent of this bvh node, NULL if it is root.
+	*/
+	inline BvhNode *GetParent() {return mParent;}
+	/** Number of leaves in the subtree induced by this node.
+	*/
+	inline int GetNumLeaves() const {return mNumLeaves;}
+	/** Box used for visualization.
+	*/
+	Box *GetOrCreateVizBox();
+	/** Reset the status of this node.
+	*/
+	virtual void ResetVisibility();
+
+	virtual int GetType() { return BVH_NODE; }
+
+
+	/////////////////////
+	//-- visibility culling related functions
+
+	inline int GetLastVisitedFrame() const;
+
+	inline void SetLastVisitedFrame(const int lastVisited);
+	/** If this node is considered visible.
+	*/
+	inline bool IsVisible() const;
+	/** Set visibility flag of the node.
+	*/
+	inline void SetVisible(const bool visible);
+	/** The assumed visible time span of this node.
+	*/
+	inline void SetAssumedVisibleFrames(const int t);
+	/** See set.
+	*/
+	inline int GetAssumedVisibleFrames() const;
+	/** The time span this node remains visible.
+	*/
+	inline void SetRemainingVisibleFrames(const int t);
+	/** See set.
+	*/
+	inline int GetRemainingVisibleFrames() const;
+	/** Decrease the time this node is considered visible.
+	*/
+	inline void DecRemainingVisibleFrames();
+	/** Returns the frame id when this node was tested
+		visible.
+	*/
+	inline void SetLastTestedVisibleFrame(const int t);
+	/** See set.
+	*/
+	inline int GetLastTestedVisibleFrame() const;
+	/** Increases the #times this node was
+		successfully tested invisible.
+	*/
+	inline void IncTimesInvisible();
+	/** If the subtree induced by this node
+		is fully visible.
+	*/
+	inline bool IsFullyVisible();
+
+	inline int GetTimesInvisible() const;
+	inline void SetTimesInvisible(const int t);
+	
+	inline int GetTurnedVisibleFrame() const;
+	inline void SetTurnedVisibleFrame(const int turnedVisibleFrame);
+
+	inline int GetLastTestedFrame();
+	inline void SetLastTestedFrame(const int lastTested);
+
+	inline void SetVisiblePixels(const int pixels);
+	inline int GetVisiblePixels();
+
+	inline void IncTimesTested();
+	inline void IncTimesChangedClassficiation();
+
+	inline int GetTimesTested();
+	inline int GetTimesChangedClassification();
+	
+	inline bool IsViewFrustumCulled() const;
+	inline void SetViewFrustumCulled(const bool frustumCulled);
+
+	inline bool IsNew() const;
+	inline void SetIsNew(const bool isNew);
+
+
+	//////////
+	//-- rendering
+	
+	/** Returns the frame in which this node was last rendered.
+	*/
+	inline int GetLastRenderedFrame() const;
+	/** See get.
+	*/
+	inline void SetLastRenderedFrame(int lastRenderedFrame);
+	
+	
+	/** Does this node contain renderable geometry?
+	*/
+	inline bool Empty() const;
+	/** Counts #geometry stored in the subtree.
+	*/
+	inline int CountPrimitives() const;
+
+
+	///////////////
+	//-- public stuff
+
+	/// Visibility measurements for the period this node was visible
+	VisibilityMeasurementBuffer mMeasurements;
+	/// some flags
+	int mFlags;
+
+	//int mNumFailedTests;
+
+protected:
+
+	/** Marks the subtree induced by this node
+		as fully visible.
+	*/
+	inline void SetFullyVisible(const bool visible);
+
+
+	/////////////
+
+	/// the depth of this node
+	unsigned char mDepth;
+	/// the split axis
+	char mAxis;
+	/// the parent node
+	BvhNode *mParent;
+	
+	/// stores the visibility related parameters
+	VisibilityInfo mVisibility;
+
+	/////////
+	// used for view frustum culling
+
+	int mPlaneMask;
+	int mPreferredPlane;
+
+
+	////////////
+	//-- rendering related options
+	
+	/// when the node was last rendered
+	int mLastRenderedFrame;
+	/// #leaves under this node
+	int mNumLeaves;
+	
+	// Indices to first and last triangle in the triangle array
+	// assumes the traingle are placed in continuous chunk of memory
+	// however this need not be a global array!
+	
+	/// the index of the first triangle
+	int mFirst;
+	/// the index of the last triangle
+	int mLast;
+
+	/// these nodes can be tested instead of the current node
+	int mTestNodesIdx;
+	int mNumTestNodes;
+	int mIndicesPtr;
+
+	/// Area of this node
+	float mArea;
+};
+
+
+/////////////////
+//-- public inline functions
+
+int BvhNode::GetLastVisitedFrame() const 
+{ 
+	return mVisibility.mLastVisitedFrame; 
+}
+
+
+void BvhNode::SetLastVisitedFrame(const int lastVisited) 
+{ 
+	mVisibility.mLastVisitedFrame = lastVisited; 
+}
+
+
+bool BvhNode::IsVisible() const 
+{ 
+	return mVisibility.mIsVisible; 
+}
+
+
+void BvhNode::SetVisible(const bool visible) 
+{ 
+	mVisibility.mIsVisible = visible; 
+}
+
+
+void BvhNode::SetLastTestedVisibleFrame(const int t) 
+{ 
+	mVisibility.mLastTestedVisibleFrame = t; 
+}
+
+
+int BvhNode::GetLastTestedVisibleFrame() const 
+{ 
+	return mVisibility.mLastTestedVisibleFrame; 
+}
+
+
+void BvhNode::IncTimesInvisible() 
+{ 
+	++ mVisibility.mTimesInvisible; 
+}
+
+
+bool BvhNode::IsFullyVisible() 
+{ 
+	return mVisibility.mIsFullyVisible; 
+}
+
+
+int BvhNode::GetTimesInvisible() const 
+{ 
+	return mVisibility.mTimesInvisible; 
+}
+
+
+void BvhNode::SetTimesInvisible(const int t) 
+{ 
+	mVisibility.mTimesInvisible = t; 
+}
+
+
+int BvhNode::GetTurnedVisibleFrame() const 
+{ 
+	return mVisibility.mTurnedVisibleFrame; 
+}
+
+
+void BvhNode::SetTurnedVisibleFrame(const int turnedVisibleFrame) 
+{ 
+	mVisibility.mTurnedVisibleFrame = turnedVisibleFrame; 
+}
+
+
+int BvhNode::GetLastTestedFrame() 
+{ 
+	return mVisibility.mLastTestedFrame; 
+}
+
+
+void BvhNode::SetLastTestedFrame(const int lastTested)
+{ 
+	mVisibility.mLastTestedFrame = lastTested; 
+}
+
+
+void BvhNode::SetVisiblePixels(const int pixels)
+{ 
+	mVisibility.mVisiblePixels = pixels; 
+}
+
+
+int BvhNode::GetVisiblePixels() 
+{
+	return mVisibility.mVisiblePixels;
+}
+
+
+void BvhNode::IncTimesTested() 
+{ 
+	++ mVisibility.mTimesTested; 
+}
+
+
+void BvhNode::IncTimesChangedClassficiation() 
+{ 
+	++ mVisibility.mTimesChangedClassification; 
+}
+
+
+int BvhNode::GetTimesTested() 
+{ 
+	return mVisibility.mTimesTested; 
+}
+
+
+int BvhNode::GetTimesChangedClassification() 
+{ 
+	return mVisibility.mTimesChangedClassification; 
+}
+
+
+bool BvhNode::IsViewFrustumCulled() const 
+{ 
+	return mVisibility.mIsFrustumCulled; 
+}
+
+
+void BvhNode::SetViewFrustumCulled(const bool frustumCulled) 
+{ 
+	mVisibility.mIsFrustumCulled = frustumCulled; 
+}
+
+
+bool BvhNode::IsNew() const 
+{ 
+	return mVisibility.mIsNew;
+}
+
+
+void BvhNode::SetIsNew(const bool isNew) 
+{ 
+	mVisibility.mIsNew = isNew; 
+}
+
+
+int BvhNode::GetLastRenderedFrame() const 
+{ 
+	return mLastRenderedFrame; 
+}
+	
+
+void BvhNode::SetLastRenderedFrame(int lastRenderedFrame)
+{
+	mLastRenderedFrame = lastRenderedFrame;
+}
+	
+	
+bool BvhNode::Empty() const 
+{ 
+	return mFirst == -1; 
+}
+
+
+int BvhNode::CountPrimitives() const 
+{
+	return mLast - mFirst + 1; 
+}
+
+
+void BvhNode::SetAssumedVisibleFrames(const int t) 
+{ 
+	mVisibility.mAssumedVisibleFrames = t; 
+}
+
+
+int BvhNode::GetAssumedVisibleFrames() const 
+{ 
+	return mVisibility.mAssumedVisibleFrames;
+}
+
+
+void BvhNode::SetRemainingVisibleFrames(const int t) 
+{ 
+	mVisibility.mRemainingVisibleFrames = t; 
+}
+
+
+int BvhNode::GetRemainingVisibleFrames() const 
+{ 
+	return mVisibility.mRemainingVisibleFrames;
+}
+
+
+void BvhNode::DecRemainingVisibleFrames() 
+{
+	-- mVisibility.mRemainingVisibleFrames; 
+}
+
+
+
+/** Internal node of the bv hierarchy.
+*/
+class BvhInterior: public BvhNode
+{
+friend class Bvh;
+
+public:
+	BvhInterior(BvhNode *parent): mBack(NULL), mFront(NULL), BvhNode(parent)
+	{}
+	virtual bool IsLeaf() {return false;}
+	/** Back child.
+	*/
+	inline BvhNode *GetBack() { return mBack;}
+	/** Front child.
+	*/
+	inline BvhNode *GetFront() { return mFront;}
+	/** recursivly delete tree.
+	*/
+	virtual ~BvhInterior() { if (mBack) delete mBack; if (mFront) delete mFront;}
+	/** Returns split axis of this interior node.
+	*/
+	inline int GetAxis() {return (int)mAxis;}
+	/** Returns position of the split axis.
+	*/
+	inline float GetPosition() {return (float)mPosition;}
+
+
+protected:
+
+	/// the position of the split plane
+	float mPosition;
+	
+	BvhNode *mBack;
+	BvhNode *mFront;
+};
+
+
+class BvhLeaf: public BvhNode
+{
+friend class Bvh;
+
+public:
+
+	BvhLeaf(BvhNode *parent): BvhNode(parent), mTriangleBvh(NULL)
+	{}
+
+	~BvhLeaf();
+
+	virtual bool IsLeaf() { return true; }
+};
+
+
+struct NodeGeometry
+{
+	virtual AxisAlignedBox3 GetBoundingBox() = 0;
+};
+
+
+struct ObjectGeometry
+{
+	virtual AxisAlignedBox3 GetBoundingBox() = 0;
+};
+
+
+struct TriangleGeometry
+{
+	virtual AxisAlignedBox3 GetBoundingBox() = 0;
+};
+
+
+/** Class representing a bounding volume hierarchy.
+*/
+class Bvh
+{
+	/** Bvh properties
+	*/
+	struct BvhStats
+	{
+		BvhStats(): 
+		mInteriorSA(0),
+		mLeafSA(0),
+		mInteriorVol(0),
+		mLeafVol(0),
+		mBoundsInteriorSA(0),
+		mBoundsLeafSA(0),
+		mBoundsInteriorVol(0),
+		mBoundsLeafVol(0),
+		mBoundsLeavesCount(0),
+		mTriangles(0),
+		mTriangleRatio(0),
+		mGeometryRatio(0),
+		mMaxGeometry(0),
+		mMaxTriangles(0)
+		{}
+		
+		float mInteriorSA;
+		float mLeafSA;
+		float mInteriorVol;
+		float mLeafVol; 
+		float mBoundsInteriorSA;
+		float mBoundsLeafSA;
+		float mBoundsInteriorVol;
+		float mBoundsLeafVol;
+		int mBoundsLeavesCount;
+		int mTriangles;
+
+		float mTriangleRatio;
+		float mGeometryRatio;
+
+		int mMaxGeometry;
+		int mMaxTriangles;
+	};
+
+public:
+
+	/** Returns number of bvh nodes.
+	*/
+	inline int GetNumNodes() const {return mNumNodes;}
+	/** Returns number of bvh leaves.
+	*/
+	inline int GetNumLeaves() const {return mNumNodes / 2 + 1;}
+	/** Constructs the bounding volume hierarchy.
+	*/
+	void Construct();
+	/** Returns root node of the bvh.
+	*/
+	BvhNode *GetRoot() {return mRoot;}
+	/** Counts the triangle in this leaf.
+	*/
+	int CountTriangles(BvhLeaf *leaf) const;
+
+
+	//////////////////////
+
+	/** Returns geometry by reference (faster).
+	*/
+	NodeGeometry **GetGeometry(BvhNode *node, int &size);
+
+
+	/////////////
+	//-- Rendering related options
+
+	/** Renders the bounding box of this node (Or the tigher bounding boxes of some subnodes).
+	*/
+	void RenderBoundingBox(BvhNode *node);
+	/** Renders bounding boxes of the collection of nodes.
+		@returns #rendered boxes
+	*/
+	int RenderBoundingBoxes(const BvhNodeContainer &nodes);
+
+
+
+	//////////////
+	//-- Traversal related options
+
+	/** Pulls up visible classification.
+	*/
+	void MakeParentsVisible(BvhNode *node);
+	/** Does the view frustum culling for this node with respect to the previous culls
+		@returns: 0 if completely outside, 1 if completely inside, -1 if intersecting (partly inside),
+	*/
+	int	IsWithinViewFrustum(BvhNode *node);
+	/** sets frame dependent values.
+	*/
+	void InitFrame(Camera *camera, const int currentFrameId, Viewer *viewer);
+	/** this gives the orthogonal distance from the viewpoint to the nearest BBox-Vertex
+		note that negative values can appear because culling is done only afterwards
+	*/
+	float CalcDistance(BvhNode *node) const;
+	/** Sets maximal depth for taking the bounding boxes to test the
+		visibility of a node. The deeper the more the box fits to the geometry.
+	*/
+	void SetMaxDepthForTestingChildren(const int maxDepth);
+	/** Resize the visibility buffers of the bvh nodes.
+	*/
+	void ResizeVisibilityBuffers(const int size);
+
+	
+	/** Pulls up the fully visible classification in the bvh.
+	*/
+	void UpdateFullVisibility(BvhNode *node) const;
+	/** Pulls up the last visited classification in the bvh.
+	*/
+	void PullUpLastVisited(BvhNode *node, const int frameId) const;
+	/** Resets the node classifications in the tree.
+	*/
+	void ResetNodeClassifications();
+	/** Helper function that renders the bounding boxes of the leaves as
+		wireframes for visualization purpose.
+	*/
+	void RenderBoundingBoxesForViz(const int mode);
+	/** Returns squared min distance to the view point.
+	*/
+	float GetSquareDistance(BvhNode *node) const;
+	/** Count triangles the node contains.
+	*/
+	int CountTriangles(BvhNode *node) const;
+	/** Returns area of the geometry contained in the node.
+	*/
+	float GetGeometryArea(BvhNode *node) const;
+
+
+	////////////
+	//-- functions that change the boundaries of the nodes
+
+	void SetUseTighterBoundsForTests(bool tighterBoundsForTests);
+
+	void SetAreaRatioThresholdForTestingChildren(const float ratio);
+
+	void SetVolRatioThresholdForTestingChildren(const float ratio);
+	/** Returns depth of bvh hierarchy.
+	*/
+	float GetAvgDepth() const;
+
+
+	const BvhStats &GetBvhStats() const {return mBvhStats;}
+
+	void SetCollectTighterBoundsWithMaxLevel(bool t);
+
+
+	/////////////
+	//-- timers
+
+	mutable PerfTimer timeSetup;
+	mutable PerfTimer timeIssueDrawElements;
+	mutable PerfTimer timeViewFrustumCulling;
+	mutable PerfTimer timeDistance;
+
+	//////////////
+
+	static unsigned int sCurrentVboId;
+	
+
+protected:
+
+	/** Small struct representing a frustum.
+	*/
+	struct Frustum
+	{
+		/// the 6 clip planes
+		float mClipPlane[6][4];
+	};
+	
+
+
+	////////////////////
+
+	/** Constructor taking the geometry and a pointer to the current render action.
+	*/
+	Bvh(const GeometryVector &geometry, 
+		DistanceSortedRenderAction *const renderer);
+	/** protected constructor: do nothing.
+	*/
+	Bvh(): mCamera(NULL), mFrameId(-1), mVertices(NULL), mRenderer(NULL) {}
+	/** Destructor.
+	*/
+	~Bvh();
+
+
+
+	//-- sorting functions
+
+	/** The method of subdividing objects into halves in some axis using spatial median
+	*/
+	int	SortTrianglesSpatialMedian(BvhLeaf *leaf, const int axis);
+	/** The method of subdividing objects into halves in some axis
+		using object median.
+	*/
+	int	SortTrianglesObjectMedian(BvhLeaf *leaf, const int axis, float &pos);
+	/** sort triangles along the axis with respect to the splitting plane 
+		given by axis/position return index of the first tiangles belong 
+		to the front of the splitting plane.
+	*/
+	int	SortTriangles(BvhLeaf *leaf, const int axis, const float position);
+	/** sort triangles by their area.
+	*/
+	int SortTrianglesSurfaceArea(BvhLeaf *leaf, const float sa);
+
+
+	/////////////
+
+
+	/** Subdivides the current leaf.
+	*/
+	BvhNode *SubdivideLeaf(BvhLeaf *leaf, const int parentAxis);	
+	/** select splitting plane using SAH - returns the position of the splitting plane.
+	*/
+	float SelectPlaneSah(BvhLeaf *leaf, int &axis, float &minCost);
+	/** evaluate the SAH cost of a particulkar splitting plane
+	*/
+	float EvaluateSahCost(BvhLeaf *leaf, const int axis, const float position);
+
+	void ComputeBvhStats();
+	void PrintBvhStats() const;
+
+	/** Update children nodes recursively.
+	*/
+	void UpdateBoxes(BvhNode *node);
+	/** Update box for a leaf node
+	*/
+	void UpdateLeafBox(BvhLeaf *leaf);
+	/** Returns true if termination criteria met.
+	*/
+	bool TerminationCriteriaMet(BvhLeaf *leaf) const;
+	/** Compute unique ids for the bvh nodes.
+	*/
+	void ComputeIds();
+	/** Helper method that updates the number of leaves in the subtree under
+		this node.
+	*/
+	void UpdateNumLeaves(BvhNode *node) const;
+
+
+	//////////
+	//-- Helper methods for bounding box rendering in immediate and vbo mode.
+	
+	void PrepareVertices();
+
+	int PrepareBoundingBoxesWithDrawArrays(const BvhNodeContainer &nodes);
+	void RenderBoundingBoxesWithDrawArrays(const int numNodes);
+
+	int RenderBoundingBoxesImmediate(const BvhNodeContainer &nodes);
+	/** Renders a bounding box in immediate mode using index restart
+		and restarts the strip only if wished.
+	*/
+	void RenderBoundingBoxImmediate(const BoundingBox &box, const bool restartStrip);
+	/** Create the indices that each node needs to use vbo rendering.
+	*/
+	void CreateIndices();
+	/** Create the list of nodes whose bounding boxes are tested instead of the
+		bounding box of the node itself.
+	*/
+	bool CreateNodeRenderList(BvhNode *node);
+	/** Recursivly updates indices so we can 
+		render also interior nodes without traversing hierarchy
+	*/
+	void UpdateInteriors(BvhNode *node);
+	/** Recomputes the boundaries of the nodes. This function is always called if
+		some boundary options are changed.
+	*/
+	void RecomputeBounds();
+	/** Does some postprocessing on the leaves.
+		@returns #leaves that were chosen for tighter bounds.
+	*/
+	int PostProcessLeaves(BvhLeafContainer &leaves);
+
+	bool CreateTighterBoundsForLeaf(BvhLeaf *leaf, Point3f *triangles, const int triangleCount);
+
+	Point3f *CollectTriangles(BvhLeaf *leaf, int &triangleCount);
+
+	int	IsWithinViewFrustumLocal(BvhNode *node);
+
+	int	IsWithinViewFrustum(const BoundingBox &box, const int planeMask, const int preferredPlane);
+	/** Compute screen space projection of a bounding box.
+	*/
+	float ComputeScreenSpaceProjection(const BoundingBox &box) const;
+
+	float GetMinSquareDistance(const BoundingBox &box) const;
+	/** Computes the area of the triangles in a leaf.
+	*/
+	float ComputeGeometryArea(BvhLeaf *leaf, Point3f *triangles, const int triangleCount) const;
+
+
+	
+	////////////////////////
+
+	/// the root of the hierarchy
+	BvhNode *mRoot;
+	/// pointers to the geometry associated with this node
+	NodeGeometry **mGeometry;
+	/// #of entities of geometry
+	size_t mGeometrySize;
+
+
+	/////////
+	//-- termination criteria
+
+	int mMaxGeometry;
+	int mMaxTriangles;
+	int mSplitType;
+	int mNumNodes;
+	int mMaxDepth;
+	float mMinArea;
+
+	////////////////
+
+
+	/// these values are valid for all nodes
+	char mClipPlaneAABBVertexIndices[6][2];
+	/// the current view frustum
+	Frustum mFrustum;
+	/// the current camera
+	Camera *mCamera;
+	// current frame id
+	int mFrameId;
+	/// a vertex array used if working with indexed arrays (without vbo)
+	Point3f *mVertices;
+	/// indices used for draw array rendering
+	unsigned int *mIndices;
+
+	/** maximal depth from which children are fetched for 
+		testing instead of the current node
+	*/
+	int mMaxDepthForTestingChildren;
+
+	float mAreaRatioThreshold;
+	float mVolRatioThreshold;
+
+	/// pointer to the renderaction
+	DistanceSortedRenderAction *const mRenderer;
+
+	BvhStats mBvhStats;
+
+	HierarchyNodeContainer mTestNodes;
+	
+	unsigned int *mTestIndices;
+
+	/// a pointer to the end of the indices array
+	int mCurrentIndicesPtr;
+
+	float mScale;
+
+	float mAvgDepth;
+
+	Vector3f mNearPlane;
+	float mNearPlaneD;
+};
+
+#endif
+#endif // __BVH_H
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.cpp	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.cpp	(revision 2746)
@@ -0,0 +1,287 @@
+#include "Matrix4x4.h"
+#include "Vector3.h"
+#include <float.h>
+
+
+
+//namespace CHCDemo
+//{
+
+// Given min a vector to minimize and a candidate vector, replace
+// elements of min whose corresponding elements in Candidate are
+// smaller.  This function is used for finding objects' bounds,
+// among other things.
+void Minimize(Vector3 &min, const Vector3 &Candidate)
+{
+	if (Candidate.x < min.x)
+		min.x = Candidate.x;
+	if (Candidate.y < min.y)
+		min.y = Candidate.y;
+	if (Candidate.z < min.z)
+		min.z = Candidate.z;
+}
+
+// Given min a vector to minimize and a candidate vector, replace
+// elements of min whose corresponding elements in Candidate are
+// larger.  This function is used for finding objects' bounds,
+// among other things.
+void Maximize(Vector3 &max, const Vector3 &Candidate)
+{
+	if (Candidate.x > max.x)
+		max.x = Candidate.x;
+	if (Candidate.y > max.y)
+		max.y = Candidate.y;
+	if (Candidate.z > max.z)
+		max.z = Candidate.z;
+}
+
+// Project the vector onto the YZ, XZ, or XY plane depending on which.
+// which      Coordinate plane to project onto
+// 0          YZ
+// 1          XZ
+// 2          XY
+// This function is used by the polygon intersection code.
+
+void Vector3::ExtractVerts(float *px, float *py, int which) const
+{
+	switch (which) {
+	case 0:
+		*px = y;
+		*py = z;
+		break;
+	case 1:
+		*px = x;
+		*py = z;
+		break;
+	case 2:
+		*px = x;
+		*py = y;
+		break;
+	}
+}
+
+// returns the axis, where the vector has the largest value
+int
+Vector3::DrivingAxis(void) const
+{
+	int axis = 0;
+	float val = fabs(x);
+
+	if (fabs(y) > val) {
+		val = fabs(y);
+		axis = 1;
+	}
+
+	if (fabs(z) > val)
+		axis = 2;
+
+	return axis;
+}
+
+// returns the axis, where the vector has the smallest value
+int
+Vector3::TinyAxis(void) const
+{
+	int axis = 0;
+	float val = fabs(x);
+
+	if (fabs(y) < val) {
+		val = fabs(y);
+		axis = 1;
+	}
+
+	if (fabs(z) < val)
+		axis = 2;
+
+	return axis;
+}
+
+
+// Construct a view vector ViewN, and the vector ViewU perpendicular
+// to ViewN and lying in the plane given by ViewNxUpl
+// the last vector of ortogonal system is ViewV, that is
+// perpendicular to both ViewN and ViewU. 
+// |ViewN| = |ViewU| = |ViewV| = 1
+// The ViewN vector pierces the center of the synthesized image
+//     ViewU vector goes from the center image rightwards
+//     ViewV vector goes from the center image upwards
+void
+ViewVectors(const Vector3 &DirAt, const Vector3 &Viewer,
+			const Vector3 &UpL, Vector3 &ViewV, Vector3 &ViewU,
+			Vector3 &ViewN)
+{
+	Vector3 U, V, N;
+	Vector3 Up = Normalize(UpL);
+
+	N = -Normalize(DirAt);
+
+	V = Normalize(Up - DirAt);
+	V -= N * DotProd(V, N);
+	V = Normalize(V);
+	U = CrossProd(V, N);
+
+	ViewU = U; // rightwards
+	ViewV = V; // upwards
+	ViewN = -N; // forwards
+#ifdef GTP_DEBUG
+	const float eps = 1e-3f;
+	if (fabs(Magnitude(ViewU) - 1.0) > eps) {
+		Debug << "ViewU magnitude error= " << Magnitude(ViewU) << "\n";
+	}
+	if (fabs(Magnitude(ViewV) - 1.0) > eps) {
+		Debug << "ViewU magnitude error= " << Magnitude(ViewV) << "\n";
+	}
+	if (fabs(Magnitude(ViewN) - 1.0) > eps) {
+		Debug << "ViewU magnitude error= " << Magnitude(ViewN) << "\n";
+	}
+#endif // GTP_DEBUG
+
+	return;
+}
+
+// Given the intersection point `P', you have available normal `N'
+// of unit length. Let us suppose the incoming ray has direction `D'.
+// Then we can construct such two vectors `U' and `V' that
+// `U',`N', and `D' are coplanar, and `V' is perpendicular
+// to the vectors `N','D', and `V'. Then 'N', 'U', and 'V' create
+// the orthonormal base in space R3.
+void
+TangentVectors(Vector3 &U,
+			   Vector3 &V, // output
+			   const Vector3 &normal, // input
+			   const Vector3 &dirIncoming)
+{
+#ifdef GTP_DEBUG
+	float d = Magnitude(normal); 
+	if ( (d < 0.99) ||
+		(d > 1.01) ) {
+			Debug << " The normal has not unit length = " << d << endl;
+	}
+	d = Magnitude(dirIncoming); 
+	if ( (d < 0.99) ||
+		(d > 1.01) ) {
+			Debug << " The incoming dir has not unit length = " << d << endl;
+	}
+#endif
+
+	V = CrossProd(normal, dirIncoming);
+
+	if (SqrMagnitude(V) < 1e-3) {
+		// the normal and dirIncoming are colinear
+		// we can/have to generate arbitrary perpendicular vector to normal.
+		if (fabs(normal.x) < 0.6)
+			V.SetValue(0.0, -normal.z, normal.y);
+		else {
+			if (fabs(normal.y) < 0.6)
+				V.SetValue(-normal.z, 0.0, normal.x);
+			else
+				V.SetValue(-normal.y, normal.x, 0.0);
+		}
+	}
+	V = Normalize(V);
+
+	U = CrossProd(normal, V);
+#ifdef GTP_DEBUG
+	d = SqrMagnitude(U);
+	if ( (d < 0.99) ||
+		(d > 1.01) ) {
+			Debug << "The size of U vector incorrect\n";
+	}
+#endif
+	return;  
+}
+
+
+Vector3 UniformRandomVector()
+{
+	float r1 = RandomValue(0.0f, 1.0f);
+	float r2 = RandomValue(0.0f, 1.0f);
+
+	return UniformRandomVector(r1, r2);
+}
+
+
+Vector3 UniformRandomVector(const float r1, const float r2)
+{
+	float cosTheta = 1.0f - 2*r1;
+	float sinTheta = sqrt(1 - sqrt(cosTheta));
+	float fi = 2.0f * M_PI * r2;
+
+	Vector3 dir(sinTheta*sin(fi), cosTheta,	sinTheta * cos(fi));
+
+	return dir;
+}
+
+
+Vector3 CosineRandomVector(const Vector3 &normal)
+{
+	float r1, r2;
+
+	r1 = RandomValue(0.0f, 1.0f);
+	r2 = RandomValue(0.0f, 1.0f);
+
+	return CosineRandomVector(r1, r2, normal);
+}
+
+
+Vector3 CosineRandomVector(const float r1,
+						   const float r2,
+						   const Vector3 &normal)
+{
+	float theta = 2.0f * M_PI * r2;
+	float radius = sqrt(r1);
+
+	float x = radius * sin(theta);
+	float z = radius * cos(theta);
+	float y = sqrt(1.0f - x*x - z*z);
+
+	Vector3 dir(x, y, z);
+
+	//  return Normalize(dir);
+
+	Matrix4x4 m = RotationVectorsMatrix(normal, Vector3(0,1,0));
+
+	Matrix4x4 mi = Invert(m);
+
+	m = m * RotationVectorsMatrix(Vector3(0,1,0), Normalize(dir)) * mi;
+
+	return TransformNormal(m, normal);
+}
+
+
+Vector3 UniformRandomVector(const Vector3 &normal)
+{
+	float r1 = RandomValue(0.0f, 1.0f);
+	float r2 = RandomValue(0.0f, 1.0f);
+	
+	float cosTheta = 1.0f - r1;
+	float sinTheta = sqrt(1 - sqr(cosTheta));
+	float fi = 2.0f*M_PI*r2;
+
+	Vector3 dir(sinTheta*sin(fi),
+		cosTheta,
+		sinTheta*cos(fi));
+
+	Matrix4x4 m = RotationVectorsMatrix(
+		normal,
+		Vector3(0,1,0));
+	Matrix4x4 mi = Invert(m);
+	m = m*RotationVectorsMatrix(
+		Vector3(0,1,0),
+		Normalize(dir))*mi;
+
+	return TransformNormal(m, normal);
+}
+
+
+bool Vector3::CheckValidity() const
+{
+#ifdef _WIN32
+	return !(_isnan(x) || _isnan(y) || _isnan(z));
+#else
+	return !(isnan(x) || isnan(y) || isnan(z));
+#endif  
+	//return ((x != x) || (y != y) || (z != z));
+}
+
+//}
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/common.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/common.cpp	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/common.cpp	(revision 2746)
@@ -0,0 +1,339 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#include "common.h"
+
+#ifdef _X_WINDOW_
+#include "xwindow.h"
+#endif
+
+#ifndef _MSC_VER
+#include <unistd.h>
+#include <sys/types.h>
+#include <sys/time.h>
+#include <sys/times.h>
+#include <sys/resource.h>
+#else
+#include <windows.h>
+#include <sys/types.h>
+#include <sys/timeb.h>
+#endif
+
+
+using namespace std;
+
+//namespace CHCDemo {
+
+
+Real Limits::Threshold = 1e-6f;
+Real Limits::Small = 0.01f;
+Real Limits::Infinity = 1e20f;
+
+
+Real Random(Real max)
+{
+#ifdef __UNIX__
+  return (random()*max/0x7FFFFFFF);
+#else
+  return (((Real)rand()) / RAND_MAX) * max;
+#endif
+}
+
+
+int Random(int max)
+{
+	return rand() % max;
+}
+
+
+void Randomize()
+{
+	time_t t;
+	srand((unsigned) time(&t));
+}
+
+void
+Randomize(const unsigned int seed)
+{
+  srand(seed);
+}
+
+
+string
+ReplaceSuffix(const string &filename, const string &a, const string &b)
+{
+	string result = filename;
+
+	int pos = (int)filename.rfind(a, (int)filename.size() - 1);
+	if (pos == filename.size() - a.size()) 
+	{
+		result.replace(pos, a.size(), b);
+	}
+
+	return result;
+}
+
+
+int
+SplitFilenames(const string &str, vector<string> &filenames)
+{
+	int pos = 0;
+
+	while(1) {
+		int npos = (int)str.find(';', pos);
+		
+		if (npos < 0 || npos - pos < 1)
+			break;
+		filenames.push_back(string(str, pos, npos - pos));
+		pos = npos + 1;
+	}
+	
+	filenames.push_back(string(str, pos, str.size() - pos));
+	return (int)filenames.size();
+}
+
+void GetKey(char *s)
+{
+#ifdef _X_WINDOW_
+XWaitForButton(s);
+#else
+getch();
+#endif
+}
+
+#ifdef  _MSC_VER
+static bool hasHRTimer = false;
+static LARGE_INTEGER hrFreq;
+#endif
+
+
+void InitTiming()
+{
+#ifdef  _MSC_VER
+	hasHRTimer = (bool) QueryPerformanceFrequency(&hrFreq);
+#endif
+}
+
+long GetTime()
+{
+#ifndef  _MSC_VER
+	static struct rusage r;
+	getrusage(RUSAGE_SELF,&r);
+	return r.ru_utime.tv_usec+1000000*r.ru_utime.tv_sec;
+#else
+	if (hasHRTimer) {
+		LARGE_INTEGER counter;
+		QueryPerformanceCounter(&counter);
+		// return in usec
+		//return (long) (1000000*counter.QuadPart/(hrFreq.QuadPart));
+		// $$
+		// tmp store time in ms
+		return (long) (1000*counter.QuadPart/(hrFreq.QuadPart));
+	} else {
+		static struct _timeb mtime;
+		_ftime(&mtime);
+
+		return 1000*(1000*mtime.time + mtime.millitm);
+	}
+#endif
+}
+
+
+long
+GetRealTime()
+{
+#ifndef  _MSC_VER
+	static struct timeval _tstart;
+	static struct timezone tz;
+	gettimeofday(&_tstart,&tz);
+	return (long)(1000000*_tstart.tv_sec + _tstart.tv_usec);
+#else
+	if (hasHRTimer) {
+		LARGE_INTEGER counter;
+		QueryPerformanceCounter(&counter);
+		// return in usec
+		return (long) (1000000*counter.QuadPart/(hrFreq.QuadPart));
+	} else {
+		static struct _timeb mtime;
+		_ftime(&mtime);
+		return 1000*(1000*mtime.time + mtime.millitm);
+	}
+#endif
+}
+
+// return time diff. in ms
+Real TimeDiff(long time1, long time2) // in ms
+{
+	//const Real clk=1.0e-3f; // ticks per second
+	// $$ tmp store time in ms
+	const Real clk=1.0f; // ticks per second
+	long t=time2-time1;
+
+	return ((t<0)?-t:t)*clk;
+}
+
+char *
+TimeString()
+{
+	time_t t;
+	time(&t);
+	return ctime(&t);
+}
+
+char *
+GetAbsPath(char *name)
+{
+	//  char *url = new char[strlen(name)+1];
+	//  strcpy(url,name);
+
+	// get the absolute path
+	int c = strlen(name);
+	int i = c-1;
+	bool wasDot = false;
+	// search for double dot
+	for (;i>=0;i--) 
+		if (name[i]=='.') {
+			if (wasDot)
+				break;
+			wasDot = true;
+		} else
+			wasDot = false;
+
+		if (i>0)
+			i+=3;
+		if (i<0)
+			i=0;
+
+		char *url = new char[c-i+1];
+		int j=0;
+		for (;i<c;i++,j++)
+			url[j] = name[i];
+		url[j]=0;
+
+		return url;
+}
+
+
+void indent(ostream &app, int ind)
+{
+	int i;
+	
+	for (i = 0; i < ind; i++)
+		app << ' ';
+}
+
+
+char *GetPath(const char *s)
+{
+	int i=strlen(s);
+	for (; i>0; i--) {
+		if (s[i]=='/' || s[i]=='\\')
+			break;
+	}
+
+	char *path = new char[i+1];
+	int j=0;
+	for (; j<i; j++)
+		path[j] = s[j];
+	path[j] = 0;
+	return path;
+}
+
+
+char *
+StripPath(const char *s)
+{
+	const int len = strlen(s);
+
+	int i = len;
+	for (; i>0; i--) {
+		if (s[i]=='/' || s[i]=='\\')
+			break;
+	}
+
+	char *filename = new char[len - i];
+	int j = 0;
+
+	for (; j < len - i; ++ j)
+		filename[j] = s[i + j + 1];
+	//filename[j] = 0;
+
+	return filename;
+}
+
+char *
+strdup(char *a)
+{
+	if (a) {
+		char *s = new char[strlen(a)+1];
+		strcpy(s, a);
+		return s;
+	} else
+		return NULL;
+}
+
+bool
+FileExists(char *filename)
+{
+	FILE *f;
+
+	f = fopen(filename,"r");
+	if (f) {
+		fclose(f);
+		return true;
+	}
+	return false;
+}
+
+
+#if TOIMPLEMENT
+
+bool
+CreateDir(char *dir)
+{
+#ifdef _MSC_VER
+	HANDLE fFile; // File Handle
+	WIN32_FIND_DATA fileinfo; // File Information Structure
+	int n = strlen(dir) - 1;
+
+	if (dir[n] == '\\' || dir[n] == '/')
+		dir[n]=0;
+
+	fFile = FindFirstFile(dir, &fileinfo);
+
+	// if the file exists and it is a directory
+	if (fileinfo.dwFileAttributes == FILE_ATTRIBUTE_DIRECTORY) {
+		// Directory Exists close file and return
+		FindClose(fFile);
+		return true;
+	}
+
+	FindClose(fFile);
+
+	// Now lets cycle through the String Array and create each directory in turn
+	if (CreateDirectory(dir, 0)) {
+		SetFileAttributes(dir, FILE_ATTRIBUTE_NORMAL);
+		return true;
+	} else {
+		char *temp = GetPath(dir);
+
+		if (strcmp(temp, dir)!=0) {
+			CreateDir(temp);
+			if (CreateDirectory(dir, 0)) {
+				SetFileAttributes(dir, FILE_ATTRIBUTE_NORMAL);
+				return true;
+			}
+		}
+		delete temp;
+	}
+	Debug << "cannot create directory " << dir << endl;
+#endif
+	return false;
+}
+
+
+#endif
+
+//}
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.h	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.h	(revision 2746)
@@ -0,0 +1,535 @@
+#ifndef _Vector3_h__
+#define _Vector3_h__
+
+#include <iostream>
+#include "common.h"
+#include <math.h>
+
+
+// Forward-declare some other classes.
+class Matrix4x4;
+
+
+
+/** Vector class.
+*/
+class Vector3
+{
+public:
+
+	///////////
+	//-- members;
+
+	float x; 
+	float y;
+	float z;
+
+
+	/////////////////////////
+
+	void  SetX(float q) { x = q; }
+	void  SetY(float q) { y = q; }
+	void  SetZ(float q) { z = q; }
+
+	float GetX() const { return x; }
+	float GetY() const { return y; }
+	float GetZ() const { return z; }
+
+
+	//////////////
+	//-- constructors
+
+	Vector3() { }
+
+	Vector3(float x, float y, float z):
+	x(x), y(y), z(z)
+	{
+	}
+	Vector3(float v):
+	x(v), y(v), z(v)
+	{
+	}
+	/** Copy constructor.
+	*/
+	Vector3(const Vector3 &v): x(v.x), y(v.y), z(v.z)
+	{}
+
+	// Functions to get at the std::vector components
+	float& operator[] (const int inx) 
+	{
+		return (&x)[inx];
+	}
+
+	operator const float*() const 
+	{ 
+		return (const float*) this; 
+	}
+
+
+	const float &operator[] (const int inx) const 
+	{
+		return *(&x + inx);
+	}
+
+	void ExtractVerts(float *px, float *py, int which) const;
+  
+	void SetValue(float a, float b, float c)
+	{  
+		x = a; y = b; z = c; 
+	}
+
+	
+	void SetValue(float a) 
+	{ 
+		x = y = z = a; 
+	}
+ 
+	/** Returns the axis, where the std::vector has the largest value.
+	*/
+	int DrivingAxis(void) const;
+
+	/** returns the axis, where the std::vector has the smallest value
+	*/
+	int TinyAxis(void) const;
+	/** Returns largest component in this vector
+	*/
+	inline float MaxComponent(void) const 
+	{
+		return (x > y) ? ( (x > z) ? x : z) : ( (y > z) ? y : z);
+	}
+	
+	/** Returns copy of this vector where all components are positiv.
+	*/
+	inline Vector3 Abs(void) const 
+	{
+		return Vector3(fabs(x), fabs(y), fabs(z));
+	}
+
+	// normalizes the std::vector of unit size corresponding to given std::vector
+	inline void Normalize();
+
+	/** Returns false if this std::vector has a nan component.
+	*/
+	bool CheckValidity() const;
+
+	/**
+	===> Using ArbitraryNormal() for constructing coord systems 
+	===> is obsoleted by RightHandedBase() method (<JK> 12/20/03).   
+
+	Return an arbitrary normal to `v'.
+	In fact it tries v x (0,0,1) an if the result is too small,
+	it definitely does v x (0,1,0). It will always work for
+	non-degenareted std::vector and is much faster than to use
+	TangentVectors.
+
+	@param v(in) The std::vector we want to find normal for.
+	@return The normal std::vector to v.
+	*/
+	friend inline Vector3 ArbitraryNormal(const Vector3 &v); 
+
+	/**
+	Find a right handed coordinate system with (*this) being
+	the z-axis. For a right-handed system, U x V = (*this) holds.
+	This implementation is here to avoid inconsistence and confusion
+	when construction coordinate systems using ArbitraryNormal():
+	In fact:
+	V = ArbitraryNormal(N);
+	U = CrossProd(V,N);
+	constructs a right-handed coordinate system as well, BUT:
+	1) bugs can be introduced if one mistakenly constructs a 
+	left handed sytems e.g. by doing
+	U = ArbitraryNormal(N);
+	V = CrossProd(U,N);
+	2) this implementation gives non-negative base vectors
+	for (*this)==(0,0,1) |  (0,1,0) | (1,0,0), which is
+	good for debugging and is not the case with the implementation
+	using ArbitraryNormal()
+
+	===> Using ArbitraryNormal() for constructing coord systems 
+	is obsoleted by this method (<JK> 12/20/03).   
+	*/
+	void RightHandedBase(Vector3& U, Vector3& V) const;
+
+	/// Transforms a std::vector to the global coordinate frame.
+	/**
+	Given a local coordinate frame (U,V,N) (i.e. U,V,N are 
+	the x,y,z axes of the local coordinate system) and
+	a std::vector 'loc' in the local coordiante system, this
+	function returns a the coordinates of the same std::vector
+	in global frame (i.e. frame (1,0,0), (0,1,0), (0,0,1).
+	*/
+	friend inline Vector3 ToGlobalFrame(const Vector3& loc,
+		const Vector3& U,
+		const Vector3& V,
+		const Vector3& N);
+
+	/// Transforms a std::vector to a local coordinate frame.
+	/**
+	Given a local coordinate frame (U,V,N) (i.e. U,V,N are 
+	the x,y,z axes of the local coordinate system) and
+	a std::vector 'loc' in the global coordiante system, this
+	function returns a the coordinates of the same std::vector
+	in the local frame.
+	*/
+	friend inline Vector3 ToLocalFrame(const Vector3& loc,
+		const Vector3& U,
+		const Vector3& V,
+		const Vector3& N);
+
+	/// the magnitude=size of the std::vector
+	friend inline float Magnitude(const Vector3 &v);
+	/// the squared magnitude of the std::vector .. for efficiency in some cases
+	friend inline float SqrMagnitude(const Vector3 &v);
+	/// Magnitude(v1-v2)
+	friend inline float Distance(const Vector3 &v1, const Vector3 &v2);
+	/// SqrMagnitude(v1-v2)
+	friend inline float SqrDistance(const Vector3 &v1, const Vector3 &v2);
+
+	// creates the std::vector of unit size corresponding to given std::vector
+	friend inline Vector3 Normalize(const Vector3 &A);
+
+	// Rotate a normal std::vector.
+	friend Vector3 PlaneRotate(const Matrix4x4 &, const Vector3 &);
+
+	// construct view vectors .. DirAt is the main viewing direction
+	// Viewer is the coordinates of viewer location, UpL is the std::vector.
+	friend void ViewVectors(const Vector3 &DirAt, 
+							const Vector3 &Viewer,
+							const Vector3 &UpL, 
+							Vector3 &ViewV,
+							Vector3 &ViewU, 
+							Vector3 &ViewN);
+
+	// Given the intersection point `P', you have available normal `N'
+	// of unit length. Let us suppose the incoming ray has direction `D'.
+	// Then we can construct such two vectors `U' and `V' that
+	// `U',`N', and `D' are coplanar, and `V' is perpendicular
+	// to the vectors `N','D', and `V'. Then 'N', 'U', and 'V' create
+	// the orthonormal base in space R3.
+	friend void TangentVectors(Vector3 &U, 
+ 		                       Vector3 &V, // output
+							   const Vector3 &normal, // input
+							   const Vector3 &dirIncoming);
+
+	// Unary operators
+	Vector3 operator+ () const;
+	Vector3 operator- () const;
+
+	// Assignment operators
+	Vector3& operator+= (const Vector3 &A);
+	Vector3& operator-= (const Vector3 &A);
+	Vector3& operator*= (const Vector3 &A);
+	Vector3& operator*= (float A);
+	Vector3& operator/= (float A);
+
+	// Binary operators
+	friend inline Vector3 operator+ (const Vector3 &A, const Vector3 &B);
+	friend inline Vector3 operator- (const Vector3 &A, const Vector3 &B);
+	friend inline Vector3 operator* (const Vector3 &A, const Vector3 &B);
+	friend inline Vector3 operator* (const Vector3 &A, float B);
+	friend inline Vector3 operator* (float A, const Vector3 &B);
+	friend Vector3 operator* (const Matrix4x4 &, const Vector3 &);
+	friend inline Vector3 operator/ (const Vector3 &A, const Vector3 &B);
+
+	friend inline int operator< (const Vector3 &A, const Vector3 &B);
+	friend inline int operator<= (const Vector3 &A, const Vector3 &B);
+
+	friend inline Vector3 operator/ (const Vector3 &A, float B);
+	friend inline int operator== (const Vector3 &A, const Vector3 &B);
+	friend inline float DotProd(const Vector3 &A, const Vector3 &B);
+	friend inline Vector3 CrossProd (const Vector3 &A, const Vector3 &B);
+
+	friend std::ostream& operator<< (std::ostream &s, const Vector3 &A);
+	friend std::istream& operator>> (std::istream &s, Vector3 &A);
+
+	friend void Minimize(Vector3 &min, const Vector3 &candidate);
+	friend void Maximize(Vector3 &max, const Vector3 &candidate);
+
+	friend inline int EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2, float thr);
+	friend inline int EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2);
+
+	friend Vector3 CosineRandomVector(const Vector3 &normal);
+	friend Vector3 UniformRandomVector(const Vector3 &normal);
+	friend Vector3 UniformRandomVector();
+
+	friend Vector3 UniformRandomVector(float r1, float r2);
+
+	friend Vector3 CosineRandomVector(float r1, float r2, const Vector3 &normal);
+};
+
+
+
+// forward declaration
+extern Vector3 UniformRandomVector(const Vector3 &normal);
+extern Vector3 UniformRandomVector();
+extern Vector3 UniformRandomVector(const float r1, const float r2);
+
+inline Vector3 ArbitraryNormal(const Vector3 &N)
+{
+	float dist2 = N.x * N.x + N.y * N.y;
+
+	if (dist2 > 0.0001f) 
+	{
+		float inv_size = 1.0f / sqrtf(dist2);
+		return Vector3(N.y * inv_size, -N.x * inv_size, 0); // N x (0,0,1)
+	}
+
+	float inv_size = 1.0f / sqrtf(N.z * N.z + N.x * N.x);
+	return Vector3(-N.z * inv_size, 0, N.x * inv_size); // N x (0,1,0)
+}
+
+
+inline void Vector3::RightHandedBase(Vector3& U, Vector3& V) const
+{
+	// HACK
+	V = ArbitraryNormal(*this);
+	U = CrossProd(V, *this);
+}
+
+
+inline Vector3 ToGlobalFrame(const Vector3 &loc,
+							 const Vector3 &U,
+							 const Vector3 &V,
+							 const Vector3 &N)
+{
+	return loc.x * U + loc.y * V + loc.z * N;
+}
+
+
+inline Vector3 ToLocalFrame(const Vector3 &loc,
+							const Vector3 &U,
+							const Vector3 &V,
+							const Vector3 &N)
+{
+	return Vector3(loc.x * U.x + loc.y * U.y + loc.z * U.z,
+				   loc.x * V.x + loc.y * V.y + loc.z * V.z,
+				   loc.x * N.x + loc.y * N.y + loc.z * N.z);
+}
+
+
+inline float Magnitude(const Vector3 &v) 
+{
+	return sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
+}
+
+
+inline float SqrMagnitude(const Vector3 &v) 
+{
+  return v.x * v.x + v.y * v.y + v.z * v.z;
+}
+
+
+inline float Distance(const Vector3 &v1, const Vector3 &v2)
+{
+  return sqrtf(sqrt(v1.x - v2.x) + sqrt(v1.y - v2.y) + sqrt(v1.z - v2.z));
+}
+
+
+inline float SqrDistance(const Vector3 &v1, const Vector3 &v2)
+{
+  return sqrt(v1.x - v2.x) + sqrt(v1.y - v2.y) + sqrt(v1.z - v2.z);
+}
+
+
+inline Vector3 Normalize(const Vector3 &A)
+{
+  return A * (1.0f / Magnitude(A));
+}
+
+
+inline float DotProd(const Vector3 &A, const Vector3 &B)
+{
+  return A.x * B.x + A.y * B.y + A.z * B.z;
+}
+
+
+// angle between two vectors with respect to a surface normal in the
+// range [0 .. 2 * pi]
+inline float Angle(const Vector3 &A, const Vector3 &B, const Vector3 &norm)
+{
+	Vector3 cross = CrossProd(A, B);
+
+	float signedAngle;
+
+	if (DotProd(cross, norm) > 0)
+		signedAngle = atan2(-Magnitude(CrossProd(A, B)), DotProd(A, B));
+	else
+		signedAngle = atan2(Magnitude(CrossProd(A, B)), DotProd(A, B));
+
+	if (signedAngle < 0)
+		return 2 * PI + signedAngle;
+
+	return signedAngle;
+}
+
+
+inline Vector3 Vector3::operator+() const
+{
+	return *this;
+}
+
+
+inline Vector3 Vector3::operator-() const
+{
+  return Vector3(-x, -y, -z);
+}
+
+
+inline Vector3 &Vector3::operator+=(const Vector3 &A)
+{
+	x += A.x;  y += A.y;  z += A.z;
+	return *this;
+}
+
+
+inline Vector3& Vector3::operator-=(const Vector3 &A)
+{
+	x -= A.x;  y -= A.y;  z -= A.z;
+	
+	return *this;
+}
+
+
+inline Vector3& Vector3::operator*= (float A)
+{
+	x *= A;  y *= A;  z *= A;
+	return *this;
+}
+
+
+inline Vector3& Vector3::operator/=(float A)
+{
+	float a = 1.0f / A;
+	
+	x *= a;  y *= a;  z *= a;
+	
+	return *this;
+}
+
+
+inline Vector3& Vector3::operator*= (const Vector3 &A)
+{
+  x *= A.x;  y *= A.y;  z *= A.z;
+  return *this;
+}
+
+
+inline Vector3 operator+ (const Vector3 &A, const Vector3 &B)
+{
+	return Vector3(A.x + B.x, A.y + B.y, A.z + B.z);
+}
+
+
+inline Vector3 operator- (const Vector3 &A, const Vector3 &B)
+{
+	return Vector3(A.x - B.x, A.y - B.y, A.z - B.z);
+}
+
+
+inline Vector3 operator* (const Vector3 &A, const Vector3 &B)
+{
+	return Vector3(A.x * B.x, A.y * B.y, A.z * B.z);
+}
+
+
+inline Vector3 operator* (const Vector3 &A, float B)
+{
+  return Vector3(A.x * B, A.y * B, A.z * B);
+}
+
+
+inline Vector3 operator* (float A, const Vector3 &B)
+{
+	return Vector3(B.x * A, B.y * A, B.z * A);
+}
+
+
+inline Vector3 operator/ (const Vector3 &A, const Vector3 &B)
+{
+	return Vector3(A.x / B.x, A.y / B.y, A.z / B.z);
+}
+
+
+inline Vector3 operator/ (const Vector3 &A, float B)
+{
+	float b = 1.0f / B;
+	return Vector3(A.x * b, A.y * b, A.z * b);
+}
+
+
+inline int operator< (const Vector3 &A, const Vector3 &B)
+{
+	return A.x < B.x && A.y < B.y && A.z < B.z;
+}
+
+
+inline int operator<= (const Vector3 &A, const Vector3 &B)
+{
+	return A.x <= B.x && A.y <= B.y && A.z <= B.z;
+}
+
+
+// Might replace floating-point == with comparisons of
+// magnitudes, if needed.
+inline int operator== (const Vector3 &A, const Vector3 &B)
+{ 
+	return (A.x == B.x) && (A.y == B.y) && (A.z == B.z);
+}
+
+
+inline Vector3 CrossProd (const Vector3 &A, const Vector3 &B)
+{
+	return Vector3(A.y * B.z - A.z * B.y,
+		           A.z * B.x - A.x * B.z,
+				   A.x * B.y - A.y * B.x);
+}
+
+
+inline void Vector3::Normalize()
+{
+	float sqrmag = x * x + y * y + z * z;
+	
+	if (sqrmag > 0.0f)
+		(*this) *= 1.0f / sqrtf(sqrmag);
+}
+
+
+// Overload << operator for C++-style output
+inline std::ostream& operator<< (std::ostream &s, const Vector3 &A)
+{
+	return s << "(" << A.x << ", " << A.y << ", " << A.z << ")";
+}
+
+
+// Overload >> operator for C++-style input
+inline std::istream& operator>> (std::istream &s, Vector3 &A)
+{
+	char a;
+	// read "(x, y, z)"
+	return s >> a >> A.x >> a >> A.y >> a >> A.z >> a;
+}
+
+
+inline int EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2, float thr)
+{
+	if (fabsf(v1.x-v2.x) > thr)
+		return false;
+	if (fabsf(v1.y-v2.y) > thr)
+		return false;
+	if (fabsf(v1.z-v2.z) > thr)
+		return false;
+	
+	return true;
+}
+
+
+inline int EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2)
+{
+	return EpsilonEqualV3(v1, v2, Limits::Small);
+}
+
+
+
+//}
+
+#endif
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/occquery.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/occquery.cpp	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/occquery.cpp	(revision 2746)
@@ -1,17 +1,11 @@
 // occquery.cpp : Defines the entry point for the console application.
 //
-
-//#include "stdafx.h"
+#if TOIMPLEMENT
+
+#include <math.h>
+#include <time.h>
 #include "glInterface.h"
 #include "RenderTraverser.h"
 
-extern "C"
-{
-	#include "MathStuff.h"
-	#include "DataTypes.h"
-}
-
-#include <math.h>
-#include <time.h>
 
 
@@ -109,7 +103,9 @@
 void deleteGeometry();
 
+#endif
 
 int main(int argc, char* argv[])
 {
+#if 0
 	glutInitWindowSize(800,600);
 	glutInit(&argc,argv);
@@ -141,9 +137,9 @@
 	// clean up
 	cleanUp();
-	
+#endif
 	return 0;
 }
 
-
+#if 0
 void initGLstate(void) 
 {
@@ -850,2 +846,4 @@
 	}
 }
+
+#endif
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/Geometry.cpp
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/Geometry.cpp	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/Geometry.cpp	(revision 2746)
@@ -89,4 +89,5 @@
 	return true;
 }
+
 
 void Geometry::GenerateList()
@@ -112,7 +113,9 @@
 				break;
 		}
+
 		glEndList();
 	}
 }
+
 
 void Geometry::Render()
@@ -238,18 +241,21 @@
 	{
 		glBegin(GL_TRIANGLE_STRIP);
-			while(teapot_indices[i] != STRIP_END)
-			{	
-				int index = teapot_indices[i] * 3;
-				
-				glNormal3fv(teapot_normals + index);
-				glVertex3fv(teapot_vertices + index);
-				
-				i++;		
-			}
+		
+		while (teapot_indices[i] != STRIP_END)
+		{	
+			int index = teapot_indices[i] * 3;
+
+			glNormal3fv(teapot_normals + index);
+			glVertex3fv(teapot_vertices + index);
+			
+			++ i;		
+		}
+
 		glEnd();
-		
-		i++; // skip strip end flag
-	}
-}
+
+		++ i; // skip strip end flag
+	}
+}
+
 
 void Geometry::CalcBoundingVolume()
@@ -276,5 +282,5 @@
 	Vector3 *rotatedPoints = new Vector3[num_vertices];
 	
-	for(int i = 0; i < num_vertices; i++)
+	for (int i = 0; i < num_vertices; ++ i)
 	{
 		copyVector3Values(rotatedPoints[i], vertices[i * 3],
@@ -288,5 +294,5 @@
 	calcCubicHull(mBoundingBox.min, mBoundingBox.max, rotatedPoints, num_vertices);
 
-	for(int i=0; i< 3; i++)
+	for (int i = 0; i < 3; ++ i)
 	{
 		mBoundingBox.min[i] *= mScale;
@@ -325,6 +331,5 @@
 }
 
-/*
-	creates a torus with specified radii, with number of rings and ring
+/*	Creates a torus with specified radii, with number of rings and ring
 	subdivision = precision.
 */
@@ -338,5 +343,5 @@
 	torus_indices = new int[num_torus_indices];
 
-	for(int i = 0; i < precision + 1; ++ i)
+	for (int i = 0; i < precision + 1; ++ i)
 	{
 		int index = i * 3;
@@ -364,7 +369,7 @@
 	}
 
-	for(int i_rng = 1; i_rng < precision + 1; ++ i_rng)
-	{
-		for(int i=0; i<precision+1; ++i)
+	for (int i_rng = 1; i_rng < precision + 1; ++ i_rng)
+	{
+		for(int i = 0; i<precision + 1; ++ i)
 		{
 			int index = 3 * (i_rng * (precision + 1) + i);
@@ -388,7 +393,7 @@
 	}
 
-	for(int i_rng = 0; i_rng < precision; ++ i_rng)
-	{
-		for(int i = 0; i < precision; ++ i)
+	for (int i_rng = 0; i_rng < precision; ++ i_rng)
+	{
+		for (int i = 0; i < precision; ++ i)
 		{
 			int index = ((i_rng * precision + i) * 2) * 3;
@@ -407,6 +412,6 @@
 }
 
-/**
-	counts the triangles of the teapot.
+
+/** Counts the triangles of the teapot.
 	traverses through all the triangle strips.
 */
@@ -414,21 +419,23 @@
 {
 	int result = 0;
-	int i=0;
+	int i = 0;
+
 	// n - 2 triangles are drawn for a strip
-	while(i < num_teapot_indices)
-	{
-		while(teapot_indices[i] != STRIP_END)
+	while (i < num_teapot_indices)
+	{
+		while (teapot_indices[i] != STRIP_END)
 		{	
 			result ++;;
-			i++;		
+			++ i;		
 		}
 		result -= 2;
-		i++; // skip STRIP_END
-	}
+		++ i; // skip STRIP_END
+	}
+
 	return result;
 }
 
-/**
-	counts the triangles of the torus
+
+/**	Counts the triangles of the torus
 */
 int Geometry::CountTorusTriangles()
@@ -452,5 +459,5 @@
 	int result = 0;
 
-	switch(objectType)
+	switch (objectType)
 	{
 		case TEAPOT:
@@ -466,4 +473,5 @@
 			break;
 	}
+
 	return result;
 }
@@ -497,7 +505,7 @@
 }
 
-/**
-   renders a sphere with sphere_precision subdivisions.
-   note: works only for even sphere_precision
+
+/** Renders a sphere with sphere_precision subdivisions. 
+    Note: works only for even sphere_precision
 */
 void Geometry::RenderSphere()
@@ -508,5 +516,5 @@
 	glFrontFace(GL_CW);
 
-	for (int j = 0; j < sphere_precision/2; ++j) 
+	for (int j = 0; j < sphere_precision / 2; ++ j) 
 	{
 		double alpha = j * PI * 2.0 / (double)sphere_precision - PI_2;
@@ -515,5 +523,5 @@
 		glBegin(GL_TRIANGLE_STRIP);
 
-		for (int i = 0; i <= sphere_precision; ++i) 
+		for (int i = 0; i <= sphere_precision; ++ i) 
 		{
 			double gamma = i * PI * 2.0 / (double)sphere_precision;
@@ -536,4 +544,5 @@
 		glEnd();
 	}
+
 	glFrontFace(GL_CCW);
 }
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/chc_revisited.vcproj
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/chc_revisited.vcproj	(revision 2642)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/chc_revisited.vcproj	(revision 2746)
@@ -182,9 +182,9 @@
 			>
 			<File
-				RelativePath=".\DataTypes.c"
-				>
-			</File>
-			<File
-				RelativePath=".\Geometry.cpp"
+				RelativePath=".\AxisAlignedBox3.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\common.cpp"
 				>
 			</File>
@@ -194,5 +194,9 @@
 			</File>
 			<File
-				RelativePath=".\MathStuff.c"
+				RelativePath=".\Material.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\Matrix4x4.cpp"
 				>
 			</File>
@@ -216,9 +220,17 @@
 			>
 			<File
-				RelativePath=".\DataTypes.h"
-				>
-			</File>
-			<File
-				RelativePath=".\Geometry.h"
+				RelativePath=".\AxisAlignedBox3.h"
+				>
+			</File>
+			<File
+				RelativePath=".\Bvh.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\Bvh.h"
+				>
+			</File>
+			<File
+				RelativePath=".\common.h"
 				>
 			</File>
@@ -228,9 +240,9 @@
 			</File>
 			<File
-				RelativePath=".\HierarchyNode.h"
-				>
-			</File>
-			<File
-				RelativePath=".\MathStuff.h"
+				RelativePath=".\Material.h"
+				>
+			</File>
+			<File
+				RelativePath=".\Matrix4x4.h"
 				>
 			</File>
@@ -240,9 +252,13 @@
 			</File>
 			<File
-				RelativePath=".\teapot.h"
-				>
-			</File>
-			<File
 				RelativePath=".\Timers.h"
+				>
+			</File>
+			<File
+				RelativePath=".\Vector3.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\Vector3.h"
 				>
 			</File>
Index: GTP/trunk/App/Demos/Vis/CHC_revisited/common.h
===================================================================
--- GTP/trunk/App/Demos/Vis/CHC_revisited/common.h	(revision 2746)
+++ GTP/trunk/App/Demos/Vis/CHC_revisited/common.h	(revision 2746)
@@ -0,0 +1,500 @@
+/**
+
+This file contains various macros, templates and constants for the ERS system.
+
+@author Jiri Bittner
+*/
+
+#ifndef __COMMON_H
+#define __COMMON_H
+
+
+#include <math.h>
+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <limits.h>
+#include <vector>
+
+
+//namespace GtpVisibilityPreprocessor {
+
+  
+#if defined(_MSC_VER)
+// use perftimer only on msvc
+#define USE_PERFTIMER
+
+// define __SSE__ macro as it is not defined under MSVC
+#define __SSE__
+#endif
+
+// This constant should be used for the length of the array char for filenames
+// etc., for example: char filename[MaxStringLength]
+const int MaxStringLength = 256;
+
+#if defined(_MSC_VER)
+#pragma warning(disable:4018)
+#pragma warning(disable:4800)
+//#pragma warning(disable:4244)
+
+#if 0 // Note matt: comment this out because conflicts with definition in qt library!!
+typedef unsigned int uint;
+#endif
+typedef unsigned short ushort;
+typedef unsigned char uchar;
+typedef unsigned long ulong;
+#endif
+
+#if defined(__GNUC__) || defined(_MSC_VER)
+//#define DIRCAT '.'
+#endif
+
+#if !defined(__WATCOMC__) && !defined(__CYGWIN32__) && !defined(_MSC_VER)
+#include <values.h>
+#else // __WATCOMC__
+#define M_PI        3.14159265358979323846
+#define MAXFLOAT    3.40282347e+37F
+#endif // __WATCOMC__
+
+// some compilers do not define the bool yet, but it was declared by ANSI
+#if !defined(__WATCOMC__) && !defined(_MSC_VER)
+#if defined (__GNUC__) || (_BOOL)
+//#error "HAS BOOL defined"
+#define HAS_BOOL
+#endif
+#else // __WATCOMC__
+#if  (__WATCOMC__ > 1060)
+//#error "Watcom HAS BOOL defined"
+#define HAS_BOOL
+#endif
+#endif // __WATCOMC__
+
+
+#if defined(__WATCOMC__) || defined(_MSC_VER)
+#define strcasecmp stricmp
+#define strncasecmp strnicmp
+#endif // __WATCOMC__
+
+// matt
+
+#define USE_GZLIB 1
+
+#if USE_GZLIB
+
+#define OUT_BIN_MODE ios::out
+#define IN_BIN_MODE ios::in
+
+#else
+
+#ifdef sgi
+#define OUT_BIN_MODE ios::out
+#define IN_BIN_MODE ios::in
+#else // sgi
+#if defined(__WATCOMC__) || defined(_MSC_VER)
+#define OUT_BIN_MODE ios::out | ios::binary
+#define IN_BIN_MODE ios::in | ios::binary
+#else
+#define OUT_BIN_MODE ios::out | ios::bin
+#define IN_BIN_MODE ios::in | ios::bin
+#endif // __WATCOMC_
+#endif // sgi
+
+#endif
+
+//  #ifndef HAS_BOOL
+//  //enum bool {
+//  //  false = 0,
+//  //  true
+//  //};
+//  #define bool int
+//  #define false 0
+//  #define true  1
+//  #endif // HAS_BOOL
+
+typedef unsigned long dword;
+
+#ifndef NULL
+#define NULL (void *)0
+#endif // NULL
+
+//ostream& operator<<(ostream &s,const BaseC &c);
+
+typedef float Real;
+typedef unsigned char byte;
+
+//#ifndef FALSE
+//#define FALSE 0
+//#define TRUE !0
+//#endif
+
+#ifndef __GNUG__
+// typedef int bool;
+#endif
+
+#ifndef getch
+#define getch() getchar()
+#endif
+
+
+#define TRASH 1.0e-5
+
+#ifndef PI
+#define PI    3.14159265358979323846f
+#endif
+
+#define MIN_FLOAT -1e30f
+#define MAX_FLOAT  1e30f
+
+
+// tolerance value for polygon area
+#define AREA_LIMIT 0.0001f
+
+#ifndef DEL_PTR
+#define DEL_PTR(ptr) do {if (ptr) {        \
+                           delete (ptr);   \
+						   (ptr) = 0;}}    \
+                     while (0)
+#endif
+// Clears a container (i.e., a vector of pointers) and deletes the pointers
+#if 0
+#ifndef CLEAR_CONTAINER
+#define CLEAR_CONTAINER(co) do { while (!(co).empty()) {   \
+						           DEL_PTR((co).back());     \
+								   (co).pop_back();}}      \
+							while (0)
+#endif
+
+#else
+
+#ifndef CLEAR_CONTAINER
+#define CLEAR_CONTAINER(co) do { for (int _i = 0; _i < (int)(co).size(); ++ _i) { \
+	DEL_PTR((co)[_i]);} \
+	(co).clear(); } \
+while (0)
+#endif
+
+#endif
+
+
+inline
+int signum(const Real a, const Real thresh = TRASH)
+{
+  if (a>thresh)
+    return 1; 
+  else
+    if (a<-thresh)
+      return -1; 
+  return 0;
+}
+
+inline double Absd(const double a)
+{
+  return (a >= 0.0) ? a : -a;
+}
+
+inline float Abs(const float a)
+{
+  return (a >= 0.0f) ? a : -a;
+}
+
+// =======================================================
+// Comparing things
+//struct Limits {
+//  const Real thresh=TRASH;
+//  const Real small=0.1;
+//};
+
+template <class T>
+bool
+ClipValue(T &v, const T m, const T M)
+{
+  if (v<m) {
+    v = m;
+    return true;
+  }
+  if (v>M) {
+    v = M;
+    return true;
+  }
+
+  return false;
+}
+
+
+inline int eq(Real a, Real b, Real t=TRASH)
+{
+  return Abs(a-b)<t;
+}
+
+inline int leq(Real a,Real b,Real t=TRASH)
+{
+  return a - b < t;
+}
+
+inline int geq(Real a,Real b,Real t=TRASH)
+{
+  return t > b - a;
+}
+
+inline int le(Real a,Real b,Real t=TRASH)
+{
+  return !geq(a,b,t);
+}
+
+inline int ge(Real a,Real b,Real t=TRASH)
+{
+  return !leq(a,b,t);
+}
+
+// ========================================================
+
+// -------------------------------------------------------------------
+// Indents to a given stream by the number of spaces specified.
+// This routine is located in main.cpp, for lack of a better place.
+// -------------------------------------------------------------------
+void indent(std::ostream &app, int ind);
+
+// ---------------------------------------------------------
+// RandomValue
+//	Returns a random Realing-point value between the two 
+//	values.  Range is inclusive; the function should 
+//	occasionally return exactly a or b.
+// ---------------------------------------------------------
+inline Real RandomValue(Real a, Real b)
+{
+  Real range = (Real) Abs(a - b);
+  return ((Real)rand() / RAND_MAX) * range + ((a < b) ? a : b);
+}
+
+
+/*inline int
+RandomValue(int a, int b)
+{
+	int range = abs(a - b);
+	return (rand() * range) / RAND_MAX + ((a < b) ? a : b);
+}*/
+
+inline Real sqr(Real a)
+{
+  return a*a;
+}
+
+
+template <class T>
+void Swap(T &a,T &b)
+{
+  T c;                 
+  c = b;
+  b = a;
+  a = c;
+}
+
+template <class T>
+int eq(T &a, T &b, T &c, T &d) {
+  return a == b && c==d && b==c;
+}
+
+template <class T>
+T Min(T a,T b)
+{
+  return a<b ? a : b;
+}
+
+template <class T>
+T Max(T a,T b)
+{
+  return a>b ? a : b;
+}
+
+
+Real Random(Real max);
+int  Random(int max);
+void Randomize();
+void
+Randomize(const unsigned int seed);
+
+
+void GetKey(char *s=NULL);
+
+inline Real Deg2Rad(const Real a)
+{
+   return a * (PI / 180.0f);
+}
+
+inline Real Rad2Deg(const Real a) {
+  return a*(180.0f/PI);
+}
+  
+void InitTiming();
+long GetTime();
+long GetRealTime();
+Real TimeDiff(long t1, long t2);
+char *TimeString();
+
+
+
+// -------------------------------------------------------------------
+// Limits.
+//  This class encapsulates all the concessions to Realing-point
+//  error made by ray tracers.
+// -------------------------------------------------------------------
+
+class Limits {
+public:
+  // This is the number used to reject too-close intersections.
+  // The default value is 1.0.
+  static Real Threshold;
+
+  // This is a "small" number.  Less than this number is assumed to
+  // be 0, when such things matter.
+  // The default value is 0.1.
+  static Real Small;
+
+  // This is an impractically "large" number, used for intersection
+  // parameters out to infinity (e.g. the span resulting from an
+  // FindAllIntersections operation on a plane).
+  // The default value is 100000.
+  static Real Infinity;
+};
+
+// ---------------------------------------------------------
+// EpsilonEqual(x,y)
+//   Returns if two values are equal or not (by epsilon)
+// ---------------------------------------------------------
+inline int
+EpsilonEqual(const Real &x, const Real &y)
+{
+  return fabs(x-y) < Limits::Small;
+}
+
+// ---------------------------------------------------------
+// EpsilonEqual(x)
+//   Returns if a value is zero (+/- epsilon)
+// ---------------------------------------------------------
+inline int
+EpsilonEqual(const Real &x)
+{
+  return fabs(x) < Limits::Small;
+}
+
+// ---------------------------------------------------------
+// EpsilonEqual(x,y,epsilon)
+//   Returns if two values are equal or not by a given epsilon
+// ---------------------------------------------------------
+inline int
+EpsilonEqual(const Real &x, const Real &y, const Real &epsilon)
+{
+  return fabs(x-y) < epsilon;
+}
+
+// ---------------------------------------------------------
+// InRange
+//	Returns nonzero if min <= candidate <= max.
+// ---------------------------------------------------------
+template<class T>
+inline int 
+InRange(T min, T max, T candidate)
+{
+  return (candidate >= min) && (candidate <= max);
+}
+
+// --------------------------------------------------------------
+// string function with new operator
+
+inline char*
+StrDup(char *src) {
+  char *p;
+  for (p = src; *p++; );
+  char *dest = new char[p-src];
+  for ( p = dest;(*p++ = *src++) != 0; );
+  return dest;
+}
+
+
+inline char *
+StrToLower(char *src) {
+  char *p;
+  for (p = src; *p; p++)
+    *p = tolower(*p);
+  return src;
+}
+
+// return l = log2(a) if a is a power of two else -ceillog2(a)
+inline int
+GetLog2(int a)
+{
+  int i, x;
+  i = 0;
+  x = 1;
+
+  while (x < a) {
+    i++;
+    x <<= 1;
+  }
+
+  return (x==a) ? i: -i;
+}
+
+
+// return ceil(log2(a)) even if a is not a power of two
+// Example:
+// GetCeilLog2(15) = 4
+// GetCeilLog2(16) = 4
+// GetCeilLog2(17) = 5
+inline int
+GetCeilLog2(int a)
+{
+  int i, x;
+  if (a < 0)
+    return -1;
+
+  i = 0;
+  x = 1;
+
+  while (x < a) {
+    i++;
+    x <<= 1;
+  }
+
+  return i;
+}
+
+
+char *
+GetAbsPath(char *path);
+
+char *
+strdup(char *a);
+
+std::string
+ReplaceSuffix(const std::string &filename, const std::string &a, const std::string &b);
+
+
+int
+SplitFilenames(const std::string &str, std::vector<std::string> &filenames);
+
+bool
+FileExists(char *filename);
+
+bool
+CreateDir(char *dir);
+
+char *
+GetPath(const char *s);
+
+char *
+StripPath(const char *s);
+
+
+//}
+
+#endif
+
+
+
+
+
+
+
+
