source: GTP/trunk/Lib/Vis/Preprocessing/src/VssPreprocessor.cpp @ 1221

#include "SceneGraph.h"
#include "KdTree.h"
#include "VssPreprocessor.h"
#include "X3dExporter.h"
#include "Environment.h"
#include "MutualVisibility.h"
#include "Polygon3.h"
#include "ViewCell.h"
#include "VssRay.h"
#include "VssTree.h"
#include "ViewCellsManager.h"
#include "RenderSimulator.h"
#include "Beam.h"
#include "GlRenderer.h"
#include "ArchModeler2MLRT.hxx"
#include "Intersectable.h"


namespace GtpVisibilityPreprocessor {

bool use2dSampling = false;
bool useViewspacePlane = false;

VssPreprocessor::VssPreprocessor():
  mPass(0),
  mVssRays(),
  mViewSpaceBox(NULL)
{
  // this should increase coherence of the samples
  Environment::GetSingleton()->GetIntValue("VssPreprocessor.samplesPerPass", mSamplesPerPass);
  Environment::GetSingleton()->GetIntValue("VssPreprocessor.initialSamples", mInitialSamples);
  Environment::GetSingleton()->GetIntValue("VssPreprocessor.vssSamples", mVssSamples);
  Environment::GetSingleton()->GetIntValue("VssPreprocessor.vssSamplesPerPass", mVssSamplesPerPass);
  Environment::GetSingleton()->GetBoolValue("VssPreprocessor.useImportanceSampling", mUseImportanceSampling);
 
  Environment::GetSingleton()->GetBoolValue("VssPreprocessor.loadInitialSamples", mLoadInitialSamples);
  Environment::GetSingleton()->GetBoolValue("VssPreprocessor.storeInitialSamples", mStoreInitialSamples);
  Environment::GetSingleton()->GetBoolValue("VssPreprocessor.useViewSpaceBox", mUseViewSpaceBox);
  Environment::GetSingleton()->GetBoolValue("VssPreprocessor.testBeamSampling", mTestBeamSampling);
  Environment::GetSingleton()->GetBoolValue("VssPreprocessor.enlargeViewSpace", mEnlargeViewSpace);
  Environment::GetSingleton()->GetBoolValue("Preprocessor.detectEmptyViewSpace", mDetectEmptyViewSpace);
  
  useViewspacePlane = mUseViewSpaceBox; //hack

  
  Debug << "*********** vss preprocessor options **************" << endl;
  Debug << "use view space box=" << mUseViewSpaceBox << endl;
  Debug << "enlarge view space=" << mEnlargeViewSpace << endl;
  Debug << "*********** end vss preprocessor options **************" << endl;

  mStats.open("stats.log");
}


VssPreprocessor::~VssPreprocessor()
{
        CLEAR_CONTAINER(mVssRays);
        DEL_PTR(mViewSpaceBox);
}


void
VssPreprocessor::SetupRay(Ray &ray,
                                                  const Vector3 &point,
                                                  const Vector3 &direction
                                                  )
{
  ray.Clear();
  // do not store anything else then intersections at the ray
  ray.Init(point, direction, Ray::LOCAL_RAY);
}


void
VssPreprocessor::CastRays(
                                                  SimpleRayContainer &rays,
                                                  VssRayContainer &vssRays
                                                  )
{
        AxisAlignedBox3 &box =  mViewSpaceBox ? *mViewSpaceBox : mKdTree->GetBox();
        
#if 1
        for (int i=0; i < rays.size(); i++)
        {
                CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays, box);
        }
#else
        for (int i = 0; i < rays.size(); i += 16)
        {
                CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays, box);
        }

#endif
}


void VssPreprocessor::CastRays16(SimpleRayContainer &rays, 
                                                                 VssRayContainer &vssRays,
                                                                 const AxisAlignedBox3 &sbox)
{
/*
        mlrtaStoreRayAS16(const float* pforg, const float* pfdir, int j);


        // input parameters
        const RTVec3f& org = *reinterpret_cast<const RTVec3f*> (pforg);
        const RTVec3f& rd  = *reinterpret_cast<const RTVec3f*> (pfdir);

        int ri = j/4;
        int ci = j%4;
        *((float*)&ibp4.ray4X[ri].m_origin.t[0] + ci) = org[0];
        *((float*)&ibp4.ray4X[ri].m_origin.t[1] + ci) = org[1];
        *((float*)&ibp4.ray4X[ri].m_origin.t[2] + ci) = org[2];

        *((float*)&ibp4.ray4X[ri].m_direction.t[0] + ci) = rd[0];
        *((float*)&ibp4.ray4X[ri].m_direction.t[1] + ci) = rd[1];
        *((float*)&ibp4.ray4X[ri].m_direction.t[2] + ci) = rd[2];

        mlrtaTraverseGroupAS16(const float* bbmin, const float* bbmax, int* hit_triangles, float* dist);*/
}


int VssPreprocessor::CastIntelDoubleRay(
                                                                                const Vector3 &viewPoint,
                                                                                const Vector3 &direction,
                                                                                VssRayContainer &vssRays,
                                                                                const AxisAlignedBox3 &box)
{
        VssRay *vssRay  = NULL;
        int hits = 0;

        Vector3 pointA, pointB;
        
        Intersectable *objectA = 
                CastIntelSingleRay(viewPoint, direction, pointA, box);
        
        // cast ray into other direction
        Intersectable *objectB = 
                CastIntelSingleRay(viewPoint, -direction, pointB, box);

        const bool validSample = (objectA != objectB);
        
        if (validSample) 
        {       
                if (objectA) 
                {
                        vssRay = new VssRay(pointB,
                                                                pointA,
                                                                objectB,
                                                                objectA,
                                                                mPass);
                        vssRays.push_back(vssRay);
                        hits ++;
                }

                if (objectB) 
                {
                        vssRay = new VssRay(pointA,
                                                                pointB,
                                                                objectA,
                                                                objectB,
                                                                mPass);
                        vssRays.push_back(vssRay);
                        hits ++;
                }
                //Debug << "intel ray: " << *vssRay << endl;
        }
//cout << "a";
        return hits;
}


Intersectable *VssPreprocessor::CastIntelSingleRay(const Vector3 &viewPoint,
                                                                                                   const Vector3 &direction,
                                                                                                   Vector3 &tPoint,
                                                                                                   const AxisAlignedBox3 &box
                                                                                                   )
{
        AxisAlignedBox3 sbox = box;
        sbox.Enlarge(Vector3(-Limits::Small));

        if (!sbox.IsInside(viewPoint))
                return 0;
        
        float pforg[3];
        float pfdir[3];
        double pfnorm[3];

        pforg[0] = viewPoint[0]; pforg[1] = viewPoint[1]; pforg[2] = viewPoint[2];
        pfdir[0] = direction[0]; pfdir[1] = direction[1]; pfdir[2] = direction[2];

        float dist = 0;
        const int hittriangle = mlrtaIntersectAS(pforg, pfdir, pfnorm, dist);

        if (hittriangle == -1)
        {
                static Ray ray;
                SetupRay(ray, viewPoint, direction);

                float tmin = 0, tmax;
                if (box.ComputeMinMaxT(ray, &tmin, &tmax) && tmin < tmax)
                {
                        tPoint = ray.Extrap(tmax);
                }

                return NULL;
        }
        else
        {
                tPoint[0] = pforg[0] + pfdir[0] * dist;
                tPoint[1] = pforg[1] + pfdir[1] * dist;
                tPoint[2] = pforg[2] + pfdir[2] * dist;

                return mFaceParents[hittriangle];
        }
}


int
VssPreprocessor::CastInternalRay(
                                                                 const Vector3 &viewPoint,
                                                                 const Vector3 &direction,
                                                                 VssRayContainer &vssRays,
                                                                 const AxisAlignedBox3 &box 
                                                                 )
{
    int hits = 0;
        static Ray ray;

        AxisAlignedBox3 sbox = box;
        sbox.Enlarge(Vector3(-Limits::Small));

        if (!sbox.IsInside(viewPoint))
                return 0;
        
        SetupRay(ray, viewPoint, direction);
        
        // cast ray to KD tree to find intersection with other objects
        Intersectable *objectA, *objectB;
        
        Vector3 pointA;
        Vector3 pointB;

        //float bsize = Magnitude(box.Size());
        if (!mDetectEmptyViewSpace)
                ray.mFlags &= ~Ray::CULL_BACKFACES;
        else
                ray.mFlags |= Ray::CULL_BACKFACES;

        if (mKdTree->CastRay(ray)) 
        {
                objectA = ray.intersections[0].mObject;
                pointA = ray.Extrap(ray.intersections[0].mT);
        } 
        else 
        {
                objectA = NULL;
                // compute intersection with the scene bounding box
                float tmin, tmax;
                if (box.ComputeMinMaxT(ray, &tmin, &tmax) && tmin < tmax)
                        pointA = ray.Extrap(tmax);
                else
                        return 0;
        }

        // matt: point A could be undefined?

        // cast ray into opposite direction
        if (1 && mDetectEmptyViewSpace) {
                SetupRay(ray, pointA, -direction);
        } else
                SetupRay(ray, viewPoint, -direction);
  
        if (!mDetectEmptyViewSpace)
                ray.mFlags &= ~Ray::CULL_BACKFACES;
        else
                ray.mFlags |= Ray::CULL_BACKFACES;

        if (mKdTree->CastRay(ray)) 
        {
                objectB = ray.intersections[0].mObject;
                pointB = ray.Extrap(ray.intersections[0].mT);
        } 
        else 
        {
                objectB = NULL;
                float tmin, tmax;
                
                if (box.ComputeMinMaxT(ray, &tmin, &tmax) && tmin < tmax)
                        pointB = ray.Extrap(tmax);
                else
                        return 0;
        }
  
        //bool thesame = objectA && (objectA == objectB);
        
        //  if (objectA == NULL && objectB != NULL) {
        if (1 && mDetectEmptyViewSpace) 
        {
                // cast again to ensure that there is no objectA
                SetupRay(ray, pointB, direction);
                ray.mFlags |= Ray::CULL_BACKFACES;
                
                if (mKdTree->CastRay(ray)) 
                {
                        objectA = ray.intersections[0].mObject;
                        pointA = ray.Extrap(ray.intersections[0].mT);
                }
        }

  
        VssRay *vssRay  = NULL;

        bool validSample = (objectA != objectB);
        
        //if (validSample && thesame) Debug << "warning!!" << endl;

        if (0 && mDetectEmptyViewSpace) 
        {   
                // consider all samples valid
                // check if the viewpoint lies on the line segment AB
                if (Distance(pointA, pointB) <
                        Distance(viewPoint, pointA) + Distance(viewPoint, pointB) - Limits::Small) 
                {
                        validSample = false;
                }
        }
        
        if (validSample) 
        {       
                if (objectA) 
                {
                        vssRay = new VssRay(pointB,
                                                                pointA,
                                                                objectB,
                                                                objectA,
                                                                mPass);
                        vssRays.push_back(vssRay);
                        hits ++;
                }

                if (objectB) 
                {
                        vssRay = new VssRay(pointA,
                                                                pointB,
                                                                objectA,
                                                                objectB,
                                                                mPass);
                        vssRays.push_back(vssRay);
                        hits ++;
                }
                //Debug << "internal ray: " << *vssRay << endl << endl;
        }
//cout << "b";
        return hits;
}


int
VssPreprocessor::CastRay(
                                                 const Vector3 &viewPoint,
                                                 const Vector3 &direction,
                                                 VssRayContainer &vssRays,
                                                 const AxisAlignedBox3 &box
                                                 )
{
        switch (mRayCastMethod)
        {
        case INTEL_RAYCASTER:
                return CastIntelDoubleRay(viewPoint, direction, vssRays, box);
        case INTERNAL_RAYCASTER:
        default:
                return CastInternalRay(viewPoint, direction, vssRays, box);
        }
}


Vector3
VssPreprocessor::GetViewpoint(AxisAlignedBox3 *viewSpaceBox)
{
  AxisAlignedBox3 box;

  if (viewSpaceBox)
        box =*viewSpaceBox;
  else
        box = mKdTree->GetBox();

  // shrink the box in the y direction
  return box.GetRandomPoint();
}

Vector3
VssPreprocessor::GetDirection(const Vector3 &viewpoint,
                                                          AxisAlignedBox3 *viewSpaceBox
                                                          )
{
  Vector3 point;
  if (!use2dSampling) 
  {
          if (0) 
          {
                  Vector3 normal;
                  int i = Random((int)mObjects.size());
                  Intersectable *object = mObjects[i];
                  object->GetRandomSurfacePoint(point, normal);
          } 
          else
                  point = mKdTree->GetBox().GetRandomPoint();
        //        point = viewpoint + UniformRandomVector();
  } 
  else 
  {
          AxisAlignedBox3 box;

          if (viewSpaceBox)
                  box =*viewSpaceBox;
          else
                  box = mKdTree->GetBox();

          point = box.GetRandomPoint();
          point.y = viewpoint.y;
  }

  return point - viewpoint;
}

int
VssPreprocessor::GenerateImportanceRays(VssTree *vssTree,
                                                                                const int desiredSamples,
                                                                                SimpleRayContainer &rays
                                                                                )
{
  int num;
  if (0) {
        float minRayContribution;
        float maxRayContribution;
        float avgRayContribution;

        vssTree->GetRayContributionStatistics(minRayContribution,
                                                                                  maxRayContribution,
                                                                                  avgRayContribution);

        cout<<
          "#MIN_RAY_CONTRIB\n"<<minRayContribution<<endl<<
          "#MAX_RAY_CONTRIB\n"<<maxRayContribution<<endl<<
          "#AVG_RAY_CONTRIB\n"<<avgRayContribution<<endl;

        float p = desiredSamples/(float)(avgRayContribution*vssTree->stat.Leaves());
        num = vssTree->GenerateRays(p, rays);
  } else {
        int leaves = vssTree->stat.Leaves();
        num = vssTree->GenerateRays(desiredSamples, leaves, rays);
  }

  cout<<"Generated "<<num<<" rays."<<endl;

  return num;
}


bool
VssPreprocessor::ExportRays(const char *filename,
                                                        const VssRayContainer &vssRays,
                                                        const int number
                                                        )
{
  cout<<"Exporting vss rays..."<<endl<<flush;

  Exporter *exporter = NULL;
  exporter = Exporter::GetExporter(filename);
  exporter->SetWireframe();
  exporter->ExportKdTree(*mKdTree);
  exporter->SetFilled();
  exporter->ExportScene(mSceneGraph->mRoot);
  exporter->SetWireframe();

  if (mViewSpaceBox) {
        exporter->SetForcedMaterial(RgbColor(1,0,1));
        exporter->ExportBox(*mViewSpaceBox);
        exporter->ResetForcedMaterial();
  }

  VssRayContainer rays;
  vssRays.SelectRays(number, rays);
  
  //exporter->ExportRays(rays, RgbColor(1, 0, 0));

  delete exporter;

  cout<<"done."<<endl<<flush;

  return true;
}


bool
VssPreprocessor::ExportVssTree(char *filename,
                                                           VssTree *tree,
                                                           const Vector3 &dir
                                                           )
{
  Exporter *exporter = Exporter::GetExporter(filename);
  exporter->SetFilled();
  exporter->ExportScene(mSceneGraph->mRoot);
  //  exporter->SetWireframe();
  bool result = exporter->ExportVssTree2( *tree, dir );
  delete exporter;
  return result;
}

bool
VssPreprocessor::ExportVssTreeLeaf(char *filename,
                                                                   VssTree *tree,
                                                                   VssTreeLeaf *leaf)
{
  Exporter *exporter = NULL;
  exporter = Exporter::GetExporter(filename);
  exporter->SetWireframe();
  exporter->ExportKdTree(*mKdTree);

  if (mViewSpaceBox) {
        exporter->SetForcedMaterial(RgbColor(1,0,0));
        exporter->ExportBox(*mViewSpaceBox);
        exporter->ResetForcedMaterial();
  }

  exporter->SetForcedMaterial(RgbColor(0,0,1));
  exporter->ExportBox(tree->GetBBox(leaf));
  exporter->ResetForcedMaterial();

  VssRayContainer rays[4];
  for (int i=0; i < leaf->rays.size(); i++) {
        int k = leaf->rays[i].GetRayClass();
        rays[k].push_back(leaf->rays[i].mRay);
  }

  // SOURCE RAY
  exporter->ExportRays(rays[0], RgbColor(1, 0, 0));
  // TERMINATION RAY
  exporter->ExportRays(rays[1], RgbColor(1, 1, 1));
  // PASSING_RAY
  exporter->ExportRays(rays[2], RgbColor(1, 1, 0));
  // CONTAINED_RAY
  exporter->ExportRays(rays[3], RgbColor(0, 0, 1));

  delete exporter;
  return true;
}

void
VssPreprocessor::ExportVssTreeLeaves(VssTree *tree, const int number)
{
  vector<VssTreeLeaf *> leaves;
  tree->CollectLeaves(leaves);

  int num = 0;
  int i;
  float p = number / (float)leaves.size();
  for (i=0; i < leaves.size(); i++) {
        if (RandomValue(0,1) < p) {
          char filename[64];
          sprintf(filename, "vss-leaf-%04d.x3d", num);
          ExportVssTreeLeaf(filename, tree, leaves[i]);
          num++;
        }
        if (num >= number)
          break;
  }
}


void VssPreprocessor::TestBeamCasting(VssTree *tree, 
                                                                          ViewCellsManager *vm, 
                                                                          const ObjectContainer &objects)
{
        //debuggerWidget = new GlDebuggerWidget(renderer);
        //  renderer->resize(640, 480);
        //debuggerWidget->resize(640, 480);

        vector<VssTreeLeaf *> leaves;
        tree->CollectLeaves(leaves);

        Exporter *exporter = Exporter::GetExporter("shafts.x3d");

        exporter->SetWireframe();
        exporter->ExportGeometry(objects);
        exporter->SetFilled();
        //Randomize();
// §§matt
//      debuggerWidget = new GlDebuggerWidget(renderer);

        /*debuggerWidget->mBeam = beam;
        debuggerWidget->mSourceObject = sourceObj;
        debuggerWidget->mSamples = 10000;
        
        Debug << "showing window" << endl;
        debuggerWidget->show();
        
        renderer->makeCurrent();*/

        for (int i = 0; i < 10; ++ i)
        {
                Beam beam;
                Intersectable *sourceObj = mObjects[5];

                const int index = (int)RandomValue(0, (Real)((int)leaves.size() - 1));
                VssTreeLeaf *leaf = leaves[index];

                AxisAlignedBox3 dirBox = tree->GetDirBBox(leaf);
                AxisAlignedBox3 box = tree->GetBBox(leaf);
                
                beam.Construct(box, dirBox);

                // collect kd leaves and view cells
                mKdTree->CastBeam(beam);
                vm->CastBeam(beam);

                Debug << "found " << (int)beam.mViewCells.size() << " view cells and " 
                          << (int)beam.mKdNodes.size() << " kd nodes" << endl;

                BeamSampleStatistics stats;
// §§matt
/*              renderer->SampleBeamContributions(sourceObj,
                                                                                  beam,
                                                                                  200000,
                                                                                  stats);

                char s[64]; sprintf(s, "shaft%04d.png", i);

                QImage image = renderer->toImage();
                image.save(s, "PNG");
                Debug << "beam statistics: " << stats << endl << endl;
*/
                if (1)
                {
                        AxisAlignedBox3 sbox = mSceneGraph->GetBox();
                        Vector3 bmin = sbox.Min() - 150.0f;
                        Vector3 bmax = sbox.Max() + 150.0f;
                        AxisAlignedBox3 vbox(bmin, bmax);
                
                        exporter->ExportBeam(beam, vbox);
                }

                bool exportViewCells = false;
                
                if (exportViewCells)
                {
                        ViewCellContainer::const_iterator it, it_end = beam.mViewCells.end();
                        
                        for (it = beam.mViewCells.begin(); it != beam.mViewCells.end(); ++ it)
                        {
                                BspNodeGeometry geom;
                                AxisAlignedBox3 vbox;
                                vbox.Initialize();
                                vbox.Include((*it)->GetMesh());
                        
                                exporter->SetWireframe();
                                exporter->ExportBox(vbox);
                                exporter->SetFilled();
                                exporter->ExportViewCell(*it);
                        }

                        /*vector<KdNode *>::const_iterator it, it_end = beam.mKdNodes.end();
                        
                        for (it = beam.mKdNodes.begin(); it != beam.mKdNodes.end(); ++ it)
                        {
                                exporter->ExportBox(mKdTree->GetBox((*it)));
                        }*/
                }
        }
        /*while (1)
        { debuggerWidget->repaint();
        };*/
        delete exporter;
}


float
VssPreprocessor::GetAvgPvsSize(VssTree *tree,
                                                           const vector<AxisAlignedBox3> &viewcells
                                                           )
{
  vector<AxisAlignedBox3>::const_iterator it, it_end = viewcells.end();

  int sum = 0;
  for (it = viewcells.begin(); it != it_end; ++ it)
        sum += tree->GetPvsSize(*it);

  return sum/(float)viewcells.size();
}

bool
VssPreprocessor::ComputeVisibility()
{
        Debug << "type: vss" << endl;

  long startTime = GetTime();

  int totalSamples = 0;


  AxisAlignedBox3 box(mKdTree->GetBox());
  
  if (!useViewspacePlane) {
        float size = 0.05f;
        float s = 0.5f - size;
        float olds = Magnitude(box.Size());
        box.Enlarge(box.Size()*Vector3(-s));
        Vector3 translation = Vector3(-olds*0.1f, 0, 0);
        box.SetMin(box.Min() + translation);
        box.SetMax(box.Max() + translation);
  } else {

        // sample city like heights
        box.SetMin(1, box.Min(1) + box.Size(1)*0.2f);
        box.SetMax(1, box.Min(1) + box.Size(1)*0.3f);
  }

  if (use2dSampling)
        box.SetMax(1, box.Min(1));

  cout << "use view space box=" << mUseViewSpaceBox << endl;

  
  if (mUseViewSpaceBox)
  {
        //mViewSpaceBox = ConstructViewSpaceBox();

          if (!mEnlargeViewSpace)
          {
                  mViewSpaceBox = new AxisAlignedBox3(box);
          }
          else
          {
                  mViewSpaceBox = new AxisAlignedBox3(mKdTree->GetBox());

                  if (0) 
                  {
                          // HACK: enlarge in y directon
                          Vector3 size = mViewSpaceBox->Size();
                        
                          size[1] *= 1.25;
                          Vector3 enlarge(size[0] * 0.25f, size[1] * 0.0f, size[2] * 0.25f);
                          //Vector3 enlarge(size[0] * 4.0f, 0.0f, 0.0f);

                          mViewSpaceBox->Enlarge(enlarge);
                          mViewSpaceBox->SetMax(mViewSpaceBox->Max() + enlarge);
                  } 
                  else 
                  {
                          AxisAlignedBox3 tbox(*mViewSpaceBox);

                          Vector3 size = mViewSpaceBox->Size();
                          tbox.SetMax(0, mViewSpaceBox->Max(0) + size[0] * 0.5);
                          tbox.SetMin(0, mViewSpaceBox->Min(0) + size[0]);
                          *mViewSpaceBox = tbox;

                          // $$ JB temporary
                          /*AxisAlignedBox3 tempbox = *mViewSpaceBox;
                          
                          float s = tempbox.Size(0);
                          
                          tempbox.Scale(0.8f);
                          tempbox.SetMax(0, box.Max(0) + s*0.8f);
                          tempbox.SetMin(0, box.Min(0) + s*0.8f);
                          *mViewSpaceBox = tempbox;*/
                  }             
          }
          //Debug << "view space box: " << *mViewSpaceBox << endl;
  }
  else
  {
          mViewSpaceBox = NULL;
  }
  
  AxisAlignedBox3 vbox = mViewSpaceBox ? *mViewSpaceBox : mKdTree->GetBox();

  mSceneGraph->CollectObjects(&mObjects);

  //-- load view cells from file if requested
  if (!mLoadViewCells)
  {
          mViewCellsManager->SetViewSpaceBox(vbox);
          // construct view cells using it's own set of samples
          mViewCellsManager->Construct(this);

          //-- several visualizations and statistics
          Debug << "view cells construction finished: " << endl;
          mViewCellsManager->PrintStatistics(Debug);
  }
  else
  {      
          VssRayContainer dummies;
          mViewCellsManager->Visualize(mObjects, dummies);
          mViewCellsManager->ExportViewCells("test.xml", mViewCellsManager->GetExportPvs(), mObjects);
  }

  VssTree *vssTree = NULL;
  

  long initialTime = GetTime();

  if (mLoadInitialSamples)
  {
          cout << "Loading samples from file ... ";
          LoadSamples(mVssRays, mObjects);
          cout << "finished\n" << endl;
          totalSamples = (int)mVssRays.size();
  }
  else
  {
        
        while (totalSamples < mInitialSamples) {
                int passContributingSamples = 0;
                int passSampleContributions = 0;
                int passSamples = 0;

                int index = 0;

                int sampleContributions;

                int s = Min(mSamplesPerPass, mInitialSamples);
                for (int k=0; k < s; k++) {
                        // changed by matt
                        Vector3 viewpoint; 
                        //                      viewpoint = GetViewpoint(mViewSpaceBox);
                        mViewCellsManager->GetViewPoint(viewpoint);
                        Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);

                        sampleContributions = CastRay(viewpoint, direction, mVssRays, vbox);

                        if (sampleContributions) {
                                passContributingSamples ++;
                                passSampleContributions += sampleContributions;
                        }
                        passSamples++;
                        totalSamples++;
                }

                mPass++;
                int pvsSize = 0;
                float avgRayContrib = (passContributingSamples > 0) ?
                        passSampleContributions/(float)passContributingSamples : 0;

                cout << "#Pass " << mPass << " : t = " << TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
                cout << "#TotalSamples=" << totalSamples/1000
                        << "#SampleContributions=" << passSampleContributions << " ("
                        << 100*passContributingSamples/(float)passSamples<<"%)" << " avgPVS="
                        << pvsSize/(float)mObjects.size() << endl
                        << "avg ray contrib=" << avgRayContrib << endl;

                mStats <<
                        "#Pass\n" <<mPass<<endl<<
                        "#Time\n" << TimeDiff(startTime, GetTime())*1e-3 << endl<<
                        "#TotalSamples\n" << totalSamples<< endl<<
                        "#SampleContributions\n" << passSampleContributions << endl <<
                        "#PContributingSamples\n"<<100*passContributingSamples/(float)passSamples<<endl <<
                        "#AvgPVS\n"<< pvsSize/(float)mObjects.size() << endl <<
                        "#AvgRayContrib\n" << avgRayContrib << endl;
          }
 
          cout << "#totalPvsSize=" << mKdTree->CollectLeafPvs() << endl;


          
  }
 

  cout << "#totalRayStackSize=" << (int)mVssRays.size() << endl << flush;
  Debug << (int)mVssRays.size() << " rays generated in " 
            << TimeDiff(initialTime, GetTime()) * 1e-3 << " seconds" << endl;

  if (mStoreInitialSamples)
  {
          cout << "Writing " << (int)mVssRays.size() << " samples to file ... ";
          ExportSamples(mVssRays);
          cout << "finished\n" << endl;

          /*VssRayContainer dummyRays;
          LoadSamples(dummyRays, mObjects);
          Debug << "rays " << (int)mVssRays.size() << " " << dummyRays.size() << endl;

          for (int i = 0; i < (int)mVssRays.size(); ++ i)
          {
                  Debug << mVssRays[i]->GetOrigin() << " " << mVssRays[i]->GetTermination() << " " << mVssRays[i]->mOriginObject << " " << mVssRays[i]->mTerminationObject << endl;
                  Debug << dummyRays[i]->GetOrigin() << " " << dummyRays[i]->GetTermination() << " " << dummyRays[i]->mOriginObject << " " << dummyRays[i]->mTerminationObject << endl << endl;
          }*/
  }

  
  //int numExportRays = 2000;
  int numExportRays = 0;

  if (numExportRays) {
        char filename[64];
        sprintf(filename, "vss-rays-initial.x3d");
        ExportRays(filename, mVssRays, numExportRays);
  }

  vssTree = new VssTree;
  // viewcells = Construct(mVssRays);

  vssTree->Construct(mVssRays, mViewSpaceBox);
  cout<<"VssTree root PVS size = "<<vssTree->GetRootPvsSize()<<endl;

  ExportRays("kdtree.x3d", mVssRays, 10);

  if (0)
  {
          ExportVssTree("vss-tree-100.x3d", vssTree, Vector3(1,0,0));
          ExportVssTree("vss-tree-001.x3d", vssTree, Vector3(0,0,1));
          ExportVssTree("vss-tree-101.x3d", vssTree, Vector3(1,0,1));
          ExportVssTree("vss-tree-101m.x3d", vssTree, Vector3(-1,0,-1));
          ExportVssTreeLeaves(vssTree, 10);
  }

  // viewcells->UpdatePVS(newVssRays);
  // get viewcells as kd tree boxes
  vector<AxisAlignedBox3> kdViewcells;
  if (0) {
        vector<KdLeaf *> leaves;
        mKdTree->CollectLeaves(leaves);
        vector<KdLeaf *>::const_iterator it;
        int targetLeaves = 50;
        float prob = targetLeaves/(float)leaves.size();
        for (it = leaves.begin(); it != leaves.end(); ++it)
          if (RandomValue(0.0f,1.0f) < prob)
                kdViewcells.push_back(mKdTree->GetBox(*it));

        float avgPvs = GetAvgPvsSize(vssTree, kdViewcells);
        cout<<"Initial average PVS size = "<<avgPvs<<endl;
  }


  int samples = 0;
  int pass = 0;

  
  // cast view cell samples
  while (samples < mVssSamples) 
  {
        
        int num = mVssSamplesPerPass;
        SimpleRayContainer rays;
        VssRayContainer vssRays;

        if (!mUseImportanceSampling) {
          for (int j=0; j < num; j++) {
            // changed by matt
                //Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
                Vector3 viewpoint; 
                mViewCellsManager->GetViewPoint(viewpoint);
                Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
                rays.push_back(SimpleRay(viewpoint, direction));
          }
        } else {
          num = GenerateImportanceRays(vssTree, num, rays);
        }

        for (int i=0; i < rays.size(); i++)
          CastRay(rays[i].mOrigin, rays[i].mDirection, vssRays, vbox);

        vssTree->AddRays(vssRays);

        if (0) {
          int subdivided = vssTree->UpdateSubdivision();
          cout<<"subdivided leafs = "<<subdivided<<endl;
        }

        float avgPvs = GetAvgPvsSize(vssTree, kdViewcells);
        cout<<"Average PVS size = "<<avgPvs<<endl;

        /// compute view cell contribution of rays
        mViewCellsManager->ComputeSampleContributions(vssRays, true, false);
        
        if (numExportRays) {
          char filename[64];
          if (mUseImportanceSampling)
                sprintf(filename, "vss-rays-i%04d.x3d", pass);
          else
                sprintf(filename, "vss-rays-%04d.x3d", pass);

          ExportRays(filename, vssRays, numExportRays);
        }

        samples+=num;
        float pvs = vssTree->GetAvgPvsSize();
        cout<<"*****************************\n";
        cout<<samples<<" avgPVS ="<<pvs<<endl;
        cout<<"VssTree root PVS size = "<<vssTree->GetRootPvsSize()<<endl;
        cout<<"*****************************\n";
//      if (samples >= mVssSamples) break;
        pass ++;
  }

  if (mTestBeamSampling && mUseGlRenderer)
  {     
          TestBeamCasting(vssTree, mViewCellsManager, mObjects);
  }

  if (0)  Debug << vssTree->stat << endl;

  if (0)
  {
        VssRayContainer viewCellRays;
        // compute rays used for view cells construction
        const int numRays = mViewCellsManager->GetVisualizationSamples();
        vssTree->CollectRays(viewCellRays, numRays);
  }

  //-- render simulation after merge
  cout << "\nevaluating bsp view cells render time after sampling ... ";
  Debug << "\nStatistics after sampling: " << endl;

  mRenderSimulator->RenderScene();
  SimulationStatistics ss;
  mRenderSimulator->GetStatistics(ss);
  
  cout << " finished" << endl;
  cout << ss << endl;
  Debug << ss << endl;

  delete vssTree;
  
  return true;
}

}