source: GTP/trunk/Lib/Vis/Preprocessing/src/havran/testrt.cpp @ 2602

#define USE_THREADS 0

#ifdef UNICODE
#undef UNICODE
#endif

#define NOMINMAX
#ifdef __WINDOWS__
#include <windows.h>
#ifdef _CRT_SET
#include <crtdbg.h>
#endif // _CRT_SET
#endif
#include <cstdio>

#include "Camera.h"
#include "PreprocessorFactory.h"
#include "Parser.h"
#include "Environment.h"
#include "MeshKdTree.h"
#include "Preprocessor.h"
#include "common.h"
#include "PreprocessorThread.h"
#include "ObjExporter.h"
#include "SceneGraph.h"
#include "GlobalLinesRenderer.h"
#include "RayCaster.h"
#include "timer.h"
#include "raypack.h"

#include "ViewCellsManager.h"

#ifdef USE_QT  
        #include "QtPreprocessorThread.h"
        #include "QtGlViewer.h"
        #include "QtGlRenderer.h"
#else
        #if USE_THREADS
#include "BoostPreprocessorThread.h"
        #endif
#endif

#include "ResourceManager.h"
#include "GlRenderer.h"


#define USE_EXE_PATH false


using namespace GtpVisibilityPreprocessor;

//Preprocessor *preprocessor = NULL;
extern GlRendererWidget *rendererWidget;
//GlobalLinesRenderer *globalLinesRenderer = NULL;

// DLL function signature
typedef GlRendererWidget *(*importFunction)(Preprocessor *);


// Flag if the ray should be single sided (one direction)
// or double sided (shooting ray in both directions)
static bool castDoubleRays;

extern void Cleanup();

static string ReplaceSuffix(const string &filename, const string &a, const string &b)
{
        string result = filename;

        int pos = (int)filename.rfind(a, (int)filename.size() - 1);
        if (pos == filename.size() - a.size()) 
        {
                result.replace(pos, a.size(), b);
        }

        return result;
}


static int SplitFilenames(const string &str, vector<string> &filenames)
{
        int pos = 0;

        while(1) {
                int npos = (int)str.find(';', pos);
                
                if (npos < 0 || npos - pos < 1)
                        break;
                filenames.push_back(string(str, pos, npos - pos));
                pos = npos + 1;
        }
        
        filenames.push_back(string(str, pos, str.size() - pos));
        return (int)filenames.size();
}


static string GetInternFilename(const string &filename, const string newSuffix)
{
        vector<string> filenames;
        const int files = SplitFilenames(filename, filenames);

        vector<string>::const_iterator sit, sit_end = filenames.end();
        string kdFilename;

        int i = 0;
        for (sit = filenames.begin(); sit != sit_end; ++ sit, ++ i)
        {
                string currentFile = *sit;
                string strippedFilename;

                if (i == 0)
                {       
                        strippedFilename = currentFile;
                }
                else
                {
                        char *str = StripPath(currentFile.c_str());
                        strippedFilename = string(str);

                        delete [] str;
                }
                
                string suffix("_");

                if (i == (int)filenames.size() - 1)
                {
                        suffix = newSuffix;
                }

                if (strstr(strippedFilename.c_str(), ".x3d")) 
                {
                        kdFilename += ReplaceSuffix(strippedFilename, ".x3d", suffix);
                }
        else if (strstr(strippedFilename.c_str(), ".dat"))
                {
                        kdFilename += ReplaceSuffix(strippedFilename, ".dat", suffix);
                } 
                else if (strstr(strippedFilename.c_str(), ".obj"))
                {
                        kdFilename += ReplaceSuffix(strippedFilename, ".obj", suffix);
                }
                else
                {
                        cerr << "Error: Currently unsupported format for kd, filename " << currentFile << endl;
                }
        }

        //cout << "kdfilename: " << kdFilename << endl;
        return kdFilename;
}


// -------------------------------------------------------
// Written by Vlastimil Havran

// This is for testing RT implementation
void
TestRTcamera(int argc, char **argv)
{
  int returnCode = 0;

  InitTiming();
  Debug.open("debug.log");

  Environment::GetSingleton()->Parse(argc, argv, USE_EXE_PATH);
  MeshKdTree::ParseEnvironment();

  char buff[128];
  Environment::GetSingleton()->GetStringValue("Preprocessor.type", buff);
  string preprocessorType(buff);

  if (!(preprocessor = PreprocessorFactory::CreatePreprocessor(preprocessorType)))
  {
    Environment::DelSingleton();
    cerr << "Unknown preprocessor type" << endl;
    exit(1);
  }


  Environment::GetSingleton()->GetStringValue("Scene.filename", buff);
  string filename(buff);

  const string externKdTree = ReplaceSuffix(filename, ".obj", ".kdf");
  const string internKdTree = GetInternFilename(filename, preprocessor->mLoadMeshes ?
                                                ".kdm" : ".kdt");

  if (preprocessor->InitRayCast(externKdTree, internKdTree))
  {
    cout << "ray casting initialized!" << endl;
  }
  else
  {
    cout << "ray casting initialization failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  //Debug << "using pvs type " << PVS_TYPE << endl;

  /////////////
  //-- load scene 
  
  if (!preprocessor->LoadScene(filename))
  {
    cout << "loading file " << filename << " failed" << endl;
    Cleanup();
    exit(1);
  }
        
  ////////////
  //-- initialize external ray caster
  
  if (preprocessor->LoadInternKdTree(internKdTree))
  {
    cout << "intern kd tree loaded!" << endl;
  }
  else
  {
    cout << "loading intern kd tree failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  // export objects as obj
  if (preprocessor->mExportObj)
  {
    if (strstr(filename.c_str(), ".obj"))
      {
        cerr << "already in obj format" << endl;
        if (0)  preprocessor->ExportObj("test.obj", preprocessor->mObjects);
      }
    else
      {
        
        const string objname = GetInternFilename(filename, ".obj");
        
        cout << "exporting scene to " << objname << endl;
        bool success = preprocessor->ExportObj(objname, preprocessor->mObjects);
        
        
        if (success)
          cout << "finished exporting obj" << endl;
        else
          cerr << "error exporting " << objname << endl;
      }
  }

  int width = 1000;
  int height = 500;
  float fieldOfView = 115.f;
  Camera cam(width, height, fieldOfView);
  
  AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();
  AxisAlignedBox3 bbox = bboxOrig;
  int sizeDiag = Magnitude(bbox.Diagonal());
  bbox.Enlarge(sizeDiag * 1.5f);

  Vector3 origin, dir;

//#define WIEN1
#define WIEN2
//#define ARENA1
    
#ifdef WIEN1
  // 1099.9 183.0 -387
  origin = Vector3(1099.9f, 183.0f, -387.0f);
  dir = Vector3(-0.6f, 0.0001f, -0.8f);
#define DIREXISTS
#endif  
#ifdef WIEN2
  // 935.6 215.8 -1012.3
  origin = Vector3(935.6, 215.8, -1012.4);
  dir = Vector3(0.01f, 0.01f, 1.0f);
#define DIREXISTS
#endif
#ifdef ARENA1
  origin = Vector3(22.16, 19.63, -950.44);
  dir = Vector3(0.9, 0.001f, -0.5f);
#define DIREXISTS
#endif  

  cam.SetPosition(origin);

#ifndef DIREXISTS
  Vector3 center = bbox.Center();
  dir = center - origin;
#endif
  
  dir.Normalize();
  cam.SetDirection(dir);
  cout << "Computing image\n" << endl;

#if 1
  // ray by ray
  cam.SnapImage("test-rays.tga",
                preprocessor->mRayCaster,
                bboxOrig,
                preprocessor->mSceneGraph);
#endif
#if 1
  // using ray packets
  cam.SnapImagePacket("test-packet.tga",
                      preprocessor->mRayCaster,
                      bboxOrig,
                      preprocessor->mSceneGraph);
#endif
#if 1
  // using ray packets
  cam.SnapImage2("test-oneDir.tga",
                 preprocessor->mRayCaster,
                 bboxOrig,
                 preprocessor->mSceneGraph);
#endif
  
  cout << "Done\n" << endl;
  return;
}

struct RESult {
  int hitA;
  float hitAT;
  int hitB;
  float hitBT;
  RESult(int hA, float tA, int hB, float tB):
    hitA(hA), hitAT(tA), hitB(hB), hitBT(tB) { }
};

// This is for testing RT implementation
void
TestRTfromFile(int argc, char **argv)
{
  int returnCode = 0;

  InitTiming();
  Debug.open("debug.log");

  Environment::GetSingleton()->Parse(argc, argv, USE_EXE_PATH);
  MeshKdTree::ParseEnvironment();

  char buff[128];
  Environment::GetSingleton()->GetStringValue("Preprocessor.type", buff);
  string preprocessorType(buff);

  if (!(preprocessor = PreprocessorFactory::CreatePreprocessor(preprocessorType)))
  {
    Environment::DelSingleton();
    cerr << "Unknown preprocessor type" << endl;
    exit(1);
  }


  Environment::GetSingleton()->GetStringValue("Scene.filename", buff);
  string filename(buff);

  const string externKdTree = ReplaceSuffix(filename, ".obj", ".kdf");
  const string internKdTree = GetInternFilename(filename, preprocessor->mLoadMeshes ?
                                                ".kdm" : ".kdt");

  if (preprocessor->InitRayCast(externKdTree, internKdTree))
  {
    cout << "ray casting initialized!" << endl;
  }
  else
  {
    cout << "ray casting initialization failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  //Debug << "using pvs type " << PVS_TYPE << endl;

  /////////////
  //-- load scene 
  
  if (!preprocessor->LoadScene(filename))
  {
    cout << "loading file " << filename << " failed" << endl;
    Cleanup();
    exit(1);
  }
        
  ////////////
  //-- initialize external ray caster
  
  if (preprocessor->LoadInternKdTree(internKdTree))
  {
    cout << "intern kd tree loaded!" << endl;
  }
  else
  {
    cout << "loading intern kd tree failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  // export objects as obj
  if (preprocessor->mExportObj)
  {
    if (strstr(filename.c_str(), ".obj"))
    {
      cerr << "already in obj format" << endl;
      if (0)    preprocessor->ExportObj("test.obj", preprocessor->mObjects);
    }
    else
    {
      const string objname = GetInternFilename(filename, ".obj");
        
      cout << "exporting scene to " << objname << endl;
      bool success = preprocessor->ExportObj(objname, preprocessor->mObjects);  
        
      if (success)
        cout << "finished exporting obj" << endl;
      else
        cerr << "error exporting " << objname << endl;
    }
  }

  Environment::GetSingleton()->GetStringValue("Rays.file", buff);
  FILE *fp = fopen(buff, "rb");
  if (!fp) {
    cerr << "ERROR: file " << buff << " cannot be opened for reading" << endl;
    cerr << "EXITING" << endl;
    exit(3);
  }
  cout << "File " << buff << " was opened successfully" << endl;
  int cntMaxRays = 100000;
  Environment::GetSingleton()->GetIntValue("Rays.cnt", cntMaxRays);
  vector<SimpleRay> rays;
  SimpleRay rayTest;
  vector<RESult> results;
  
  int cntRays = 0;
  for (int i = 0; cntRays < cntMaxRays; i++) {
    int ch = fgetc(fp);
    switch (ch) {
    case 'G': { // two-sided rays of 16
      for (int j = 0; j < 16; j++) {
        int order;
        float ox, oy, oz, dx, dy, dz;
        int hitA; float hitAT;
        int hitB; float hitBT;  
        if (fscanf(fp, "%d %f %f %f %f %f %f %d %f %d %f\n",
                   &order, &ox, &oy, &oz, &dx, &dy, &dz, &hitA, &hitAT, &hitB, &hitBT) != 11) {
          cerr << "Problem parsing the ray file" << endl;
          cerr << "EXITING for ray order" << order << endl;
          goto FINISH;
        }
        Vector3 mPosition(ox, oy, oz);
        Vector3 mDirection(dx, dy, dz);
        rayTest.Set(mPosition, mDirection, 0, 1.0, true);
        rays.push_back(rayTest);
        results.push_back(RESult(hitA, hitAT, hitB, hitBT));
        cntRays++;
      }
      break;
    }
    case 'H': // one-sided rays of 16
      cerr << "Not yet implemented " << endl; abort();
      
    case 'D': // two-sided ray
      cerr << "Not yet implemented " << endl; abort();
      
    case 'S': // one-sided ray
      cerr << "Not yet implemented " << endl; abort();
    } // switch
  } // for i

FINISH:
  fclose(fp);

  Environment::GetSingleton()->GetBoolValue("TestDoubleRays", castDoubleRays);

  double mult = 1.0;
  if (castDoubleRays)
    mult = 2.0;
  
  cout << "Starting to shoot " << cntRays * mult << " rays" << endl;
  
  RayCaster *raycaster = preprocessor->mRayCaster;
  VssRayContainer vssRays;
  AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();

//#define DEBUGRESULTS

#ifdef DEBUGRESULTS
  int cntIntelYesWeNo = 0;
  int cntIntelNoWeYes = 0;
  int cntIntelOurDistMismatch = 0;
  bool printOut = false;
  double sumDistIntel = 0.f;
  double sumDistOurAlg = 0.f;
  double sumDistAbsDiff = 0.f;
#endif

  cout << "Press a key to start ray shooting" << endl;
  getchar();
  cout << "Ray shooting " << cntRays << " rays started - "
       << (castDoubleRays ? " double " : " single ")
       << " dir " << endl;

  long t1 = GetTime();
  CTimer timer;
  timer.Reset();
  timer.Start();
  
  SimpleRayContainer raysToTest;
  for (int i = 0; i < cntRays - 16; i++) {
#if 0
    int res = raycaster->CastRay(rays[i],
                                 vssRays,
                                 bboxOrig,
                                 castDoubleRays, // castDoubleRay,
                                 false); // pruneInvalidRays
#else
    raysToTest.erase(raysToTest.begin(), raysToTest.end());
    for (int j = 0; j < 16; j++, i++) {
      raysToTest.push_back(rays[i]);
    }
#if 0
    for (int j = 0; j < 16; j++) {
      cout << "orig = " << raysToTest[j].mOrigin << " dir = "
           << raysToTest[j].mDirection << endl;
    }
#endif    
        
    raycaster->CastRays16(raysToTest,
                          vssRays,
                          bboxOrig,
                          castDoubleRays,
                          false);
#endif
    
#ifdef DEBUGRESULTS
    if (!castDoubleRays) {
      float T = rays[i].IntersectionRes[0].tdist;
      if ( (res != results[i].hitA) ||
           ((res)&&(results[i].hitA)&& (fabs(T - results[i].hitAT)) > 0.1f)
           )
        {
          if (printOut)
            cout << " i = " << i << " ";
          if ((res != 0) && (results[i].hitA == 0)) {
            cntIntelNoWeYes++;
            if (printOut)
              cout << "Intel no intersection, our alg intersection at t = " << T << endl;
          }
          else {
            if ((res == 0) && (results[i].hitA)) {
              cntIntelYesWeNo++;
              if (printOut)
                cout << "Intel intersection at = " << results[i].hitAT
                     << " , our alg no intersection = " << endl;
              
            }
            else {
              cntIntelOurDistMismatch++;
              sumDistOurAlg += T;
              sumDistIntel += results[i].hitAT;
              sumDistAbsDiff += fabs(T - results[i].hitAT);
              if (printOut)
                cout << "Intel intersection at = " << results[i].hitAT
                     << " , our alg intersection at = " << T << endl;       
            }
          }
        }
    } // double rays
    else {
      // checking for results of double rays - not yet implemented
    }
#endif    
  }

  timer.Stop();
  long t2 = GetTime();
  cout<<"\n#RAY_CAST_TIME = ";
  cout << TimeDiff(t1, t2)<<" [mikrosec]"
       << " userTime = " << timer.UserTime() << " realTime = "
       << timer.RealTime() << endl;
  cout << "Rays shot per milisecond [userTimer] = "
       << ((double)cntRays * mult/(double)timer.UserTime()) / 1000.f << endl;

#ifdef DEBUGRESULTS  
  cout << "cntIntelYesWeNo = " << cntIntelYesWeNo << endl;
  cout << "cntIntelNoWeYes = " << cntIntelNoWeYes << endl;
  cout << "cntIntelOurDistMismatch = " << cntIntelOurDistMismatch++ << endl;
  cout << " sumDistIntel = " << sumDistIntel/(double)cntIntelOurDistMismatch << endl;
  cout << " sumDistOur = " << sumDistOurAlg/(double)cntIntelOurDistMismatch << endl;
  cout << " sumAbsDiffDist = " << sumDistAbsDiff/(double)cntIntelOurDistMismatch
       << endl;
#endif  
  cout << "Done\n" << endl;
  return;
}

//-------------------------------------------------------------------------
// This is for testing RT implementation using ray packets
void
TestRTfromFilePackets(int argc, char **argv)
{
  int returnCode = 0;

  InitTiming();
  Debug.open("debug.log");

  Environment::GetSingleton()->Parse(argc, argv, USE_EXE_PATH);
  MeshKdTree::ParseEnvironment();

  char buff[128];
  Environment::GetSingleton()->GetStringValue("Preprocessor.type", buff);
  string preprocessorType(buff);

  if (!(preprocessor = PreprocessorFactory::CreatePreprocessor(preprocessorType)))
  {
    Environment::DelSingleton();
    cerr << "Unknown preprocessor type" << endl;
    exit(1);
  }


  Environment::GetSingleton()->GetStringValue("Scene.filename", buff);
  string filename(buff);

  const string externKdTree = ReplaceSuffix(filename, ".obj", ".kdf");
  const string internKdTree = GetInternFilename(filename, preprocessor->mLoadMeshes ?
                                                ".kdm" : ".kdt");

  if (preprocessor->InitRayCast(externKdTree, internKdTree))
  {
    cout << "ray casting initialized!" << endl;
  }
  else
  {
    cout << "ray casting initialization failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  //Debug << "using pvs type " << PVS_TYPE << endl;

  /////////////
  //-- load scene 
  
  if (!preprocessor->LoadScene(filename))
  {
    cout << "loading file " << filename << " failed" << endl;
    Cleanup();
    exit(1);
  }
        
  ////////////
  //-- initialize external ray caster
  
  if (preprocessor->LoadInternKdTree(internKdTree))
  {
    cout << "intern kd tree loaded!" << endl;
  }
  else
  {
    cout << "loading intern kd tree failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  // export objects as obj
  if (preprocessor->mExportObj)
  {
    if (strstr(filename.c_str(), ".obj"))
    {
      cerr << "already in obj format" << endl;
      if (0)    preprocessor->ExportObj("test.obj", preprocessor->mObjects);
    }
    else
    {
      const string objname = GetInternFilename(filename, ".obj");
        
      cout << "exporting scene to " << objname << endl;
      bool success = preprocessor->ExportObj(objname, preprocessor->mObjects);  
        
      if (success)
        cout << "finished exporting obj" << endl;
      else
        cerr << "error exporting " << objname << endl;
    }
  }

  Environment::GetSingleton()->GetStringValue("Rays.file", buff);
  FILE *fp = fopen(buff, "rb");
  if (!fp) {
    cerr << "ERROR: file " << buff << " cannot be opened for reading" << endl;
    cerr << "EXITING" << endl;
    exit(3);
  }
  cout << "File " << buff << " was opened successfully" << endl;
  int cntMaxRays = 100000;
  Environment::GetSingleton()->GetIntValue("Rays.cnt", cntMaxRays);
  vector<SimpleRay> rays;
  SimpleRay rayTest;
  vector<RESult> results;
  
  int cntRays = 0;
  for (int i = 0; cntRays < cntMaxRays; i++) {
    int ch = fgetc(fp);
    switch (ch) {
    case 'G': { // two-sided rays of 16
      for (int j = 0; j < 16; j++) {
        int order;
        float ox, oy, oz, dx, dy, dz;
        int hitA; float hitAT;
        int hitB; float hitBT;  
        if (fscanf(fp, "%d %f %f %f %f %f %f %d %f %d %f\n",
                   &order, &ox, &oy, &oz, &dx, &dy, &dz, &hitA, &hitAT, &hitB, &hitBT) != 11) {
          cerr << "Problem parsing the ray file" << endl;
          cerr << "EXITING for ray order" << order << endl;
          goto FINISH;
        }
        Vector3 mPosition(ox, oy, oz);
        Vector3 mDirection(dx, dy, dz);
        rayTest.Set(mPosition, mDirection, 0, 1.0, true);
        rays.push_back(rayTest);
        results.push_back(RESult(hitA, hitAT, hitB, hitBT));
        cntRays++;
      }
      break;
    }
    case 'H': // one-sided rays of 16
      cerr << "Not yet implemented " << endl; abort();
      
    case 'D': // two-sided ray
      cerr << "Not yet implemented " << endl; abort();
      
    case 'S': // one-sided ray
      cerr << "Not yet implemented " << endl; abort();
    } // switch
  } // for i
FINISH:
  fclose(fp);
  
  Environment::GetSingleton()->GetBoolValue("TestDoubleRays", castDoubleRays);

  double mult = 1.0;
  if (castDoubleRays)
    mult = 2.0;

  cout << "Starting to shoot " << cntRays * mult << " rays" << endl;
  
  RayCaster *raycaster = preprocessor->mRayCaster;
  VssRayContainer vssRays;
  AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();

#undef DEBUGRESULTS
//#define DEBUGRESULTS

#ifdef DEBUGRESULTS
  int cntIntelYesWeNo = 0;
  int cntIntelNoWeYes = 0;
  int cntIntelOurDistMismatch = 0;
  bool printOut = false;
  double sumDistIntel = 0.;
  double sumDistOurAlg = 0.;
  double sumDistAbsDiff = 0.f;
#endif

  cout << "Press a key to start ray shooting" << endl;
  getchar();
  cout << "Ray packs shooting " << cntRays << " rays started - "
       << (castDoubleRays ? " double " : " single ")
       << " dir " << endl;

  long t1 = GetTime();
  CTimer timer;
  timer.Reset();
  timer.Start();

#ifdef _USE_HAVRAN_SSE
#ifdef __SSE__  
  RayPacket2x2 raysPack;
  
  for (int i = 0; i < cntRays - 16; i++) {
    for (int j = 0; j < 4; j++, i++) {
      raysPack.SetLoc(j, rays[i].mOrigin);
      raysPack.SetDir(j, rays[i].mDirection);
    } // for
#if 0
    for (int j = 0; j < 16; j++) {
      cout << "orig = " << raysToTest[j].mOrigin << " dir = "
           << raysToTest[j].mDirection << endl;
    }
#endif
    raysPack.ComputeDirSign();
    raycaster->CastRaysPacket2x2(raysPack, castDoubleRays);
    
#ifdef DEBUGRESULTS
    if (!castDoubleRays) {
      float T = rays[i].IntersectionRes[0].tdist;
      if ( (res != results[i].hitA) ||
           ((res)&&(results[i].hitA)&& (fabs(T - results[i].hitAT)) > 2.0f)
           )
        {
          if (printOut)
            cout << " i = " << i << " ";
          if ((res != 0) && (results[i].hitA == 0)) {
            cntIntelNoWeYes++;
            if (printOut)
              cout << "Intel no intersection, our alg intersection at t = " << T << endl;
          }
          else {
            if ((res = 0) && (results[i].hitA)) {
              cntIntelYesWeNo++;
              if (printOut)
                cout << "Intel intersection at = " << results[i].hitAT
                     << " , our alg no intersection = " << endl;
              
            }
            else {
              cntIntelOurDistMismatch++;
              sumDistOurAlg += T;
              sumDistIntel += results[i].hitAT;
              sumDistAbsDiff += fabs(T - results[i].hitAT);
              if (printOut)
                cout << "Intel intersection at = " << results[i].hitAT
                     << " , our alg intersection at = " << T << endl;       
            }
          }
        }
    } // double rays
    else {
      // checking for results of double rays - not yet implemented
    }
#endif    
  }

#endif // __SSE__ 
#endif // _USE_HAVRAN_SSE

  timer.Stop();
  long t2 = GetTime();
  cout<<"\n#RAY_CAST_TIME = ";
  cout << TimeDiff(t1, t2)<<" [mikrosec]"
       << " userTime = " << timer.UserTime() << " realTime = "
       << timer.RealTime() << endl;
    ;
  cout << "Rays shot per milisecond [userTimer] = "
       << ((double)cntRays * mult/(double)timer.UserTime()) / 1000.f << endl;

#ifdef DEBUGRESULTS  
  cout << "cntIntelYesWeNo = " << cntIntelYesWeNo << endl;
  cout << "cntIntelNoWeYes = " << cntIntelNoWeYes << endl;
  cout << "cntIntelOurDistMismatch = " << cntIntelOurDistMismatch++ << endl;
  cout << " sumDistIntel = " << sumDistIntel/(double)cntIntelOurDistMismatch << endl;
  cout << " sumDistOur = " << sumDistOurAlg/(double)cntIntelOurDistMismatch << endl;
  cout << " sumAbsDiffDist = " << sumDistAbsDiff/(double)cntIntelOurDistMismatch
       << endl;
#endif  
  cout << "Done\n" << endl;
  return;
}

// This is for testing RT implementation
void
TestRT_4_fromFile(int argc, char **argv)
{
  int returnCode = 0;

  InitTiming();
  Debug.open("debug.log");

  Environment::GetSingleton()->Parse(argc, argv, USE_EXE_PATH);
  MeshKdTree::ParseEnvironment();

  char buff[128];
  Environment::GetSingleton()->GetStringValue("Preprocessor.type", buff);
  string preprocessorType(buff);

  if (!(preprocessor = PreprocessorFactory::CreatePreprocessor(preprocessorType)))
  {
    Environment::DelSingleton();
    cerr << "Unknown preprocessor type" << endl;
    exit(1);
  }


  Environment::GetSingleton()->GetStringValue("Scene.filename", buff);
  string filename(buff);

  const string externKdTree = ReplaceSuffix(filename, ".obj", ".kdf");
  const string internKdTree = GetInternFilename(filename, preprocessor->mLoadMeshes ?
                                                ".kdm" : ".kdt");

  if (preprocessor->InitRayCast(externKdTree, internKdTree))
  {
    cout << "ray casting initialized!" << endl;
  }
  else
  {
    cout << "ray casting initialization failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  //Debug << "using pvs type " << PVS_TYPE << endl;

  /////////////
  //-- load scene 
  
  if (!preprocessor->LoadScene(filename))
  {
    cout << "loading file " << filename << " failed" << endl;
    Cleanup();
    exit(1);
  }
        
  ////////////
  //-- initialize external ray caster
  
  if (preprocessor->LoadInternKdTree(internKdTree))
  {
    cout << "intern kd tree loaded!" << endl;
  }
  else
  {
    cout << "loading intern kd tree failed!" << endl;
    Cleanup();
    exit(1);
  }
  
  // export objects as obj
  if (preprocessor->mExportObj)
  {
    if (strstr(filename.c_str(), ".obj"))
    {
      cerr << "already in obj format" << endl;
      if (0)    preprocessor->ExportObj("test.obj", preprocessor->mObjects);
    }
    else
    {
      const string objname = GetInternFilename(filename, ".obj");
        
      cout << "exporting scene to " << objname << endl;
      bool success = preprocessor->ExportObj(objname, preprocessor->mObjects);  
        
      if (success)
        cout << "finished exporting obj" << endl;
      else
        cerr << "error exporting " << objname << endl;
    }
  }

  Environment::GetSingleton()->GetStringValue("Rays.file", buff);
  FILE *fp = fopen(buff, "rb");
  if (!fp) {
    cerr << "ERROR: file " << buff << " cannot be opened for reading" << endl;
    cerr << "EXITING" << endl;
    exit(3);
  }
  cout << "File " << buff << " was opened successfully" << endl;
  int cntMaxRays = 100000;
  Environment::GetSingleton()->GetIntValue("Rays.cnt", cntMaxRays);
  vector<SimpleRay> rays;
  SimpleRay rayTest;
  vector<RESult> results;
  vector<AxisAlignedBox3> boxes;

  int cntRays = 0;
  for (int i = 0; cntRays < cntMaxRays;) {
    int ch = fgetc(fp);
    switch (ch) {
    case 'I': { // two-sided rays of 16
      Vector3 minv, maxv;
      if (fscanf(fp, "%f %f %f %f %f %f\n",
                 &minv.x, &minv.y, &minv.z,
                 &maxv.x, &maxv.y, &maxv.z) != 6) {
        cerr << "Problem parsing the ray file" << endl;
        cerr << "EXITING for ray order" << i << endl;
        goto FINISH;
      }
      boxes.push_back(AxisAlignedBox3(minv, maxv));      

      for (int j = 0; j < 4; j++) {
        int order;
        float ox, oy, oz, dx, dy, dz;
        int hitA; float hitAT;
        int hitB; float hitBT;  
        if (fscanf(fp, "%d %f %f %f %f %f %f %d %f\n",
                   &order, &ox, &oy, &oz, &dx, &dy, &dz,
                   &hitA, &hitAT) != 9) {
          cerr << "Problem parsing the ray file" << endl;
          cerr << "EXITING for ray order" << order << endl;
          goto FINISH;
        }
        Vector3 mPosition(ox, oy, oz);
        Vector3 mDirection(dx, dy, dz);
        rayTest.Set(mPosition, mDirection, 0, 1.0, true);
        rays.push_back(rayTest);
        results.push_back(RESult(hitA, hitAT, 0, -1.f));
        cntRays++;
        i++;
      }
      break;
    }
    default: {
      cerr << "Not yet implemented or end of file" << endl;
      goto FINISH;
    }
    } // switch
  } // for i

FINISH:
  
  fclose(fp);
  cout << "Starting to shoot " << cntRays << " rays" << endl;
  
  RayCaster *raycaster = preprocessor->mRayCaster;
  VssRayContainer vssRays;
  AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();

//#define DEBUGRESULTS

#ifdef DEBUGRESULTS
  int cntIntelYesWeNo = 0;
  int cntIntelNoWeYes = 0;
  int cntIntelOurDistMismatch = 0;
  bool printOut = false;
  double sumDistIntel = 0.f;
  double sumDistOurAlg = 0.f;
  double sumDistAbsDiff = 0.f;
#endif

  //cout << "Press a key to start ray shooting" << endl;
  //getchar();
  cout << "Ray shooting " << cntRays << " rays started - "
       << (castDoubleRays ? " double " : " single ")
       << " dir " << endl;

  long t1 = GetTime();
  CTimer timer;
  timer.Reset();
  timer.Start();
  Vector3 boxMin, boxMax;
  
  SimpleRayContainer raysToTest;
  int boxI = 0;
  for (int i = 0; i < cntRays - 4; i+= 4, boxI++) {
    Vector3 origin4[4];
    Vector3 direction4[4];
    int     result4[4];
    float   dist4[4];
    boxMin = boxes[boxI].Min();
    boxMax = boxes[boxI].Max();
    for (int j = 0; j < 4; j++) {
      int o = i+j;
      origin4[j] = rays[o].mOrigin;
      direction4[j] = rays[o].mDirection;
    }
    raycaster->CastRaysPacket4(boxMin, boxMax,
                               origin4, direction4,
                               result4, dist4);
    printf("I %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f\n",
           boxMin.x, boxMin.y, boxMin.z, boxMax.x, boxMax.y, boxMax.z);

    for (int j = 0; j < 4; j++) {
      printf("%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
             i+j,
             origin4[j].x,
             origin4[j].y,
             origin4[j].z,
             direction4[j].x,
             direction4[j].y,
             direction4[j].z,
             (result4[j] != -1) ? 1 : 0,
             (result4[j] != -1) ? dist4[j] : 0);      
    } // for j    
  } // for i

  timer.Stop();
  long t2 = GetTime();
  cout<<"\n#RAY_CAST_TIME = ";
  cout << TimeDiff(t1, t2)<<" [mikrosec]"
       << " userTime = " << timer.UserTime() << " realTime = "
       << timer.RealTime() << endl;
  cout << "Rays shot per milisecond [userTimer] = "
       << ((double)cntRays/(double)timer.UserTime()) / 1000.f << endl;

#ifdef DEBUGRESULTS  
  cout << "cntIntelYesWeNo = " << cntIntelYesWeNo << endl;
  cout << "cntIntelNoWeYes = " << cntIntelNoWeYes << endl;
  cout << "cntIntelOurDistMismatch = " << cntIntelOurDistMismatch++ << endl;
  cout << " sumDistIntel = " << sumDistIntel/(double)cntIntelOurDistMismatch << endl;
  cout << " sumDistOur = " << sumDistOurAlg/(double)cntIntelOurDistMismatch << endl;
  cout << " sumAbsDiffDist = " << sumDistAbsDiff/(double)cntIntelOurDistMismatch
       << endl;
#endif  
  cout << "Done\n" << endl;
  return;
}