source: trunk/VUT/GtpVisibilityPreprocessor/src/VssPreprocessor.cpp @ 436

#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"


bool useViewSpaceBox = true;//true;
bool use2dSampling = false;
bool useViewspacePlane = true;//true;

VssPreprocessor::VssPreprocessor():
  mPass(0),
  mVssRays()
{
  // this should increase coherence of the samples
  environment->GetIntValue("VssPreprocessor.samplesPerPass", mSamplesPerPass);
  environment->GetIntValue("VssPreprocessor.initialSamples", mInitialSamples);
  environment->GetIntValue("VssPreprocessor.vssSamples", mVssSamples);
  environment->GetIntValue("VssPreprocessor.vssSamplesPerPass", mVssSamplesPerPass);
  environment->GetBoolValue("VssPreprocessor.useImportanceSampling", mUseImportanceSampling);
  environment->GetIntValue("BspTree.Construction.samples", mBspConstructionSamples);
        
  mStats.open("stats.log");
}

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

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

int
VssPreprocessor::CastRay(
                                                 Vector3 &viewPoint,
                                                 Vector3 &direction,
                                                 VssRayContainer &vssRays
                                                 )
{
  int hits = 0;
  static Ray ray;
  AxisAlignedBox3 box = mKdTree->GetBox();

  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, pointB;
  float bsize = Magnitude(box.Size());
  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;
        box.ComputeMinMaxT(ray, &tmin, &tmax);
        if (tmax > bsize) {
          //                    cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
          //                    cerr<<"ray"<<ray<<endl;
        }
        pointA = ray.Extrap(tmax);
                
  }

  bool detectEmptyViewSpace = true;
        
  if (detectEmptyViewSpace) {
        SetupRay(ray, pointA, -direction);
  } else
        SetupRay(ray, viewPoint, -direction);
        
        
  if (mKdTree->CastRay(ray)) {
                
        objectB = ray.intersections[0].mObject;
        pointB = ray.Extrap(ray.intersections[0].mT);

  } else {
        objectB = NULL;
        float tmin, tmax;
        box.ComputeMinMaxT(ray, &tmin, &tmax);
        if (tmax > bsize) {
          //                    cerr<<"Warning: tmax > box size tmax="<<tmax<<" tmin="<<tmin<<" size="<<bsize<<endl;
          //                    cerr<<"ray"<<ray<<endl;
        }
                
        pointB = ray.Extrap(tmax);
  }

  VssRay *vssRay  = NULL;

  bool validSample = true;
  if (detectEmptyViewSpace) {
        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);
          vssRays.push_back(vssRay);
          hits ++;
        }
                
        if (objectB) {
          vssRay = new VssRay(pointA,
                                                  pointB,
                                                  objectA,
                                                  objectB);
          vssRays.push_back(vssRay);
          hits ++;
        }
  }
        
  return hits;
}


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) {
        Vector3 normal;
        int i = RandomValue(0, mObjects.size()-1);
        Intersectable *object = mObjects[i];
        object->GetRandomSurfacePoint(point, normal);
  } 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()/2;
        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;
        
  float prob = number/(float)vssRays.size();


  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,0));
        exporter->ExportBox(*mViewSpaceBox);
        exporter->ResetForcedMaterial();
  }
        
  VssRayContainer rays; for (int i=0; i < vssRays.size(); i++)
        if (RandomValue(0,1) < prob)
          rays.push_back(vssRays[i]);

  exporter->ExportRays(rays, RgbColor(1, 0, 0));
        
  delete exporter;

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

  return true;
}


bool
VssPreprocessor::ExportVssTree(char *filename,
                                                           VssTree *tree)
{
  Exporter *exporter = Exporter::GetExporter(filename);
  exporter->SetFilled();
  exporter->ExportScene(mSceneGraph->mRoot);
  exporter->SetWireframe();
  bool result = exporter->ExportVssTree(*tree);
  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;
  }
}


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()
{
  
  mSceneGraph->CollectObjects(&mObjects);
        
  long startTime = GetTime();
  
  int totalSamples = 0;


  AxisAlignedBox3 *box = new AxisAlignedBox3(mKdTree->GetBox());

  if (!useViewspacePlane) {
        float size = 0.01f;
        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.1);
        box->SetMax(1, box->Min(1) + box->Size(1)*0.2);
  }

  if (use2dSampling)
        box->SetMax(1, box->Min(1));
        
  if (useViewSpaceBox)
        mViewSpaceBox = box;
  else
        mViewSpaceBox = NULL;
                

  VssTree *vssTree = NULL;

  RayContainer bspRays;

  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++) {
                        
          Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
          Vector3 direction = GetDirection(viewpoint, mViewSpaceBox);
                        
          sampleContributions = CastRay(viewpoint, direction, mVssRays);
                        
                        
          //-- CORR matt: put block inside loop
          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 
                 << "k   #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=" << mVssRays.size() << endl <<flush;
        
  //int numExportRays = 10000;
  int numExportRays = 0;

  if (numExportRays) {
        char filename[64];
        sprintf(filename, "vss-rays-initial.x3d");
        ExportRays(filename, mVssRays, numExportRays);
  }
        
  //-- construct BSP view cells
  if (ViewCell::sHierarchy == ViewCell::BSP) 
  {
          const int bspSamples = min((int)mVssRays.size(), mBspConstructionSamples);
                
          Debug << "bpssamples: " << bspSamples << endl;
          for (int i = 0; i < bspSamples; ++ i)
                  bspRays.push_back(new Ray(*mVssRays[i]));
         
          //-- construct BSP tree using the samples
          mBspTree = new BspTree(&mUnbounded);  

          ObjectContainer objects;
          mSceneGraph->CollectObjects(&objects);
          
          mBspTree->SetGenerateViewCells(true);
          mBspTree->Construct(objects, bspRays);
        }

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

  ExportVssTree("vss-tree.x3d", vssTree);

  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;
  while (1) {
        int num = mVssSamplesPerPass;
        SimpleRayContainer rays;
        VssRayContainer vssRays;
                
        if (!mUseImportanceSampling) {
          for (int j=0; j < num; j++) {
                Vector3 viewpoint = GetViewpoint(mViewSpaceBox);
                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);
                
        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;

        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);
        }

        // cast rays into BSP tree
        if (ViewCell::sHierarchy == ViewCell::BSP) 
          {
                for (int i = 0; i < (int)vssRays.size(); ++ i)
                  {
                        CastRay(*mBspTree, *vssRays[i]);
                  }
          }

        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++;
  }
        
  delete vssTree;

  if (ViewCell::sHierarchy == ViewCell::BSP) 
        {
          Debug << mBspTree->GetStatistics();

          ObjectContainer objects;
          ExportSplits(objects, bspRays, 10000);
          ExportBspPvs(objects, bspRays, 10000);

          BspViewCellsStatistics stat;
          mBspTree->EvaluateViewCellsStats(stat);
          Debug << "original view cell partition:\n" << stat << endl;

          // clear BSP samples
          CLEAR_CONTAINER(bspRays);
        }

  return true;
}

void VssPreprocessor::CastRay(const BspTree &tree, const VssRay & vssRay)
{
  //-- cast ray to BSP tree to get intersection with view cells
  Ray ray(vssRay);
  mBspTree->CastRay(ray);
                                
  //if (0 && ray.sourceObject.mObject)
  //sampleContributions += 
  //AddObjectSamples(ray.sourceObject.mObject, ray);
                        
  if (!ray.intersections.empty()) // second intersection found
  {
          //sampleContributions += 
          AddObjectSamples(ray.intersections[0].mObject, ray);
  }
}

int VssPreprocessor::AddObjectSamples(Intersectable *obj, const Ray &ray)
{
  int contributingSamples = 0;
  int j;
 
  // object can be seen from the view cell => add to view cell pvs
  for (j=0; j < ray.bspIntersections.size(); ++ j)
        {       
          BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
          // if ray not in unbounded space
          if (leaf->GetViewCell() != &mUnbounded)
                contributingSamples += 
                  leaf->GetViewCell()->GetPvs().AddSample(obj);
        }
 
  // rays passing through this viewcell
  if (mPass > 1)
        for (j=1; j < ((int)ray.bspIntersections.size() - 1); ++ j) 
          {
                BspLeaf *leaf = ray.bspIntersections[j].mLeaf;
                
                if (leaf->GetViewCell() != &mUnbounded)
                  leaf->GetViewCell()->
                        AddPassingRay(ray, contributingSamples ? 1 : 0);
          }
  
  return contributingSamples;
}