Index: trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h	(revision 177)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.h	(revision 191)
@@ -5,4 +5,5 @@
 
 class KdNode;
+class Ray;
 
 struct LtKdNode
@@ -29,4 +30,7 @@
   int mSamples;
   KdPvsMap mEntries;
+  
+  KdPvs():mSamples(0),mEntries() {}
+  
   KdPvsData *Find(KdNode *node);
   int AddNodeSample(KdNode *node);
@@ -40,4 +44,6 @@
 };
 
+
+
 #endif
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.h	(revision 177)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.h	(revision 191)
@@ -8,4 +8,5 @@
 #include "AxisAlignedBox3.h"
 #include "Ray.h"
+#include "Pvs.h"
 
   
@@ -159,8 +160,14 @@
     mObjects.reserve(objects);
   }
-
+  
+  void AddPassingRay(const Ray &ray, const int contributions) {
+    mPassingRays.AddRay(ray, contributions);
+  }
   
   /** \sa KdNode::IsLeaf() */
   virtual bool IsLeaf() const { return true; }
+  
+
+  
 
   /** pointers to occluders contained in this node */
@@ -169,4 +176,11 @@
   /** pointers to viewcells contained in this node */
   //  ViewCellContainer mViewCells;
+
+  /** Ray set description of the rays passing through this node */
+  PassingRaySet mPassingRays;
+
+  /** PVS consisting of visible KdTree nodes */
+  KdPvs mKdPvs;
+  
 };
   
@@ -268,6 +282,9 @@
   void
   CollectObjects(KdNode *n, ObjectContainer &objects);
-  
-    AxisAlignedBox3 GetBox(const KdNode *node) const {
+
+  void
+  CollectLeaves(vector<KdLeaf *> &leaves);
+
+  AxisAlignedBox3 GetBox(const KdNode *node) const {
     KdInterior *parent = node->mParent;
     if (parent == NULL)
@@ -295,5 +312,8 @@
 		bool onlyUnmailed
 		);
-  
+
+  int
+  CollectLeafPvs();
+
 protected:
 
@@ -469,5 +489,5 @@
   AxisAlignedBox3 mBox;
   KdTreeStatistics mStat;
-
+  
 };
   
Index: trunk/VUT/GtpVisibilityPreprocessor/src/MutualVisibility.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/MutualVisibility.h	(revision 191)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/MutualVisibility.h	(revision 191)
@@ -0,0 +1,138 @@
+#ifndef __MUTUAL_VISIBILITY_H
+#define __MUTUAL_VISIBILITY_H
+
+#include "Vector3.h"
+#include "Ray.h"
+#include "Triangle3.h"
+
+class Intersectable;
+class AxisAlignedBox3;
+
+
+
+
+struct SimpleRay {
+  Vector3 mOrigin;
+  Vector3 mDirection;
+};
+
+
+struct TriangleSample {
+public:
+  /// evaluted sampling error
+  float mError;
+  /// depth in recursion
+  int mDepth;
+  /// The source triangle
+  Triangle3 mSource;
+  /// The target triangle
+  Triangle3 mTarget;
+  
+  /// intersections with the scene
+  vector<Ray::Intersection> mIntersections[9];
+  
+  void ComputeError();
+  
+  TriangleSample () {}
+  
+  TriangleSample (
+		  const Triangle3 &source,
+		  const Triangle3 &target)
+  {
+    Init(source, target);
+  }
+
+  // initial triangle sample
+  void Init(
+	    const Triangle3 &source,
+	    const Triangle3 &target);
+  
+  Vector3
+  GetIntersectionPoint(const int rayIndex,
+		       const int depth) const;
+
+  void GetRay(const int rayIndex,
+	      Vector3 &origin,
+	      Vector3 &direction) const;
+
+};
+
+enum { VISIBLE, INVISIBLE };
+
+  
+class MutualVisibilitySampler {
+public:
+  KdTree *mKdTree;
+  AxisAlignedBox3 mSource;
+  AxisAlignedBox3 mTarget;
+  float mSolidAngleThreshold;
+
+
+  MutualVisibilitySampler(KdTree *kdTree,
+			  AxisAlignedBox3 &source,
+			  AxisAlignedBox3 &target,
+			  const float solidAngleThreshold);
+  int
+  ComputeVisibility();
+
+  void
+  ConstructInitialSamples(
+			  const AxisAlignedBox3 &source,
+			  const AxisAlignedBox3 &target,
+			  vector<TriangleSample *> &samples
+			  );
+
+  void
+  AddInitialSamples(
+		    const Rectangle3 &sourceRect,
+		    const Rectangle3 &targetRect,
+		    vector<TriangleSample *> &samples
+		    );
+
+  // the split sample method contains a methodology to create new samples
+  // or terminate the sampling
+  bool
+  SplitSample(
+	      const TriangleSample &source,
+	      TriangleSample &sample1,
+	      TriangleSample &sample2
+	      );
+  void
+  PerformSplit(
+	       const TriangleSample &sample,
+	       const bool splitSource,
+	       const int edge,
+	       const Vector3 &splitPoint,
+	       const Ray &ray,
+	       TriangleSample &sample1,
+	       TriangleSample &sample2
+	       );
+  
+
+  bool
+  SampleTerminationCriteriaMet(
+			       const TriangleSample &sample);
+  
+  float
+  GetSpatialAngle(const TriangleSample &sample,
+		  const Vector3 &point
+		  );
+
+  void
+  ComputeError(TriangleSample &sample);
+
+
+};
+  
+int
+ComputeBoxVisibility(KdTree *kdTree,
+		     AxisAlignedBox3 &source,
+		     AxisAlignedBox3 &target,
+		     float solidAngleThreshold);
+
+
+
+
+
+
+#endif
Index: trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/main.cpp	(revision 191)
@@ -45,5 +45,6 @@
     Camera camera;
     camera.LookAtBox(p->mKdTree->GetBox());
-    camera.SnapImage("camera.png", p->mKdTree);
+    camera.SnapImage("camera.jpg", p->mKdTree);
+
     
     camera.LookInBox(p->mKdTree->GetBox());
Index: trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.h	(revision 191)
@@ -24,5 +24,8 @@
 
   int
-  AddNodeSamples(Intersectable *object, const Ray &ray);
+  AddNodeSamples(Intersectable *object, const Ray &ray, const int pass);
+
+  void
+  HoleSamplingPass();
 
 };
Index: trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.cpp	(revision 191)
@@ -145,5 +145,5 @@
 
 
-  if (wireframe)
+  if (mWireframe)
     stream<<"<IndexedLineSet ccw=\"TRUE\" coordIndex=\""<<endl;
   else
@@ -174,5 +174,5 @@
   stream<<"</Coordinate>"<<endl;
 
-  if (wireframe)
+  if (mWireframe)
     stream<<"</IndexedLineSet>"<<endl;
   else
@@ -211,4 +211,8 @@
 X3dExporter::ExportKdTree(const KdTree &tree)
 {
+  if (mExportRayDensity) {
+    return ExportKdTreeRayDensity(tree);
+  }
+  
   stack<KdNode *> tStack;
 
@@ -249,2 +253,59 @@
 
 
+bool
+X3dExporter::ExportKdTreeRayDensity(const KdTree &tree)
+{
+  stack<KdNode *> tStack;
+
+  tStack.push(tree.GetRoot());
+
+  bool fm = mUseForcedMaterial;
+  mUseForcedMaterial = true;
+  mForcedMaterial.mDiffuseColor.g = 1.0f;
+  mForcedMaterial.mDiffuseColor.b = 1.0f;
+  while (!tStack.empty()) {
+    KdNode *node = tStack.top();
+    tStack.pop();
+    if (node->IsLeaf()) {
+      AxisAlignedBox3 box = tree.GetBox(node);
+      Mesh *mesh = new Mesh;
+      
+      // add 6 vertices of the box
+      int index = mesh->mVertices.size();
+      for (int i=0; i < 8; i++) {
+	Vector3 v;
+	box.GetVertex(i, v);
+	mesh->mVertices.push_back(v);
+      }
+      mesh->AddFace(new Face(index + 0, index + 1, index + 3, index + 2) );
+      mesh->AddFace(new Face(index + 0, index + 2, index + 6, index + 4) );
+      mesh->AddFace(new Face(index + 4, index + 6, index + 7, index + 5) );
+      
+      mesh->AddFace(new Face(index + 3, index + 1, index + 5, index + 7) );
+      mesh->AddFace(new Face(index + 0, index + 4, index + 5, index + 1) );
+      mesh->AddFace(new Face(index + 2, index + 3, index + 7, index + 6) );
+
+
+      // set the mesh material according to the ray density
+      KdLeaf *leaf = (KdLeaf *) node;
+      if (leaf->mPassingRays.mRays) {
+	float importance = leaf->mPassingRays.mContributions/(float)leaf->mPassingRays.mRays;
+	//	float importance = leaf->mPassingRays.mContributions/1000.0f;
+	//	float importance = leaf->mPassingRays.mRays/1000.0f;
+	///(float)leaf->mPassingRays.mRays;
+	// mForcedMaterial.mDiffuseColor.r = log10(leaf->mPassingRays.mRays)/3.0f;
+	mForcedMaterial.mDiffuseColor.r = importance;
+	mForcedMaterial.mDiffuseColor.g = 1.0f - mForcedMaterial.mDiffuseColor.r;
+	ExportMesh(mesh);
+      }
+      delete mesh;
+    } else {
+      KdInterior *interior = (KdInterior *)node;
+      tStack.push(interior->mFront);
+      tStack.push(interior->mBack);
+    }
+  }
+  // restore the state of forced material
+  mUseForcedMaterial = fm;
+  return true;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h	(revision 177)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.h	(revision 191)
@@ -2,5 +2,5 @@
 #define _AxisAlignedBox3_H__
 
-
+#include "Rectangle3.h"
 #include "Matrix4x4.h"
 #include "Vector3.h"
@@ -232,9 +232,21 @@
   // x,y,z are either 0 or 1; (0 .. lower coordinate, 1 .. large coordinate)
   // (xmin,ymin, zmin) .. N = 0, (xmax, ymax, zmax) .. N= 7
-  void GetVertex(int N, Vector3 &vertex) const;
-  
+  void GetVertex(const int N, Vector3 &vertex) const;
+
+  Vector3 GetVertex(const int N) const {
+    Vector3 v;
+    GetVertex(N, v);
+    return v;
+  }
+
   // Returns 1, if the box includes on arbitrary face a given box
   int IsPiercedByBox(const AxisAlignedBox3 &box, int &axis) const;
 
+
+  int GetFaceVisibilityMask(const Vector3 &position) const;
+  int GetFaceVisibilityMask(const Rectangle3 &rectangle) const;
+
+  Rectangle3 GetFace(const int face) const;
+  
   // For a given point returns the region, where the point is located
   // there are 27 regions (0..26) .. determined by the planes embedding in the
@@ -243,5 +255,5 @@
   // R = 9*x + 3*y + z  ; e.g. region .. inside the box is 13.
   int GetRegionID(const Vector3 &point) const;
-
+  
   // Set the corner point of rectangle on the face of bounding box
   // given by the index number and the rectangle lying on this face
Index: trunk/VUT/GtpVisibilityPreprocessor/src/default.env
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/default.env	(revision 177)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/default.env	(revision 191)
@@ -8,6 +8,7 @@
 #	filename glasgow1.x3d
 #	filename vienna.x3d
-	filename atlanta2.x3d
+#	filename atlanta2.x3d
 #	filename soda.dat
+	filename soda5.dat
 
 }
@@ -53,5 +54,5 @@
 
 Sampling {
-	totalSamples	500000
-	samplesPerPass  2
+	totalSamples	1000000
+	samplesPerPass  5
 }
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Rectangle3.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Rectangle3.h	(revision 191)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Rectangle3.h	(revision 191)
@@ -0,0 +1,34 @@
+#ifndef __RECTANGLE3_H
+#define __RECTANGLE3_H
+
+#include "Vector3.h"
+
+class Rectangle3 {
+public:
+  Vector3 mVertices[4];
+
+  Rectangle3(const Vector3 &v0,
+	     const Vector3 &v1,
+	     const Vector3 &v2,
+	     const Vector3 &v3) {
+    mVertices[0] = v0;
+    mVertices[1] = v1;
+    mVertices[2] = v2;
+    mVertices[3] = v3;
+  }
+  
+  Vector3 GetNormal() {
+    return Normalize(CrossProd(mVertices[0]-mVertices[1],
+			       mVertices[2]-mVertices[1]
+			       ));
+  }
+  
+  Vector3 GetCenter() {
+    return (mVertices[0] + mVertices[1] + mVertices[2] + mVertices[3])/4.0f;
+  }
+
+  
+};
+
+
+#endif
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Makefile
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Makefile	(revision 181)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Makefile	(revision 191)
@@ -1,5 +1,5 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (1.07a) (Qt 3.3.2) on: Mon Jul 18 20:55:14 2005
+# Generated by qmake (1.07a) (Qt 3.3.2) on: Tue Aug 02 22:47:09 2005
 # Project:  preprocessor.pro
 # Template: app
@@ -13,12 +13,12 @@
 LEX		= flex
 YACC		= byacc
-CFLAGS	=	-nologo -Zm200 -W0 -MDd -Zi -GX  -DUNICODE -DQT_DLL -DQT_THREAD_SUPPORT
-CXXFLAGS	=	-nologo -Zm200 -W0 -MDd -Zi -GX  -DUNICODE -DQT_DLL -DQT_THREAD_SUPPORT
+CFLAGS	=	-nologo -Zm200 -W0 -MD -O2 -GX  -DUNICODE -DQT_DLL -DQT_THREAD_SUPPORT -DQT_NO_DEBUG
+CXXFLAGS	=	-nologo -Zm200 -W0 -MD -O2 -GX  -DUNICODE -DQT_DLL -DQT_THREAD_SUPPORT -DQT_NO_DEBUG
 LEXFLAGS	=
 YACCFLAGS	=-d
-INCPATH	=	 -I"..\src" -I"..\support\xerces\include" -I"..\support\zlib\include" -I"..\support\boost" -I"$(QTDIR)\include" -I"d:\gametools\svn\trunk\VUT\GtpVisibilityPreprocessor\src" -I"C:\Qt\3.3.2\mkspecs\win32-msvc"
+INCPATH	=	 -I"..\src" -I"..\support\xerces\include" -I"..\support\zlib\include" -I"..\support\boost" -I"..\support\devil\include" -I"$(QTDIR)\include" -I"d:\gametools\svn\trunk\VUT\GtpVisibilityPreprocessor\src" -I"C:\Qt\3.3.2\mkspecs\win32-msvc"
 LINK	=	link
-LFLAGS	=	/NOLOGO /DEBUG /SUBSYSTEM:console /LIBPATH:"../support/xerces/lib" /LIBPATH:"$(QTDIR)\lib"
-LIBS	=	 "qt-mt332.lib" "qtmain.lib" "kernel32.lib" "user32.lib" "gdi32.lib" "comdlg32.lib" "advapi32.lib" "shell32.lib" "ole32.lib" "oleaut32.lib" "uuid.lib" "imm32.lib" "winmm.lib" "wsock32.lib" "winspool.lib" "xerces-c_2.lib" "opengl32.lib" "glu32.lib" "delayimp.lib"
+LFLAGS	=	/NOLOGO delayimp.lib /DELAYLOAD:comdlg32.dll /DELAYLOAD:oleaut32.dll /DELAYLOAD:winmm.dll /DELAYLOAD:wsock32.dll /DELAYLOAD:winspool.dll /SUBSYSTEM:console /LIBPATH:"../support/xerces/lib" /LIBPATH:"../support/devil/lib" /LIBPATH:"$(QTDIR)\lib"
+LIBS	=	 "qt-mt332.lib" "qtmain.lib" "kernel32.lib" "user32.lib" "gdi32.lib" "comdlg32.lib" "advapi32.lib" "shell32.lib" "ole32.lib" "oleaut32.lib" "uuid.lib" "imm32.lib" "winmm.lib" "wsock32.lib" "winspool.lib" "xerces-c_2.lib" "devil.lib" "ilu.lib" "ilut.lib" "opengl32.lib" "glu32.lib" "delayimp.lib"
 MOC		=	$(QTDIR)\bin\moc.exe
 UIC		=	$(QTDIR)\bin\uic.exe
@@ -61,5 +61,7 @@
 		X3dParser.cpp \
 		MeshKdTree.cpp \
-		Pvs.cpp
+		Pvs.cpp \
+		MutualVisibility.cpp \
+		Triangle3.cpp
 OBJECTS =	Preprocessor.obj \
 		SamplingPreprocessor.obj \
@@ -83,5 +85,7 @@
 		X3dParser.obj \
 		MeshKdTree.obj \
-		Pvs.obj
+		Pvs.obj \
+		MutualVisibility.obj \
+		Triangle3.obj
 FORMS =	
 UICDECLS =	
@@ -167,9 +171,6 @@
 	-$(DEL_FILE) MeshKdTree.obj
 	-$(DEL_FILE) Pvs.obj
-	-$(DEL_FILE) preprocessor.pdb
-	-$(DEL_FILE) preprocessor.ilk
-	-$(DEL_FILE) vc*.pdb
-	-$(DEL_FILE) vc*.idb
-
+	-$(DEL_FILE) MutualVisibility.obj
+	-$(DEL_FILE) Triangle3.obj
 
 
@@ -190,4 +191,5 @@
 		Containers.h \
 		AxisAlignedBox3.h \
+		Rectangle3.h \
 		Matrix4x4.h \
 		Vector3.h \
@@ -210,6 +212,8 @@
 		X3dExporter.h \
 		Environment.h \
-		Containers.h \
-		AxisAlignedBox3.h \
+		MutualVisibility.h \
+		Containers.h \
+		AxisAlignedBox3.h \
+		Rectangle3.h \
 		Matrix4x4.h \
 		Vector3.h \
@@ -239,4 +243,5 @@
 		Material.h \
 		Pvs.h \
+		Rectangle3.h \
 		Vector3.h \
 		common.h \
@@ -255,4 +260,5 @@
 		Material.h \
 		Pvs.h \
+		Rectangle3.h \
 		Vector3.h \
 		Containers.h \
@@ -288,4 +294,5 @@
 		AxisAlignedBox3.h \
 		Pvs.h \
+		Rectangle3.h \
 		Containers.h \
 		Parser.h \
@@ -300,4 +307,5 @@
 		Containers.h \
 		AxisAlignedBox3.h \
+		Rectangle3.h \
 		Matrix4x4.h \
 		Vector3.h \
@@ -315,4 +323,5 @@
 		Containers.h \
 		AxisAlignedBox3.h \
+		Rectangle3.h \
 		Matrix4x4.h \
 		Vector3.h \
@@ -343,4 +352,5 @@
 		AxisAlignedBox3.h \
 		Ray.h \
+		Rectangle3.h \
 		Matrix4x4.h \
 		Vector3.h \
@@ -376,4 +386,5 @@
 		Material.h \
 		Pvs.h \
+		Rectangle3.h \
 		Vector3.h \
 		common.h \
@@ -393,4 +404,5 @@
 		Material.h \
 		Pvs.h \
+		Rectangle3.h \
 		Containers.h \
 		
@@ -411,4 +423,5 @@
 		AxisAlignedBox3.h \
 		common.h \
+		Rectangle3.h \
 		Matrix4x4.h \
 		Plane3.h \
@@ -431,4 +444,5 @@
 		Material.h \
 		Pvs.h \
+		Rectangle3.h \
 		Vector3.h \
 		common.h \
@@ -447,4 +461,5 @@
 		Material.h \
 		Pvs.h \
+		Rectangle3.h \
 		Vector3.h \
 		Containers.h \
@@ -456,4 +471,23 @@
 		
 
+MutualVisibility.obj: MutualVisibility.cpp  \
+		KdTree.h \
+		AxisAlignedBox3.h \
+		Ray.h \
+		MutualVisibility.h \
+		Containers.h \
+		Rectangle3.h \
+		Matrix4x4.h \
+		Vector3.h \
+		Plane3.h \
+		common.h \
+		
+
+Triangle3.obj: Triangle3.cpp  \
+		Triangle3.h \
+		Vector3.h \
+		common.h \
+		
+
 ####### Install
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Triangle3.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Triangle3.h	(revision 191)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Triangle3.h	(revision 191)
@@ -0,0 +1,43 @@
+#ifndef __TRIANGLE3_H
+#define __TRIANGLE3_H
+
+#include "Vector3.h"
+
+struct Triangle3 {
+  
+  Vector3 mVertices[3];
+  
+  Triangle3() {}
+  Triangle3(const Vector3 &a,
+		 const Vector3 &b,
+		 const Vector3 &c
+		 ) {
+    Init(a, b, c);
+  }
+
+  void
+  Init(const Vector3 &a,
+       const Vector3 &b,
+       const Vector3 &c
+       ) {
+    mVertices[0] = a;
+    mVertices[1] = b;
+    mVertices[2] = c;
+    
+  }
+
+  Vector3 GetNormal() const {
+    Vector3 v1=mVertices[0]-mVertices[1], v2=mVertices[2]-mVertices[1];
+    return Normalize(CrossProd(v2,v1));
+  }
+
+  Vector3 GetCenter() const {
+    return (mVertices[0] + mVertices[1] + mVertices[2])/3.0f;
+  }
+
+  float GetSpatialAngle(const Vector3 &point) const;
+};
+
+
+#endif
+
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Ray.h	(revision 191)
@@ -25,5 +25,5 @@
   enum { NO_INTERSECTION=0, INTERSECTION_OUT_OF_LIMITS, INTERSECTION };
 
-  struct RayIntersection {
+  struct Intersection {
     // the point of intersection
     float mT;
@@ -35,12 +35,12 @@
     int mFace;
 
-    RayIntersection(const float t,
-		    Intersectable *object,
-		    const int face):mT(t), mObject(object), mFace(face) {}
-
-    RayIntersection() {}
+    Intersection(const float t,
+		 Intersectable *object,
+		 const int face):mT(t), mObject(object), mFace(face) {}
+
+    Intersection() {}
 
     bool operator<(
-		   const RayIntersection &b) const {
+		   const Intersection &b) const {
       return 
 	mT
@@ -54,7 +54,7 @@
   // corresponds to the spatial elementary cell
   /** intersection with the source object if any */
-  RayIntersection sourceObject;
-  
-  vector<RayIntersection> intersections;
+  Intersection sourceObject;
+  
+  vector<Intersection> intersections;
   vector<KdLeaf *> leaves;
   vector<MeshInstance *> meshes;
@@ -74,4 +74,11 @@
   Ray() {}
 
+  Intersectable *GetIntersectionObject(const int i) const {
+    return intersections[i].mObject;
+  }
+  
+  Vector3 GetIntersectionPoint(const int i) const {
+    return Extrap(intersections[i].mT);
+  }
   
   // Inititalize the ray again when already constructed
@@ -151,6 +158,6 @@
 
   // the object on which the ray starts at
-  const RayIntersection* GetStartObject() const { return &intersections[0]; }
-  const RayIntersection* GetStopObject() const { return &intersections[intersections.size()-1]; }
+  const Intersection* GetStartObject() const { return &intersections[0]; }
+  const Intersection* GetStopObject() const { return &intersections[intersections.size()-1]; }
   
   
@@ -220,4 +227,22 @@
 
 
+class PassingRaySet {
+public:
+  enum { Resolution = 2 };
+  int mDirectionalContributions[3*Resolution*Resolution];
+  int mRays;
+  int mContributions;
+  PassingRaySet() {
+    Reset();
+  }
+  void
+  Reset();
+  
+  void AddRay(const Ray &ray, const int contributions);
+  int GetEntryIndex(const Vector3 &direction) const;
+
+  friend ostream &operator<<(ostream &s, const PassingRaySet &set);
+
+};
 
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.h	(revision 191)
@@ -10,5 +10,5 @@
 #include "Material.h"
 
-
+class Triangle3;
 class MeshInstance;
 class MeshKdTree;
@@ -73,4 +73,6 @@
       delete mFaces[i];
   }
+
+  void AddTriangle(const Triangle3 &triangle);
   
   void AddFace(Face *face) 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro	(revision 181)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/preprocessor.pro	(revision 191)
@@ -6,12 +6,12 @@
 #win32:INCLUDEPATH += c:/STLport-4.6.2/stlport
 
-INCLUDEPATH += ../src ../support/xerces/include  ../support/zlib/include ../support/boost
+INCLUDEPATH += ../src ../support/xerces/include  ../support/zlib/include ../support/boost ../support/devil/include
 
-win32:LIBPATH += ../support/xerces/lib
+win32:LIBPATH += ../support/xerces/lib ../support/devil/lib
 
 #win32:LIBPATH += c:/STLport-4.6.2/lib
 
 # debuc config
-CONFIG		+= console warn_off thread debug
+CONFIG		+= console warn_off thread release
 
 # RELEASE CONFIG
@@ -20,5 +20,5 @@
 # DEPENDPATH	= ../../include
 
-win32:LIBS += xerces-c_2.lib
+win32:LIBS += xerces-c_2.lib devil.lib ilu.lib ilut.lib
 
 # Input
@@ -31,5 +31,6 @@
 UnigraphicsParser.cpp X3dExporter.cpp SceneGraph.cpp Material.cpp \
 Matrix4x4.cpp Vector3.cpp AxisAlignedBox3.cpp Ray.cpp main.cpp Mesh.cpp \
-Exporter.cpp Camera.cpp X3dParser.cpp MeshKdTree.cpp Pvs.cpp
+Exporter.cpp Camera.cpp X3dParser.cpp MeshKdTree.cpp Pvs.cpp \
+MutualVisibility.cpp Triangle3.cpp
 
 
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.cpp	(revision 177)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Pvs.cpp	(revision 191)
@@ -39,5 +39,4 @@
   KdPvsMap::iterator i = mEntries.begin();
   for (int k = 0; k!=index && i != mEntries.end(); i++, k++);
-
   node = (*i).first;
   data = (*i).second;
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Exporter.h	(revision 191)
@@ -18,13 +18,16 @@
 protected:
   string mFilename;
-  bool wireframe;
+  bool mWireframe;
   bool mUseForcedMaterial;
   Material mForcedMaterial;
 
+  bool mExportRayDensity;
+  
 public:
   
   Exporter(const string filename):mFilename(filename),
-				  wireframe(false),
-				  mUseForcedMaterial(false)
+				  mWireframe(false),
+				  mUseForcedMaterial(false),
+				  mExportRayDensity(false)
   {
   }
@@ -49,6 +52,8 @@
   ExportIntersectable(Intersectable *object) = 0;
 
-  void SetWireframe() { wireframe = true; }
-  void SetFilled() { wireframe = false; }
+  void SetExportRayDensity(const bool d) { mExportRayDensity = d; }
+  
+  void SetWireframe() { mWireframe = true; }
+  void SetFilled() { mWireframe = false; }
 
   void SetForcedMaterial(const Material &m) {
Index: trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.cpp	(revision 181)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.cpp	(revision 191)
@@ -710,2 +710,61 @@
   return neighbors.size();
 }
+
+void
+KdTree::CollectLeaves(vector<KdLeaf *> &leaves)
+{
+  stack<KdNode *> nodeStack;
+  nodeStack.push(mRoot);
+
+  while (!nodeStack.empty()) {
+    KdNode *node = nodeStack.top();
+    nodeStack.pop();
+    if (node->IsLeaf()) {
+      KdLeaf *leaf = (KdLeaf *)node;
+      leaves.push_back(leaf);
+    } else {
+      KdInterior *interior = (KdInterior *)node;
+      nodeStack.push(interior->mBack);
+      nodeStack.push(interior->mFront);
+    }
+  }
+}
+
+
+int
+KdTree::CollectLeafPvs()
+{
+  int totalPvsSize = 0;
+  stack<KdNode *> nodeStack;
+  
+  nodeStack.push(mRoot);
+  
+  while (!nodeStack.empty()) {
+    KdNode *node = nodeStack.top();
+    nodeStack.pop();
+    if (node->IsLeaf()) {
+      KdLeaf *leaf = (KdLeaf *)node;
+      for (int j=0; j < leaf->mObjects.size(); j++) {
+	Intersectable *object = leaf->mObjects[j];
+	if (!object->Mailed()) {
+	  object->Mail();
+	  // add this node to pvs of all nodes it can see
+	  KdPvsMap::iterator ni = object->mKdPvs.mEntries.begin();
+	  for (; ni != object->mKdPvs.mEntries.end(); ni++) {
+	    KdNode *node = (*ni).first;
+	    // $$ JB TEMPORARY solution -> should add object PVS or explictly computed
+	    // kd tree PVS
+	    if (leaf->mKdPvs.AddNodeSample(node))
+	      totalPvsSize++;
+	  }
+	}
+      }
+    } else {
+      KdInterior *interior = (KdInterior *)node;
+      nodeStack.push(interior->mFront);
+      nodeStack.push(interior->mBack);
+    }
+  }
+
+  return totalPvsSize;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/MutualVisibility.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/MutualVisibility.cpp	(revision 191)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/MutualVisibility.cpp	(revision 191)
@@ -0,0 +1,288 @@
+#include <stack>
+using namespace std;
+#include "KdTree.h"
+#include "AxisAlignedBox3.h"
+#include "Ray.h"
+#include "MutualVisibility.h"
+
+
+void
+TriangleSample::Init(
+		     const Triangle3 &source,
+		     const Triangle3 &target)
+{
+  mDepth = 0;
+  mSource = source;
+  mTarget = target;
+}
+
+
+Vector3
+TriangleSample::GetIntersectionPoint(const int rayIndex,
+				     const int depth) const
+{
+  Vector3 origin, direction;
+  GetRay(rayIndex, origin, direction);
+  return origin + direction*mIntersections[rayIndex][depth].mT;
+}
+
+void
+TriangleSample::GetRay(const int rayIndex,
+		       Vector3 &origin,
+		       Vector3 &direction) const
+{
+  const int indices[][2] = {
+    {0,0}, {0,1}, {0,2},
+    {1,0}, {1,1}, {1,2},
+    {2,0}, {2,1}, {2,2}
+  };
+  
+  origin = mSource.mVertices[indices[rayIndex][0]];
+  direction = mTarget.mVertices[indices[rayIndex][1]] - origin;
+}
+
+void
+MutualVisibilitySampler::PerformSplit(
+				      const TriangleSample &sample,
+				      const bool splitSource,
+				      const int edge,
+				      const Vector3 &splitPoint,
+				      const Ray &ray,
+				      TriangleSample &sample1,
+				      TriangleSample &sample2
+				      )
+{
+
+
+  const Triangle3 *triangle;
+  Triangle3 *triangle1, *triangle2;
+  // split the triangles
+  if (splitSource) {
+    sample1.mTarget = sample.mTarget;
+    sample2.mTarget = sample.mTarget;
+    triangle1 = &sample1.mSource;
+    triangle2 = &sample2.mSource;
+    triangle = &sample.mSource;
+  } else {
+    sample1.mSource = sample.mSource;
+    sample2.mSource = sample.mSource;
+    triangle1 = &sample1.mTarget;
+    triangle2 = &sample2.mTarget;
+    triangle = &sample.mTarget;
+  }
+  
+  // split the intersections
+  switch (edge) {
+  case 0:
+    sample1.mIntersections[0] = sample.mIntersections[0];
+    sample1.mIntersections[1] = sample.mIntersections[1];
+    triangle1->mVertices[0] = triangle->mVertices[0];
+    triangle1->mVertices[1] = triangle->mVertices[1];
+    
+    sample2.mIntersections[0] = sample.mIntersections[2];
+    sample2.mIntersections[1] = sample.mIntersections[0];
+    triangle2->mVertices[0] = triangle->mVertices[2];
+    triangle2->mVertices[1] = triangle->mVertices[0];
+    break;
+
+  case 1:
+    sample1.mIntersections[0] = sample.mIntersections[0];
+    sample1.mIntersections[1] = sample.mIntersections[1];
+    triangle1->mVertices[0] = triangle->mVertices[0];
+    triangle1->mVertices[1] = triangle->mVertices[1];
+
+    sample2.mIntersections[0] = sample.mIntersections[1];
+    sample2.mIntersections[1] = sample.mIntersections[2];
+    triangle2->mVertices[0] = triangle->mVertices[1];
+    triangle2->mVertices[1] = triangle->mVertices[2];
+    break;
+    
+  case 2:
+    sample1.mIntersections[0] = sample.mIntersections[2];
+    sample1.mIntersections[1] = sample.mIntersections[0];
+    triangle1->mVertices[0] = triangle->mVertices[2];
+    triangle1->mVertices[1] = triangle->mVertices[0];
+
+    sample2.mIntersections[0] = sample.mIntersections[1];
+    sample2.mIntersections[1] = sample.mIntersections[2];
+    triangle2->mVertices[0] = triangle->mVertices[1];
+    triangle2->mVertices[1] = triangle->mVertices[2];
+    break;
+  }
+  
+  // the new sample
+  sample1.mIntersections[2] = ray.intersections;
+  sample2.mIntersections[2] = ray.intersections;
+  triangle1->mVertices[2] = splitPoint;
+  triangle2->mVertices[2] = splitPoint;
+  
+  // establish the new samples
+  sample1.mDepth = sample2.mDepth = sample.mDepth+1;
+}
+
+float
+MutualVisibilitySampler::GetSpatialAngle(const TriangleSample &sample,
+					 const Vector3 &point
+					 )
+{
+  const int sampleIndices[][3]={
+    (0,0,0),
+    (0,1,0),
+    (-1,-1,-1)
+  };
+  
+  float sum = 0.0f;
+  int i;
+  
+  for (i=0; i < 10; i++) {
+    Triangle3 triangle(sample.GetIntersectionPoint(sampleIndices[i][0], 0),
+			    sample.GetIntersectionPoint(sampleIndices[i][1], 0),
+			    sample.GetIntersectionPoint(sampleIndices[i][2], 0));
+    
+    sum += triangle.GetSpatialAngle(point);
+  }
+  return sum;
+}
+
+
+void
+MutualVisibilitySampler::ComputeError(TriangleSample &sample)
+{
+  // evaluate minimal error which can be achieved by more precise evaluation
+  // if this is above the threshold do not proceed further
+
+  float maxAngle = GetSpatialAngle(sample,
+				   sample.mSource.GetCenter());
+  
+  for (int i=0; i < 3; i++) {
+    float angle = GetSpatialAngle(sample,
+				  sample.mSource.mVertices[i]);
+    if (angle > maxAngle)
+      maxAngle = angle;
+  }
+  
+  sample.mError = maxAngle;
+}
+
+
+
+
+
+void
+MutualVisibilitySampler::ConstructInitialSamples(
+						 const AxisAlignedBox3 &source,
+						 const AxisAlignedBox3 &target,
+						 vector<TriangleSample *> &samples
+						 )
+{
+  // get all rectangles potentially visible from the source box
+  int i;
+  int sourceMask = 0;
+  for (i=0; i < 8; i++)
+    sourceMask |= source.GetFaceVisibilityMask(target.GetVertex(i));
+
+  // now for each visble source face find all visible target faces
+  for (i=0; i < 6; i++) {
+    Rectangle3 sourceFace = source.GetFace(i);
+    if ( sourceMask &(1<<i) ) {
+      int mask = target.GetFaceVisibilityMask(sourceFace);
+      // construct triangle samples for all visible rectangles
+      for (int j=0; j < 6; j++)
+	if (mask & (1<<j)) {
+	  AddInitialSamples(sourceFace, target.GetFace(j), samples);
+	}
+    }
+  }
+}
+
+void
+MutualVisibilitySampler::AddInitialSamples(
+					   const Rectangle3 &sourceRect,
+					   const Rectangle3 &targetRect,
+					   vector<TriangleSample *> &samples
+					   )
+{
+  Triangle3 sourceTriangles[2];
+  Triangle3 targetTriangles[2];
+
+  sourceTriangles[0].Init(sourceRect.mVertices[0],
+			  sourceRect.mVertices[1],
+			  sourceRect.mVertices[2]);
+
+  sourceTriangles[1].Init(sourceRect.mVertices[0],
+			  sourceRect.mVertices[2],
+			  sourceRect.mVertices[3]);
+
+  targetTriangles[0].Init(targetRect.mVertices[0],
+			  targetRect.mVertices[1],
+			  targetRect.mVertices[2]);
+
+  targetTriangles[1].Init(targetRect.mVertices[0],
+			  targetRect.mVertices[2],
+			  targetRect.mVertices[3]);
+  int i, j;
+  for (i=0; i < 2; i++)
+    for (j=0; j < 2; j++) {
+      TriangleSample *sample = new TriangleSample(sourceTriangles[i], 
+						  targetTriangles[j]);
+      samples.push_back(sample);
+    }
+}
+
+int
+MutualVisibilitySampler::ComputeVisibility()
+{
+  stack<TriangleSample *> sampleStack;
+  
+  Ray ray;
+  // now process the samples as long as we have something to do
+  while (!sampleStack.empty()) {
+    TriangleSample *sample = sampleStack.top();
+    sampleStack.pop();
+    
+//      // cast a new ray
+//      int triangleSplitEdge = SetupExtremalRay(sample, source, ray);
+    
+//      if (!kdTree->CastRay(ray))
+//        return VISIBLE;
+    
+//      // generate 2 new samples
+//      TriangleSample newSamples[2];
+//      sample.Split(triangleSplitEdge, ray, newSamples[0], newSamples[1]);
+//      for (i=0; i < 2; i++) {
+//        newSamples[i].ComputeError();
+//        if (newSamples[i].mError > solidAngleThreshold) {
+//  	sampleStack.push(newSamples[i]);
+//        }
+//      }
+//    }
+  }
+  return INVISIBLE;
+}
+
+MutualVisibilitySampler::MutualVisibilitySampler(KdTree *kdTree,
+						 AxisAlignedBox3 &source,
+						 AxisAlignedBox3 &target,
+						 const float solidAngleThreshold)
+{
+  mKdTree = kdTree;
+  mSource = source;
+  mTarget = target;
+  mSolidAngleThreshold = solidAngleThreshold;
+}
+  
+
+int
+ComputeBoxVisibility(KdTree *kdTree,
+		     AxisAlignedBox3 &source,
+		     AxisAlignedBox3 &target,
+		     const float solidAngleThreshold)
+{
+  MutualVisibilitySampler sampler(kdTree, source, target, solidAngleThreshold);
+
+  int visibility = sampler.ComputeVisibility();
+
+  return visibility;
+  
+}
+  
Index: trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 181)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/SamplingPreprocessor.cpp	(revision 191)
@@ -4,4 +4,6 @@
 #include "X3dExporter.h"
 #include "Environment.h"
+#include "MutualVisibility.h"
+
 SamplingPreprocessor::SamplingPreprocessor()
 {
@@ -33,14 +35,35 @@
 
 int
-SamplingPreprocessor::AddNodeSamples(Intersectable *object, const Ray &ray)
+SamplingPreprocessor::AddNodeSamples(Intersectable *object,
+				     const Ray &ray,
+				     const int pass)
 {
   int contributingSamples = 0;
-
-  for (int j=0; j < ray.leaves.size(); j++) {
+  int j;
+  for (j=0; j < ray.leaves.size(); j++) {
     KdNode *node = GetNodeForPvs( ray.leaves[j] );
     contributingSamples += object->mKdPvs.AddNodeSample(node);
   }
   
+  if (pass > 10)
+    for (j=1; j < ray.leaves.size() - 1; j++) {
+      ray.leaves[j]->AddPassingRay(ray, contributingSamples ? 1 : 0);
+    }
+  
   return contributingSamples;
+}
+
+
+void
+SamplingPreprocessor::HoleSamplingPass()
+{
+  vector<KdLeaf *> leaves;
+  mKdTree->CollectLeaves(leaves);
+  
+  // go through all the leaves and evaluate their passing contribution
+  for (int i=0 ; i < leaves.size(); i++) {
+    KdLeaf *leaf = leaves[i];
+    cout<<leaf->mPassingRays<<endl;
+  }
 }
 
@@ -63,5 +86,4 @@
   int totalSamples = 0;
 
-
   int pvsOut = Min((int)objects.size(), 10);
   
@@ -73,108 +95,118 @@
     int passSampleContributions = 0;
     int passSamples = 0;
-	int index = 0;
-    KdNode *visibleNode = NULL;
+    int index = 0;
 
     for (i =0; i < objects.size(); i++) {
       KdNode *nodeToSample = NULL;
       Intersectable *object = objects[i];
-
+      
       int pvsSize = object->mKdPvs.GetSize();
-
-      
-    if (pvsSize) {
-		// mail all nodes from the pvs
-		Intersectable::NewMail();
-		KdPvsMap::iterator i = object->mKdPvs.mEntries.begin();
-		for (; i != object->mKdPvs.mEntries.end(); i++) {
-			KdNode *node = (*i).first;
-			node->Mail();
-		}
-		int maxTries = 2*pvsSize;
-	
-		for (int tries = 0; tries < 10; tries++) {
-			index = RandomValue(0, pvsSize - 1);
-			KdPvsData data;
-			object->mKdPvs.GetData(index, visibleNode, data);
-			nodeToSample = mKdTree->FindRandomNeighbor(visibleNode, true);
-			if (nodeToSample)
-			break;
-		}
-	}
-
-    if (0 && pvsSize) {
-		// mail all nodes from the pvs
-		Intersectable::NewMail();
-		KdPvsMap::iterator i = object->mKdPvs.mEntries.begin();
-		for (; i != object->mKdPvs.mEntries.end(); i++) {
-			KdNode *node = (*i).first;
-			node->Mail();
-		}
-	
-		vector<KdNode *> invisibleNeighbors;
-		// get all neighbors of all PVS nodes
-		i = object->mKdPvs.mEntries.begin();
-		for (; i != object->mKdPvs.mEntries.end(); i++) {
-			KdNode *node = (*i).first;
-			mKdTree->FindNeighbors(visibleNode, invisibleNeighbors, true);
-		}
-
-		KdPvsData data;
-		KdNode *visibleNode;
-		object->mKdPvs.GetData(index, visibleNode, data);
-		nodeToSample = mKdTree->FindRandomNeighbor(visibleNode, true);
-		if (nodeToSample)
-			break;
-      }
-      
+      
+      if (0 && pvsSize) {
+	// mail all nodes from the pvs
+	Intersectable::NewMail();
+	KdPvsMap::iterator i = object->mKdPvs.mEntries.begin();
+	for (; i != object->mKdPvs.mEntries.end(); i++) {
+	  KdNode *node = (*i).first;
+	  node->Mail();
+	}
+	int maxTries = 2*pvsSize;
+	
+	for (int tries = 0; tries < 10; tries++) {
+	  index = RandomValue(0, pvsSize - 1);
+	  KdPvsData data;
+	  KdNode *node;
+	  object->mKdPvs.GetData(index, node, data);
+	  nodeToSample = mKdTree->FindRandomNeighbor(node, true);
+	  if (nodeToSample)
+	    break;
+	}
+      }
+      
+      if (0 && pvsSize && pass == 100 ) {
+	// mail all nodes from the pvs
+	Intersectable::NewMail();
+	KdPvsMap::iterator i = object->mKdPvs.mEntries.begin();
+	for (; i != object->mKdPvs.mEntries.end(); i++) {
+	  KdNode *node = (*i).first;
+	  node->Mail();
+	}
+	
+	vector<KdNode *> invisibleNeighbors;
+	// get all neighbors of all PVS nodes
+	i = object->mKdPvs.mEntries.begin();
+	for (; i != object->mKdPvs.mEntries.end(); i++) {
+	  KdNode *node = (*i).first;
+	  mKdTree->FindNeighbors(node, invisibleNeighbors, true);
+	  AxisAlignedBox3 box = object->GetBox();
+	  for (int j=0; j < invisibleNeighbors.size(); j++) {
+	    int visibility = ComputeBoxVisibility(mKdTree,
+						  box,
+						  mKdTree->GetBox(invisibleNeighbors[j]),
+						  1.0f);
+	    
+	  }
+	  
+	
+	  // now rank all the neighbors according to probability that a new
+	  // sample creates some contribution
+	}
+      }
 
       for (int k=0; k < mSamplesPerPass; k++) {
-		  object->GetRandomSurfacePoint(point, normal);
-		  if (nodeToSample) {
-			  int maxTries = 5;
-			  for (int tries = 0; tries < maxTries; tries++) {
-				  direction = mKdTree->GetBox(nodeToSample).GetRandomPoint() - point;
-				  if (DotProd(direction, normal) > Limits::Small)
-					break;
-			  }
+	object->GetRandomSurfacePoint(point, normal);
+	if (nodeToSample) {
+	  int maxTries = 5;
+	  for (int tries = 0; tries < maxTries; tries++) {
+	    direction = mKdTree->GetBox(nodeToSample).GetRandomPoint() - point;
+	    if (DotProd(direction, normal) > Limits::Small)
+	      break;
+	  }
 		
-			  if (tries == maxTries)
-				  direction = UniformRandomVector(normal);
-		  } else
-			  direction = UniformRandomVector(normal);
-	
-		  // construct a ray
-		  SetupRay(ray, point, direction);
-		  mKdTree->CastRay(ray);
-
-		  if (i < pvsOut)
-			  rays[i].push_back(ray);
-	
-		  int sampleContributions = 0;
+	  if (tries == maxTries)
+	    direction = UniformRandomVector(normal);
+	} else
+	  direction = UniformRandomVector(normal);
+	
+	// construct a ray
+	SetupRay(ray, point, direction);
+	mKdTree->CastRay(ray);
+
+	if (i < pvsOut)
+	  rays[i].push_back(ray);
+	
+	int sampleContributions = 0;
 	
 		  
-		  if (ray.leaves.size()) {
-			  sampleContributions += AddNodeSamples(object, ray);
+	if (ray.leaves.size()) {
+	  sampleContributions += AddNodeSamples(object, ray, pass);
 	  
-			  if (ray.intersections.size()) {
-				  sampleContributions += AddNodeSamples(ray.intersections[0].mObject, ray);
-				  // check whether we can add this to the rays
-				  for (int j = 0; j < pvsOut; j++) {
-					  if (objects[j] == ray.intersections[0].mObject) {
-						  rays[j].push_back(ray);
-					  }
-				  }
-			  }
+	  if (ray.intersections.size()) {
+	    sampleContributions += AddNodeSamples(ray.intersections[0].mObject, ray, pass);
+	    // check whether we can add this to the rays
+	    for (int j = 0; j < pvsOut; j++) {
+	      if (objects[j] == ray.intersections[0].mObject) {
+		rays[j].push_back(ray);
+	      }
+	    }
+	  }
 		
-			  passSamples++;
+	  passSamples++;
 			
-			  if (sampleContributions) {
-				  passContributingSamples++;
-				  passSampleContributions += sampleContributions;
-			  }
-		  }
-	  }
-	}
+	  if (sampleContributions) {
+	    passContributingSamples++;
+	    passSampleContributions += sampleContributions;
+	  }
+	}
+      }
+    }
+
+
     totalSamples += passSamples;
+
+    //    if (pass>10)
+    //      HoleSamplingPass();
+    
+
     pass++;
     
@@ -193,4 +225,16 @@
       endl;
   }
+
+  int totalPvsSize = mKdTree->CollectLeafPvs();
+  cout<<"#totalPvsSize="<<totalPvsSize<<endl;
+  
+  //  HoleSamplingPass();
+  if (1) {
+    Exporter *exporter = Exporter::GetExporter("ray-density.x3d");
+    exporter->SetExportRayDensity(true);
+    exporter->ExportKdTree(*mKdTree);
+    delete exporter;
+  }
+  
   bool exportRays = false;
   if (exportRays) {
@@ -214,5 +258,5 @@
       KdPvsMap::iterator i = object->mKdPvs.mEntries.begin();
       Intersectable::NewMail();
-	// avoid adding the object to the list
+      // avoid adding the object to the list
       object->Mail();
       ObjectContainer visibleObjects;
Index: trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp	(revision 177)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp	(revision 191)
@@ -549,5 +549,5 @@
 // x,y,z are either 0 or 1; (0 .. min coordinate, 1 .. max coordinate)
 void
-AxisAlignedBox3::GetVertex(int N, Vector3 &vertex) const
+AxisAlignedBox3::GetVertex(const int N, Vector3 &vertex) const
 {
   switch (N) {
@@ -1526,2 +1526,91 @@
   return (m == M) ? m : m + M;
 }
+
+int
+AxisAlignedBox3::GetFaceVisibilityMask(const Rectangle3 &rectangle) const
+{
+  int mask = 0;
+  for (int i=0; i < 4; i++)
+    mask |= GetFaceVisibilityMask(rectangle.mVertices[i]);
+  return mask;
+}
+
+int
+AxisAlignedBox3::GetFaceVisibilityMask(const Vector3 &position) const {
+  
+  // assume that we are not inside the box
+  int c=0;
+  
+  if (position.x<(mMin.x-Limits::Small))
+    c|=1;
+  else
+    if (position.x>(mMax.x+Limits::Small))
+      c|=2;
+  
+  if (position.y<(mMin.y-Limits::Small))
+    c|=4;
+  else
+    if (position.y>(mMax.y+Limits::Small))
+      c|=8;
+  
+  if (position.z<(mMin.z-Limits::Small))
+    c|=16;
+  else
+    if (position.z>(mMax.z+Limits::Small))
+      c|=32;
+  
+  return c;
+}
+
+
+Rectangle3
+AxisAlignedBox3::GetFace(const int face) const
+{
+  Vector3 v[4];
+  switch (face) {
+    
+  case 0:
+    v[3].SetValue(mMin.x,mMin.y,mMin.z);
+    v[2].SetValue(mMin.x,mMax.y,mMin.z);
+    v[1].SetValue(mMin.x,mMax.y,mMax.z);
+    v[0].SetValue(mMin.x,mMin.y,mMax.z);
+    break;
+    
+  case 1: 
+    v[0].SetValue(mMax.x,mMin.y,mMin.z);
+    v[1].SetValue(mMax.x,mMax.y,mMin.z);
+    v[2].SetValue(mMax.x,mMax.y,mMax.z);
+    v[3].SetValue(mMax.x,mMin.y,mMax.z);
+    break;
+    
+  case 2:
+    v[3].SetValue(mMin.x,mMin.y,mMin.z);
+    v[2].SetValue(mMin.x,mMin.y,mMax.z);
+    v[1].SetValue(mMax.x,mMin.y,mMax.z);
+    v[0].SetValue(mMax.x,mMin.y,mMin.z);
+    break;
+  
+  case 3:
+    v[0].SetValue(mMin.x,mMax.y,mMin.z);
+    v[1].SetValue(mMin.x,mMax.y,mMax.z);
+    v[2].SetValue(mMax.x,mMax.y,mMax.z);
+    v[3].SetValue(mMax.x,mMax.y,mMin.z);
+    break;
+
+  case 4:
+    v[3].SetValue(mMin.x,mMin.y,mMin.z);
+    v[2].SetValue(mMax.x,mMin.y,mMin.z);
+    v[1].SetValue(mMax.x,mMax.y,mMin.z);
+    v[0].SetValue(mMin.x,mMax.y,mMin.z);
+    break;
+  
+  case 5:
+    v[0].SetValue(mMin.x,mMin.y,mMax.z);
+    v[1].SetValue(mMax.x,mMin.y,mMax.z);
+    v[2].SetValue(mMax.x,mMax.y,mMax.z);
+    v[3].SetValue(mMin.x,mMax.y,mMax.z);
+    break;
+  }
+  
+  return Rectangle3(v[0], v[1], v[2], v[3]);
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Camera.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Camera.cpp	(revision 182)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Camera.cpp	(revision 191)
@@ -1,5 +1,9 @@
 #include <algorithm>
-#include <qstring.h>
-#include <qimage.h>
+
+// the devil library
+#include <IL/il.h>
+#include <IL/ilu.h>
+#include <IL/ilut.h>
+
 #include "Camera.h"
 #include "Ray.h"
@@ -7,4 +11,23 @@
 #include "Mesh.h"
 #include "Exporter.h"
+
+
+
+
+void
+InitDevIl()
+{
+  ilInit();
+  ILuint ImageName;
+  ilGenImages(1, &ImageName);
+  ilBindImage(ImageName);
+  ilEnable(IL_FILE_OVERWRITE);
+
+  //  ilRegisterFormat(IL_RGBA);
+  //  ilRegisterType(IL_FLOAT);
+
+  //  ilEnable(IL_ORIGIN_SET);
+  //  ilOriginFunc(IL_ORIGIN_UPPER_LEFT);
+}
 
 bool
@@ -19,9 +42,9 @@
 
   
-#if QT_VERSION < 0x040000
-  QImage image(mWidth, mHeight, 32);
-#else
-  //QImage image(mWidth, mHeight, QImage::Format_RGB32);
-#endif
+  InitDevIl();
+  int components = 4;
+  float *buffer = new float[components*mWidth*mHeight];
+  float *pbuffer = buffer;
+  //  - components*mWidth;
   
   vector<Ray> rays;
@@ -30,8 +53,9 @@
 
   Ray ray;
+
   for (y = 0; y < mHeight; y++) {
     cout<<"+";
     for (x = 0; x < mWidth; x++) {
-      SetupRay(ray, x, y);
+      SetupRay(ray, x, mHeight - (y + 1));
       if (tree->CastRay(ray))
 	if (ray.intersections.size()) {
@@ -41,16 +65,17 @@
 	  if (mesh->GetMesh()->mMaterial)
 	    color = mesh->GetMesh()->mMaterial->mDiffuseColor;
-#if QT_VERSION < 0x040000
-	  image.setPixel(x, y, qRgb(color.r*255,
-				    color.g*255,
-				    color.b*255
-				    ));
-#endif
+
+	  pbuffer[0] = color.r;
+	  pbuffer[1] = color.g;
+	  pbuffer[2] = color.b;
+	  pbuffer[3] = 1.0f;
 
 	}
+      pbuffer+=components;
       
       if (exportRays && (x==222) && (y==97))
 	rays.push_back(ray);
     }
+    //    pbuffer-=2*components*mWidth;
   }
   
@@ -58,9 +83,13 @@
   cout<<"#RAY_CAST_TIME\n";
   cout<<TimeDiff(t1, t2)<<"\n";
+
   
   cout<<"Saving image"<<endl;
-#if QT_VERSION < 0x040000
-  image.save(filename.c_str(), "PNG");
-#endif
+
+  ilRegisterType(IL_FLOAT);
+  ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
+  ilSaveImage((char *const)filename.c_str());
+  delete buffer;
+  
   Exporter *exporter = NULL;
   if (exportRays) {
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Triangle3.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Triangle3.cpp	(revision 191)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Triangle3.cpp	(revision 191)
@@ -0,0 +1,7 @@
+#include "Triangle3.h"
+
+float
+Triangle3::GetSpatialAngle(const Vector3 &point) const
+{
+  return 0.0f;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 162)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Ray.cpp	(revision 191)
@@ -121,2 +121,48 @@
 }
 
+void
+PassingRaySet::Reset()
+{
+  for (int i=0; i < 3*Resolution*Resolution; i++)
+    mDirectionalContributions[i] = 0;
+  mRays = 0;
+  mContributions = 0;
+}
+  
+void
+PassingRaySet::AddRay(const Ray &ray, const int contributions)
+{
+  int i = GetEntryIndex(ray.GetDir());
+  mRays++;
+  mContributions += contributions;
+  mDirectionalContributions[i] += contributions;
+}
+
+int
+PassingRaySet::GetEntryIndex(const Vector3 &direction) const
+{
+  // get face
+  int axis = direction.DrivingAxis();
+  Vector3 dir;
+  float k = direction[axis];
+  if ( k < 0.0f)
+    k = -k;
+
+  dir = direction/k;
+  float x, y;
+  dir.ExtractVerts(&x, &y, axis);
+  int ix = (x + 1.0f)*0.5f*Resolution;
+  int iy = (y + 1.0f)*0.5f*Resolution;
+
+  return Resolution*(Resolution*axis + iy) + ix;
+}
+
+ostream &
+operator<<(ostream &s, const PassingRaySet &set)
+{
+  s<<"Ray Set #rays="<<set.mRays<<" #contributions="<<set.mContributions<<endl;
+  for (int i=0; i < 3*sqr(PassingRaySet::Resolution); i++)
+    s<<set.mDirectionalContributions[i]<<" ";
+  s<<endl;
+  return s;
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/Mesh.cpp	(revision 191)
@@ -2,4 +2,5 @@
 #include "Mesh.h"
 #include "MeshKdTree.h"
+#include "Triangle3.h"
 
 int MeshInstance::mailID = 21843194198;
@@ -60,5 +61,5 @@
     switch (ray.GetType()) {
     case Ray::GLOBAL_RAY:
-      ray.intersections.push_back(Ray::RayIntersection(t, instance, faceIndex));
+      ray.intersections.push_back(Ray::Intersection(t, instance, faceIndex));
       hit++;
       break;
@@ -101,7 +102,7 @@
   if ( hits && ray.GetType() == Ray::LOCAL_RAY ) {
     if (ray.intersections.size())
-      ray.intersections[0] = Ray::RayIntersection(nearestT, instance, nearestFace);
+      ray.intersections[0] = Ray::Intersection(nearestT, instance, nearestFace);
     else
-      ray.intersections.push_back(Ray::RayIntersection(nearestT, instance, nearestFace));
+      ray.intersections.push_back(Ray::Intersection(nearestT, instance, nearestFace));
   }
   
@@ -135,7 +136,7 @@
   if ( hits && ray.GetType() == Ray::LOCAL_RAY ) {
     if (ray.intersections.size())
-      ray.intersections[0] = Ray::RayIntersection(nearestT, instance, nearestFace);
+      ray.intersections[0] = Ray::Intersection(nearestT, instance, nearestFace);
     else
-      ray.intersections.push_back(Ray::RayIntersection(nearestT, instance, nearestFace));
+      ray.intersections.push_back(Ray::Intersection(nearestT, instance, nearestFace));
   }
   
@@ -333,2 +334,13 @@
 }
 
+void
+Mesh::AddTriangle(const Triangle3 &triangle)
+{
+  int index = mVertices.size();
+
+  for (int i=0; i < 3; i++) {
+    mVertices.push_back(triangle.mVertices[i]);
+  }
+  
+  AddFace(new Face(index + 0, index + 1, index + 2) );
+}
Index: trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h
===================================================================
--- trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h	(revision 176)
+++ trunk/VUT/GtpVisibilityPreprocessor/src/X3dExporter.h	(revision 191)
@@ -57,5 +57,7 @@
   virtual void
   ExportSceneNode(SceneGraphNode *node);
-  
+
+  bool
+  ExportKdTreeRayDensity(const KdTree &tree);
 
   
