BVH_SAH

#define DEFINE_VARS
#include "OnlineHelpKeys.h"
#include <list>
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <stdio.h>

#include <iostream>
#include <algorithm>
#include <vector>
#include <cfloat>
#include <string>

#include <cstdio>
#include <cstdlib>
#include <sstream>
#include <fstream>
#include <cfloat>
//using SDL
#include <SDL.h>

//using LIB3DS
#include <types.h>
#include <material.h>
#include <mesh.h>
#include <file.h>
#include <map>
using namespace std;

#define AMBIENT     96.f
#define BVH_STACK_SIZE 32
#define MAX_RAY_DEPTH       3

#define DIFFUSE     128.f
#define SPECULAR    192.f

// Ray intersections of a distance <=NUDGE_FACTOR (from the origin) don't count
#define NUDGE_FACTOR     1e-5f
// Maximum allowed depth of BVH
// Checked at BVH build time, no runtime crash possible, see below
#define BVH_STACK_SIZE 32
#define TRI_MAGIC   0xDEADBEEF
#define TRI_MAGICNORMAL 0xDEADC0DE
#define DEBUG_LOG_BVH

using namespace std;

#define WIDTH       800
#define HEIGHT      600
#define SCREEN_DIST (HEIGHT*2)

#define PROGRESS_REPORT
#ifdef  PROGRESS_REPORT
#define REPORT(x) x
#define REPORTPRM(x) x,
#else
#define REPORT(x)
#define REPORTPRM(x)
#endif
using std::string;
unsigned g_reportCounter = 0;
typedef float coord;

struct Float3
{
   float X, Y Z;
   Float3(float x=0, float y=0, float z=0)
        :
   X(x), Y(y), Z(z){}

   Float3(const Float3& rhs)
        :
   X(rhs.X), Y(rhs.Y), Z(rhs.Z){}

   inline Float3& operator+=(const Float3& rhs)
   {    X += rhs.X; Y += rhs.Y; Z += rhs.Z; return *this; }

   inline Float3& operator-=(const Float3& rhs)
    {X -= rhs.X; Y -= rhs.Y; Z -= rhs.Z; return *this;}

   inline Float3& operator*=(const float& rhs)
    {X *= rhs; Y *= rhs; Z *= rhs; return *this;}

   inline Float3& operator/=(const float& rhs)
    {X /= rhs; Y /= rhs; Z /= rhs; return *this;}

   inline bool operator!=(const Float3& rhs)
    {return X!=rhs.X || Y!=rhs.Y || Z!=rhs.Z; }

   inline Float3 operator*(float rhs) const
    {return Float3(X*rhs, Y*rhs, Z*rhs);   }

   inline Float3 operator+(Float3& rhs) const
    {return Float3(X+rhs.X, Y+rhs.Y, Z+rhs.Z);}

   inline Float3 operator-(Float3& rhs) const
    {return Float3(X-rhs.X, Y-rhs.Y, Z-rhs.Z);   }

   inline void assignSmaller(const Float3& rhs)
   {  X = std::min(X, rhs.X); Y = std::min(Y, rhs.Y); Z = std::min(Z, rhs.Z);   }

   inline void assignBigger(const Float3& rhs)
    { X = std::max(X, rhs.X); Y = std::max(Y, rhs.Y); Z = std::max(Z, rhs.Z);}

   inline float length()
    {return sqrt(X*X +Y*Y + Z*Z); }

   inline float lengthsq()
    {return X*X +Y*Y + Z*Z; }

   inline void Normalize()
    {   float norm = length();
    X/= norm; Y/= norm; Z/= norm;   }
};

struct Vertex : public Float3
{
    Float3 n;
    unsigned _ambientOcclusionCoeff;

    Vertex(float x, float y, float z, float nx, float ny, float nz, unsigned char amb=60)
    :
    Float3(x,y,z), n(nx,ny,nz), _ambientOcclusionCoeff(amb){}

};

struct Pixel {
    float _b, _g, _r;
    Pixel(float r=0.f, float g=0.f, float b=0.f)
    :
    _b(b), _g(g), _r(r) {}

    Pixel& operator+=(const Pixel& rhs) { _b += rhs._b; _g += rhs._g; _r += rhs._r; return *this; }
    Pixel& operator-=(const Pixel& rhs) { _b -= rhs._b; _g -= rhs._g; _r -= rhs._r; return *this; }
    Pixel& operator*=(const float& rhs) { _b =  rhs*_b; _g  = rhs*_g; _r  = rhs*_r; return *this; }
    Pixel& operator/=(const float& rhs) { _b =  _b/rhs; _g  = _g/rhs; _r  = _r/rhs; return *this; }

    Pixel operator+(const Pixel& rhs) {
    float r = _r+rhs._r; if (r<0.f) r=0.f; if (r>255.f) r=255.f;
    float g = _g+rhs._g; if (g<0.f) g=0.f; if (g>255.f) g=255.f;
    float b = _b+rhs._b; if (b<0.f) b=0.f; if (b>255.f) b=255.f;
    return Pixel(r,g,b);
    }
    Pixel operator*(const float& rhs) { return Pixel(rhs*_r, rhs*_g, rhs*_b); }
};

struct Screen
{
    static SDL_Surface *_surface;
    float _Zbuffer[HEIGHT][WIDTH];
    const struct Scene& _scene;

    public:

    Screen(const struct Scene& scene)
    :
    _scene(scene)
    {
    if ( SDL_Init(SDL_INIT_VIDEO) < 0 )
    {
        std::cerr << "Couldn't initialize SDL: " <<  SDL_GetError() << std::endl;
        exit(0);
    }

    // Clean up on exit
    atexit(SDL_Quit);

    // We ask for 32bit surface...
    _surface = SDL_SetVideoMode( WIDTH, HEIGHT, 32, SDL_HWSURFACE | SDL_DOUBLEBUF | SDL_HWACCEL | SDL_ASYNCBLIT);
    if (!_surface)
    // ...and if that fails, we settle for a software-emulation of 16bit
    _surface = SDL_SetVideoMode( WIDTH, HEIGHT, 16, SDL_SWSURFACE | SDL_DOUBLEBUF);

    if (!_surface)
    {
       std::cerr << "Couldn't set video mode: " << SDL_GetError() << std::endl;
       exit(0);
    }

    if ( SDL_MUSTLOCK(_surface) )
    {
       if ( SDL_LockSurface(_surface) < 0 )
       {
        std::cerr << "Couldn't lock _surface: " << SDL_GetError() << std::endl;
        exit(0);
       }
    }
ClearScreen();
ClearZbuffer();
}

~Screen()
{
    if ( SDL_MUSTLOCK(_surface) )
    SDL_UnlockSurface(_surface);
}

void ClearScreen()
{
    Uint32 color = SDL_MapRGB(_surface->format, 0, 0, 0);
    SDL_FillRect(_surface, NULL, color);
}

void ClearZbuffer()
{
    memset(reinterpret_cast<void*>(&_Zbuffer[0][0]), 0x0, sizeof(_Zbuffer));
}

    // Almost verbatim copy from www.libsdl.org "Introduction" section
void DrawPixel(int y, int x, Uint32 color)
    {
    switch (_surface->format->BytesPerPixel)
    {
        case 1:
        { /* Assuming 8-bpp */
          Uint8 *bufp;

          bufp = (Uint8 *)_surface->pixels + y*_surface->pitch + x;
          *bufp = color;
            }
        break;

        case 2:
             { /* Probably 15-bpp or 16-bpp */
        Uint16 *bufp;

        bufp = (Uint16 *)_surface->pixels + y*_surface->pitch/2 + x;
        *bufp = color;
         }
        break;

        case 3:
              { /* Slow 24-bpp mode, usually not used */
        Uint8 *bufp;

        bufp = (Uint8 *)_surface->pixels + y*_surface->pitch + x;
        *(bufp+_surface->format->Rshift/8) =
            ((color & _surface->format->Rmask) >> _surface->format->Rshift) << _surface->format->Rloss;
        *(bufp+_surface->format->Gshift/8) =
            ((color & _surface->format->Gmask) >> _surface->format->Gshift) << _surface->format->Gloss;
        *(bufp+_surface->format->Bshift/8) =
            ((color & _surface->format->Bmask) >> _surface->format->Bshift) << _surface->format->Bloss;
         }
        break;

        case 4:
         { /* Probably 32-bpp */
        Uint32 *bufp;

        bufp = (Uint32 *)_surface->pixels + y*_surface->pitch/4 + x;
        *bufp = color;
        }
        break;
    }
    }

void ShowScreen(bool raytracerOutput=false)
{
    if (!raytracerOutput)
      { //  to add the "Press H for help" when not benchmarking
    extern bool g_benchmark;
    if (!g_benchmark)
    {
      unsigned char *pData = onlineHelpKeysImage;
       for(int h=0; h<OHELPH; h++)
          for(int w=0; w<OHELPW; w++)
        {
               pData++;
           unsigned char c = *pData++;
           pData++;
           if (c<200)
           DrawPixel(h + 20,WIDTH-20-OHELPW + w,
                SDL_MapRGB(_surface->format, 255-c, 255-c, 255-c));
        }

         }
     }
     if ( SDL_MUSTLOCK(_surface) ) {SDL_UnlockSurface(_surface);}

     if (raytracerOutput)
    SDL_UpdateRect(_surface, 0, 0, 0, 0);
     else
        SDL_Flip(_surface);

     if ( SDL_MUSTLOCK(_surface) )
    {
       if ( SDL_LockSurface(_surface) < 0 )
        {
           std::cerr << "Couldn't lock _surface: " << SDL_GetError() << std::endl;
           exit(0);
        }
    }
}

};
SDL_Surface *Screen::_surface = NULL;

struct Triangle
{
    const Vertex *_vertexA, *_vertexB, *_vertexC;
    Float3 centroid, n;

    // Color:
    Pixel _colorf;
    // precomputed for SDL surface
    Uint32 _color;

    // Should we backface cull this triangle?
    bool _twoSided;

    // Raytracing intersection pre-computed cache:
    float _d, _d1, _d2, _d3;
    Float3 _e1, _e2, _e3, _bottom, _top;

    Triangle(
    const Vertex *vertexA,
        const Vertex *vertexB,
    const Vertex *vertexC,
    unsigned r, unsigned g, unsigned b,
    bool twosided = false, bool triNormalProvided = false,
    Float3 triNormal =Float3(0.,0.,0.));
};

Triangle::Triangle(
    const Vertex *vertexA,
        const Vertex *vertexB,
    const Vertex *vertexC,
    unsigned r, unsigned g, unsigned b,
    bool twosided , bool triNormalProvided,
    Float3 triNormal)
        :
        _vertexA(vertexA), _vertexB(vertexB), _vertexC(vertexC),

        centroid((vertexA->X + vertexB->X + vertexC->X)/3.0f,
         (vertexA->Y + vertexB->Y + vertexC->Y)/3.0f,
         (vertexA->Z + vertexB->Z + vertexC->Z)/3.0f),

    _colorf((float)r,(float)g,(float)b), // For use in all other cases
        _color(SDL_MapRGB(Screen::_surface->format, r,g,b)), // For use with DrawPixel
        _twoSided(twosided),
    _bottom(FLT_MAX,FLT_MAX, FLT_MAX), // Will be updated after centering in Loader
        _top(-FLT_MAX, -FLT_MAX,-FLT_MAX) // Will be updated after centering in Loader
    {// this is where the debugger takes me to//
        if (!triNormalProvided)
        {
           n = Float3((vertexA->n.X + vertexB->n.X + vertexC->n.X)/3.0f,
                  (vertexA->n.Y + vertexB->n.Y + vertexC->n.Y)/3.0f,
                  (vertexA->n.Z + vertexB->n.Z + vertexC->n.Z)/3.0f);
           n.Normalize();
        }
        else {n = triNormal;}
}

struct Bound{
    Float3 minBound;
    Float3 maxBound;
    Float3 center;
    float dia;
}b;

struct Camera {
    Float3 From;
    Float3 At;
    Float3 Up;
};

Camera Cam =  {Float3(0.0f,1.0f/3.0f,2.0f),
           Float3(0.0f,1.0f/3.0f,0.0f),
           Float3(0.0f,1.0f,0.0f)};

struct  Ray{
    Float3 ray_d;
    Float3 ray_o;
    float tMin, tMax;
};

struct Co_ordinateFrame
{
    Float3 x_axis;
    Float3 y_axis;
    Float3 z_axis;
    Float3 origin;
};
Co_ordinateFrame CF;

Float3 Scale(Float3 vec1,float s)
{
    Float3 res;
    res.X = vec1.X*s;
    res.Y = vec1.Y*s;
    res.Z = vec1.Z*s;

    return res;
}
float Dot(Float3 vec1,Float3 vec2)
{
    float dot = vec1.X*vec2.X+vec1.Y*vec2.Y +vec1.Z*vec2.Z;
    return dot;
}

Float3 Cross(Float3 vec1,Float3 vec2)
{
     Float3 cross;

     cross.X = (vec1.Y*vec2.Z) - (vec1.Z*vec2.Y);
     cross.Y = (vec1.Z*vec2.X) - (vec1.X*vec2.Z);
     cross.Z = (vec1.X*vec2.Y) - (vec1.Y*vec2.X);

     return cross;
}

 Float3 normalize(Float3 vec)
 {

     float mag = sqrt(Dot(vec,vec));
     Float3 result;
     result.X = vec.X /mag;
     result.Y = vec.Y/mag;
     result.Z = vec.Z/mag;

     return result;
 }
 void  ComputeCameraframe()
 {
    CF.origin = Cam.From;
    CF.z_axis = normalize(Cam.From - Cam.At);
    CF.x_axis = normalize(Cross(Cam.Up,CF.z_axis));
        CF.y_axis = normalize(Cross(CF.z_axis,CF.x_axis));
 }
float distancesq(const Float3& a, const Float3& b)
{

    float dx=a.X - b.X;
    float dy=a.Y - b.Y;
    float dz=a.Z - b.Z;
    return dx*dx + dy*dy + dz*dz;
}

float Distance(const Float3& a, const Float3& b)
{

    float dx=a.X - b.X;
    float dy=a.Y - b.Y;
    float dz=a.Z - b.Z;
    return sqrt(dx*dx + dy*dy + dz*dz);
}

struct BVHNode {
    Float3 bottom;
    Float3 top;
    virtual bool IsLeaf() = 0;
};


struct BVHInner : BVHNode {
    BVHNode *left;
    BVHNode *right;
    virtual bool IsLeaf() { return false; }
};

struct BVHLeaf : BVHNode {
    Float3 leafbottom;
    Float3 leaftop;
        std::list<const Triangle*> _triangles;
        virtual bool IsLeaf() { return true; }
};

struct CacheFriendlyBVHNode
{
    Float3 _bottom;
    Float3 _top;
    union
    {
         struct {
          unsigned _idxLeft;
          unsigned _idxRight;
        } inner;
    struct {
         unsigned _count; // Top-most bit set
         unsigned _startIndexInTriIndexList;
           } leaf;
        } u;
};
// Bounding Volume Hierarchy
// Cache-friendly version of the Bounding Volume Hierarchy data
// (32 bytes per CacheFriendlyBVHNode, i.e. one CPU cache line)

struct Scene {
    static const float MaxCoordAfterRescale;

    std::vector<Vertex>    _vertices;
    std::vector<Triangle>  _triangles;

    // Bounding Volume Hierarchy
    BVHNode *_pSceneBVH;

    // Cache-friendly version of the Bounding Volume Hierarchy data
    // (32 bytes per CacheFriendlyBVHNode, i.e. one CPU cache line)
    unsigned _triIndexListNo;
    int *_triIndexList;
    unsigned _pCFBVH_No;
    CacheFriendlyBVHNode *_pCFBVH;

    Scene()
    :
    _pSceneBVH(NULL),
    _triIndexListNo(-1),
    _triIndexList(NULL),
    _pCFBVH_No(-1),
    _pCFBVH(NULL)
    {}

    // Load object
    void load(const char *filename,Bound *b, Camera *cam);

    // Update triangle normals
    void fixns(void);

    // Cache-friendly version of the Bounding Volume Hierarchy data
    // (creation functions)
    void PopulateCacheFriendlyBVH(
    Triangle *pFirstTriangle,
    BVHNode *root,
    unsigned& idxBoxes,
    unsigned& idxTriList);
    void CreateCFBVH();

    // Creates BVH and Cache-friendly version of BVH
    void UpdateBoundingVolumeHierarchy(const char *filename, bool forceRecalc=false);

    // Since it may be aborted for taking too long, this returns "boolCompletedOK"
    bool renderRaytracer(Camera& eye, Screen& canvas, bool antiAlias = false);
};

//BVHNode *CreateBVH(const Scene *pScene);
//void CreateCFBVH(Scene *);
// Work item for creation of BVH:

struct BBoxTmp {
    // Bottom point (ie minx,miny,minz)
    Float3 boxbottom;
    // Top point (ie maxx,maxy,maxz)
    Float3 boxtop;
    // Center point, ie 0.5*(top-bottom)
    Float3 boxcenter;
    // Triangle
    const Triangle *pTri;
    BBoxTmp()
    :
    boxbottom(FLT_MAX,FLT_MAX,FLT_MAX),
    boxtop(-FLT_MAX,-FLT_MAX,-FLT_MAX),
    pTri(NULL)
    {}
    };

typedef vector<BBoxTmp> BBoxEntries;

// This builds the BVH, finding optimal split planes for each depth
BVHNode *Recurse(BBoxEntries& work, REPORTPRM(float pct=0.0f) int depth=0)
{
    REPORT( float pctSpan  = 11.0f/pow(3.0f,depth); )
    if (work.size() < 4)
    {
        BVHLeaf *leaf = new BVHLeaf;
        for(BBoxEntries::iterator it=work.begin(); it!=work.end(); it++)
            leaf->_triangles.push_back(it->pTri);
        return leaf;
        }
    // Start by finding the working list's bounding box
    Float3 bottom(FLT_MAX,FLT_MAX,FLT_MAX), top(-FLT_MAX,-FLT_MAX,-FLT_MAX);

    for(unsigned i=0; i<work.size(); i++)
    {
      BBoxTmp& v = work[i];
      bottom.assignSmaller(v.boxbottom);
      top.assignBigger(v.boxtop);
        }

    float side1 = top.X-bottom.X;
    float side2 = top.Y-bottom.Y;
    float side3 = top.Z-bottom.Z;
    // The current box has a cost of (No of triangles)*surfaceArea
    float minCost = work.size() * (side1*side2 + side2*side3 + side3*side1);
    float bestSplit = FLT_MAX; // will indicate no split with better cost found (below)
    int bestAxis = -1;

    // Try all different axis
    for (int j=0; j<3; j++)
    {
      int axis = j;

    // try dividing the triangles based on the current axis,
    // and  split values from "start" to "stop", one "step" at a time.
      float start, stop, step;

      if (axis==0)
       {
        start = bottom.X;
        stop  = top.X;
       }
      else if (axis==1)
       {
        start = bottom.Y;
        stop  = top.Y;
       }
    else
      {
        start = bottom.Z;
        stop  = top.Z;
      }
    // In that axis, do the bounding boxes in the work queue "span" across?
    // Or are they all already "packed" on the axis's plane?

        if (fabsf(stop-start)<1e-4)  continue;// if so move to a different axis!

    // Try splitting at a uniform sampling that gets smaller the deeper u go:
    // size of "sampling grid": 1024 (depth 0), 512 (depth 1), etc

        step = (stop-start)/(1024.0f/(depth+1.0f));
    // track the Progress: the Progress report variables...

        #ifdef PROGRESS_REPORT
        float pctStart = pct + j*pctSpan;
        float pctStep  = pctSpan/((stop - start - 2*step)/step);
        #endif

    float testSplit;
    for(testSplit = start+step; testSplit<stop-step; testSplit+=step)
      {
        pctStart += pctStep;
        #ifdef PROGRESS_REPORT
        if ((1023&g_reportCounter++) == 0)
          {
                    printf( "\b\b\b%02d%%", int(pctStart));
            fflush(stdout);
          }
        pctStart += pctStep;
        #endif
        // The left and right bounding box
       Float3 lbottom(FLT_MAX,FLT_MAX,FLT_MAX),
              ltop(-FLT_MAX,-FLT_MAX,-FLT_MAX),
              rbottom(FLT_MAX,FLT_MAX,FLT_MAX),
              rtop(-FLT_MAX,-FLT_MAX,-FLT_MAX);

// The number of triangles in the left and right bboxes

   int countLeft=0, countRight=0;
// For each test split, allocate triangles based on their bounding boxes centers

     for(unsigned i=0; i<work.size(); i++)
    {
      BBoxTmp& v = work[i];
          float value;
      if (axis==0) value=v.boxcenter.X;
      else if (axis==1) value=v.boxcenter.Y;
      else value=v.boxcenter.Z;

        if (value < testSplit)
        {
            lbottom.assignSmaller(v.boxbottom);
            ltop.assignBigger(v.boxtop);
                    countLeft++;
        }
        else
        {
            rbottom.assignSmaller(v.boxbottom);
            rtop.assignBigger(v.boxtop);
            countRight++;
        }
        }
        // Now use the Surface Area Heuristic to see if this split has a better "cost"
        // First, check for flase partitionings
        if (countLeft<=1 || countRight<=1) continue;
        // It's a real partitioning, calculate the surface areas
        float lside1 = ltop.X-lbottom.X;
        float lside2 = ltop.Y-lbottom.Y;
        float lside3 = ltop.Z-lbottom.Z;
        float rside1 = rtop.X-rbottom.X;
        float rside2 = rtop.Y-rbottom.Y;
        float rside3 = rtop.Z-rbottom.Z;
        float surfaceLeft  = lside1*lside2 + lside2*lside3 + lside3*lside1;
        float surfaceRight = rside1*rside2 + rside2*rside3 + rside3*rside1;
        float totalCost = surfaceLeft*countLeft + surfaceRight*countRight;

        if (totalCost < minCost)
        {
            minCost = totalCost;
            bestSplit = testSplit;
            bestAxis = axis;
             }
    }
    }
    // We found no split to improve the cost :( create a BVH leaf
    if (bestAxis == -1)
    {
        BVHLeaf *leaf = new BVHLeaf;
        for(BBoxEntries::iterator it=work.begin(); it!=work.end(); it++)
        {
          leaf->_triangles.push_back(it->pTri);
          return leaf;
        }
    }
// Create a BVH inner node, split with the optimal value we found above
BBoxEntries left, right;
Float3 lbottom(FLT_MAX,FLT_MAX,FLT_MAX),
       ltop(-FLT_MAX,-FLT_MAX,-FLT_MAX),
       rbottom(FLT_MAX,FLT_MAX,FLT_MAX),
       rtop(-FLT_MAX,-FLT_MAX,-FLT_MAX);

    for(int i=0; i<(int)work.size(); i++)
    {
       BBoxTmp& v = work[i];

       float value;
       if (bestAxis==0) value=v.boxcenter.X;
       else if (bestAxis==1) value=v.boxcenter.Y;
       else value=v.boxcenter.Z;

      if (value < bestSplit)
       {
        #ifdef DEBUG_LOG_BVH
        printf("LADD: B(%f %f %f) T(%f %f %f) C(%f %f %f)\n",
            v.boxbottom.X,
            v.boxbottom.Y,
            v.boxbottom.Z,

            v.boxtop.X,
            v.boxtop.Y,
            v.boxtop.Z,

            v.boxcenter.X,
            v.boxcenter.Y,
            v.boxcenter.Z);
        #endif
            left.push_back(v);
            lbottom.assignSmaller(v.boxbottom);
            ltop.assignBigger(v.boxtop);
      }
       else
      {
        #ifdef DEBUG_LOG_BVH
        printf("RADD: B(%f %f %f) T(%f %f %f) C(%f %f %f)\n",
            v.boxbottom.X,
            v.boxbottom.Y,
            v.boxbottom.Z,

            v.boxtop.X,
            v.boxtop.Y,
            v.boxtop.Z,

            v.boxcenter.X,
            v.boxcenter.Y,
            v.boxcenter.Z);
        #endif
            right.push_back(v);
            rbottom.assignSmaller(v.boxbottom);
            rtop.assignBigger(v.boxtop);
        }
    }

    BVHInner *inner = new BVHInner;
    #ifdef PROGRESS_REPORT
     if ((1023&g_reportCounter++) == 0)
    {
       printf("\b\b\b%2d%%", int(pct+3.f*pctSpan)); // Update progress indicator
       fflush(stdout);
        }
    #endif
    inner->left = Recurse(left, REPORTPRM(pct+3.f*pctSpan) depth+1);
    inner->left->bottom = lbottom;
    inner->left->top = ltop;
    #ifdef PROGRESS_REPORT
    if ((1023&g_reportCounter++) == 0)
    {
        printf("\b\b\b%2d%%", int(pct+6.f*pctSpan)); // Update progress indicator
        fflush(stdout);
    }
    #endif
    inner->right = Recurse(right, REPORTPRM(pct+6.f*pctSpan) depth+1);
    inner->right->bottom = rbottom;
    inner->right->top = rtop;

    // Used for debugging only - dump the BVH in stdout
    #ifdef DEBUG_LOG_BVH
    float origRange[2], newRangeL[2], newRangeR[2];
    if (bestAxis==0)
    {
        origRange[0] = bottom.X;
        origRange[1] = top.X;
        newRangeL[0] = lbottom.X;
        newRangeL[1] = ltop.X;
        newRangeR[0] = rbottom.X;
        newRangeR[1] = rtop.X;
        }
   else if (bestAxis==1)
    {
        origRange[0] = bottom.Y;
        origRange[1] = top.Y;
        newRangeL[0] = lbottom.Y;
        newRangeL[1] = ltop.Y;
        newRangeR[0] = rbottom.Y;
        newRangeR[1] = rtop.Y;
        }
   else
    {
        origRange[0] = bottom.Z;
        origRange[1] = top.Z;
        newRangeL[0] = lbottom.Z;
        newRangeL[1] = ltop.Z;
        newRangeR[0] = rbottom.Z;
        newRangeR[1] = rtop.Z;
         }
    printf("(%9f,%9f) => (%9f,%9f) and (%9f,%9f)\n", origRange[0], origRange[1],
       newRangeL[0], newRangeL[1], newRangeR[0], newRangeR[1]);
    #endif

    return inner;
}

BVHNode *CreateBVH(const Scene *pScene)//Triangle *Tri
{
    vector<BBoxTmp> work;
    Float3 bottom(FLT_MAX,FLT_MAX,FLT_MAX), top(-FLT_MAX,-FLT_MAX,-FLT_MAX);

    puts("Gathering bounding box info from all triangles...");

    for(unsigned j=0; j<pScene->_triangles.size(); j++)
    {
      const Triangle& tri = pScene->_triangles[j];

      BBoxTmp box;
      box.pTri = &tri;
      box.boxbottom.assignSmaller(*tri._vertexA);
      box.boxbottom.assignSmaller(*tri._vertexB);
      box.boxbottom.assignSmaller(*tri._vertexC);
      box.boxtop.assignBigger(*tri._vertexA);
      box.boxtop.assignBigger(*tri._vertexB);
      box.boxtop.assignBigger(*tri._vertexC);

      bottom.assignSmaller(box.boxbottom);
      top.assignBigger(box.boxtop);

      box.boxcenter = box.boxtop;
      box.boxcenter += box.boxbottom;

      box.boxcenter *= 0.5f;

          work.push_back(box);
        #ifdef DEBUG_LOG_BVH
        printf("ADD: B(%f %f %f) T(%f %f %f) C(%f %f %f)\n",
            box.boxbottom.X,
            box.boxbottom.Y,
            box.boxbottom.Z,
            box.boxtop.X,
            box.boxtop.Y,
            box.boxtop.Z,
            box.boxcenter.X,
            box.boxcenter.Y,
            box.boxcenter.Z);
        #endif
    }

    // ...and pass it to the recursive function that creates the SAH AABB BVH
    // (Surface Area Heuristic, Axis-Aligned Bounding Boxes, Bounding Volume Hierarchy)
    printf("Creating Bounding Volume Hierarchy data...    ");
    fflush(stdout);
    BVHNode *root = Recurse(work);
    cout<< "\b\b\b100%%\n";
    root->bottom = bottom;
    root->top = top;

    return root;
}

int CountBoxes(BVHNode *root)
{
    if (!root->IsLeaf())
    {
      BVHInner *p = dynamic_cast<BVHInner*>(root);
      return 1 + CountBoxes(p->left) + CountBoxes(p->right);
        } else  return 1;
}

unsigned CountTriangles(BVHNode *root)
{
    if (!root->IsLeaf())
    {
       BVHInner *p = dynamic_cast<BVHInner*>(root);
       return CountTriangles(p->left) + CountTriangles(p->right);
        }
    else
    {
       BVHLeaf *p = dynamic_cast<BVHLeaf*>(root);
       return (unsigned) p->_triangles.size();
        }
}

void CountDepth(BVHNode *root, int depth, int& maxDepth)
{
    if (maxDepth<depth)
    maxDepth=depth;
    if (!root->IsLeaf())
    {
      BVHInner *p = dynamic_cast<BVHInner*>(root);
      CountDepth(p->left, depth+1, maxDepth);
      CountDepth(p->right, depth+1, maxDepth);
        }
}

void Scene::PopulateCacheFriendlyBVH(Triangle *pFirstTriangle,
    BVHNode *root, unsigned& idxBoxes,unsigned& idxTriList)
{
    unsigned currIdxBoxes = idxBoxes;
    _pCFBVH[currIdxBoxes]._bottom = root->bottom;
    _pCFBVH[currIdxBoxes]._top    = root->top;

    if (!root->IsLeaf())
    {
      BVHInner *p = dynamic_cast<BVHInner*>(root);
      int idxLeft = ++idxBoxes;
      PopulateCacheFriendlyBVH(pFirstTriangle, p->left, idxBoxes, idxTriList);
      int idxRight = ++idxBoxes;
      PopulateCacheFriendlyBVH(pFirstTriangle, p->right, idxBoxes, idxTriList);
      _pCFBVH[currIdxBoxes].u.inner._idxLeft  = idxLeft;
      _pCFBVH[currIdxBoxes].u.inner._idxRight = idxRight;
        }
    else
    {
      BVHLeaf *p = dynamic_cast<BVHLeaf*>(root);
      unsigned count = (unsigned) p->_triangles.size();
      _pCFBVH[currIdxBoxes].u.leaf._count = 0x80000000 | count;
      _pCFBVH[currIdxBoxes].u.leaf._startIndexInTriIndexList = idxTriList;
      for(std::list<const Triangle*>::iterator it= p->_triangles.begin();
            it != p->_triangles.end();it++)
        {_triIndexList[idxTriList++] = *it - pFirstTriangle;}
        }
}

void Scene::CreateCFBVH()
{
    if (!_pSceneBVH)
    {
      puts("Internal bug in CreateCFBVH, please report it..."); fflush(stdout);
      exit(1);
        }

    unsigned idxTriList=0;
    unsigned idxBoxes=0;

    _triIndexListNo = CountTriangles(_pSceneBVH);
    _triIndexList = new int[_triIndexListNo];

    _pCFBVH_No = CountBoxes(_pSceneBVH);
    _pCFBVH = new CacheFriendlyBVHNode[_pCFBVH_No];

    PopulateCacheFriendlyBVH(&_triangles[0], _pSceneBVH,idxBoxes,idxTriList);

    if ((idxBoxes != _pCFBVH_No-1) || (idxTriList != _triIndexListNo))
    {
      puts("Internal bug in CreateCFBVH, please report it..."); fflush(stdout);
      exit(1);
        }

    int maxDepth = 0;
    CountDepth(_pSceneBVH, 0, maxDepth);
    if (maxDepth>=BVH_STACK_SIZE)
    {
      printf("Max depth of BVH was %d\n", maxDepth);
      puts("Recompile with BVH_STACK_SIZE set to more than that..."); fflush(stdout);
      exit(1);
        }
}

void Scene::UpdateBoundingVolumeHierarchy(const char *filename, bool forceRecalc)
{
    if (!_pCFBVH)
    {
      std::string BVHcacheFilename(filename);
      BVHcacheFilename += ".bvh";
      FILE *fp = fopen(BVHcacheFilename.c_str(), "rb");
      if (forceRecalc || !fp)
       {
         // No cached BVH data - we need to calculate them

        _pSceneBVH = CreateBVH(this);
        printf("Building the BVH%s took %.2f seconds\n");

            // Now that the BVH has been created, copy its data into a more cache-friendly format
       // (CacheFriendlyBVHNode occupies exactly 32 bytes, i.e. a cache-line)

           CreateCFBVH();

      // Now store the results, if possible...
      fp = fopen(BVHcacheFilename.c_str(), "wb");
      if (!fp) return;
      if (1 != fwrite(&_pCFBVH_No, sizeof(unsigned), 1, fp)) return;
      if (1 != fwrite(&_triIndexListNo, sizeof(unsigned), 1, fp)) return;
      if (_pCFBVH_No != fwrite(_pCFBVH, sizeof(CacheFriendlyBVHNode), _pCFBVH_No, fp)) return;
      if (_triIndexListNo != fwrite(_triIndexList, sizeof(int), _triIndexListNo, fp)) return;
      fclose(fp);
    }
    else
    {
      puts("Cache exists, reading the pre-calculated BVH data...");
      if (1 != fread(&_pCFBVH_No, sizeof(unsigned), 1, fp))
        {
           UpdateBoundingVolumeHierarchy(filename, true);
           return;
        }
      if (1 != fread(&_triIndexListNo, sizeof(unsigned), 1, fp))
        {
           UpdateBoundingVolumeHierarchy(filename, true);
           return;
        }
    _pCFBVH = new CacheFriendlyBVHNode[_pCFBVH_No];
    _triIndexList = new int[_triIndexListNo];
    if (_pCFBVH_No != fread(_pCFBVH, sizeof(CacheFriendlyBVHNode), _pCFBVH_No, fp))
        {
          delete [] _pCFBVH;
          delete [] _triIndexList;
          _pCFBVH = NULL;
          _triIndexList = NULL;
          UpdateBoundingVolumeHierarchy(filename, true);
          return;
        }
    if (_triIndexListNo != fread(_triIndexList, sizeof(int), _triIndexListNo, fp))
        {
            delete [] _pCFBVH;
            delete [] _triIndexList;
            _pCFBVH = NULL;
            _triIndexList = NULL;
            UpdateBoundingVolumeHierarchy(filename, true);
            return;
        }
    fclose(fp);
    }
    }
}

void Scene:: fixns()
{
    for(unsigned  j=0; j<_triangles.size(); j++)
    {
       Float3 worldPointA = *_triangles[j]._vertexA;
       Float3 worldPointB = *_triangles[j]._vertexB;
       Float3 worldPointC = *_triangles[j]._vertexC;
       Float3 AB = worldPointB;
       AB -= worldPointA;
       Float3 AC = worldPointC;
       AC -= worldPointA;
       Float3 cr = normalize(Cross(AB, AC));
       _triangles[j].n = cr;

       const_cast<Vertex*>(_triangles[j]._vertexA)->n += cr;
       const_cast<Vertex*>(_triangles[j]._vertexB)->n += cr;
       const_cast<Vertex*>(_triangles[j]._vertexC)->n += cr;

    }
      for(unsigned j=0; j<_triangles.size(); j++)
      {
         const_cast<Vertex*>(_triangles[j]._vertexA)->n.Normalize();
         const_cast<Vertex*>(_triangles[j]._vertexB)->n.Normalize();
         const_cast<Vertex*>(_triangles[j]._vertexC)->n.Normalize();
      }
}

namespace enums
{
    enum ColorComponent {Red = 0, Green = 1, Blue = 2 };
}
using namespace enums;
struct MaterialColors
{

    MaterialColors(
    const Lib3dsRgba& ambient, const Lib3dsRgba& diffuse, const Lib3dsRgba& specular, bool twoSided)
    :
    _ambient(ambient), _diffuse(diffuse), _specular(specular), _twoSided(twoSided)
    {}

    template <ColorComponent c>
    unsigned getByteAmbient() const { return unsigned(255.0*_ambient[c]); }

    template <ColorComponent c>
    unsigned getByteDiffuse() const { return unsigned(255.0*_diffuse[c]); }

    template <ColorComponent c>
    unsigned getByteSpecularRed() const { return unsigned(255.0*_specular[c]); }

    const Lib3dsRgba& _ambient;
    const Lib3dsRgba& _diffuse;
    const Lib3dsRgba& _specular;
    bool _twoSided;
};

const float Scene::MaxCoordAfterRescale = 1.2f;
void Scene::load(const char *filename, Bound *b, Camera *cam)
{
    FILE *fp = fopen(filename, "rb");
    if (filename[0] == '@' && filename[1] == 'p')
    {   // Platform
       _vertices.reserve(4);
       _vertices.push_back(Vertex( 0.5, -0.5, 0.0,  0.0, 0.0, 1.0));
       _vertices.push_back(Vertex( 0.5,  0.5, 0.0,  0.0, 0.0, 1.0));
       _vertices.push_back(Vertex(-0.5,  0.5, 0.0,  0.0, 0.0, 1.0));
       _vertices.push_back(Vertex(-0.5, -0.5, 0.0,  0.0, 0.0, 1.0));
       _triangles.reserve(2);
       _triangles.push_back(Triangle(&_vertices[0], &_vertices[1], &_vertices[2],255,0,0));
       _triangles.push_back(Triangle(&_vertices[0], &_vertices[2], &_vertices[3],255,0,0));
       return;
    }

    const char *dt = strrchr(filename, '.');
    if (dt)
    {
      dt++;
        if (!strcmp(dt, "tri"))
        {

            // Simple binary format:
            //
            // Basically:
            //    <magic (uint32)>  TRI_MAGIC, TRI_MAGICNORMAL or missing
            //    one or more of these blocks:
            //        no_of_vertices (uint32)
            //        x1,y1,z1 (float coordinates of vertex)
            //        ...
            //        xn,yn,zn (float coordinates of vertex)
            //        (magic == TRI_MAGICNORMAL)? nx1, ny1, nz1 (float coordinates of normal)
            //        no_of_triangles (uint32)
            //        idx,idx,idx (uint32 indices into the vertex array)
            //        (magic == TRI_MAGIC | TRI_MAGICNORMAL)? r,g,b (floats)
            //        ...
            //        idx,idx,idx (uint32 indices into the vertex array)
            //        (magic == TRI_MAGIC | TRI_MAGICNORMAL)? r,g,b (floats)
            FILE *fp = fopen(filename, "rb");
            if (!fp)
            cout<<(string("File '") + string(filename) + string("' not found!")).c_str();
            //THROW((string("File '") + string(filename) + string("' not found!")).c_str());
            Uint32 totalPoints = 0, totalTris = 0;
            Uint32 magic;

            fread(&magic, 1, sizeof(Uint32), fp);
            if (magic != TRI_MAGIC && magic != TRI_MAGICNORMAL)
            {
                //// No magic, just vertices and points (no normals, no colors)
                fseek(fp, 0, SEEK_SET);
                cout<<"file not found";
            }

            // Calculate total number of vertices in order to reserve the vectors memory
            unsigned currentOffset = ftell(fp);
            unsigned currentTotalPoints = 0;
            unsigned currentTotalTris = 0;
            while(1)
            {
                unsigned temp;
                fread(&temp, 1, sizeof(Uint32), fp);
                if (feof(fp))
                    break;
                currentTotalPoints += temp;
                fseek(fp, temp*(magic==TRI_MAGICNORMAL?24:12), SEEK_CUR);
                fread(&temp, 1, sizeof(Uint32), fp);
                if (feof(fp))
                    break;
                currentTotalTris += temp;
                fseek(fp, temp*24, SEEK_CUR);
            }

            // Reserve the space, now that you know
            if (magic == TRI_MAGICNORMAL || magic == TRI_MAGIC)
                _vertices.reserve(currentTotalPoints);
            else
            // when we have no normals, we'll need extra space for fixns()...
            _vertices.reserve(currentTotalTris*3);
            _triangles.reserve(currentTotalTris);

            // Now load them inside the std::vectors...
            fseek(fp, currentOffset, SEEK_SET);
            do
            {
                Uint32 noOfPoints;
                fread(&noOfPoints, 1, sizeof(Uint32), fp);
                if (feof(fp))
                    break;
                for(Uint32 i=0; i<noOfPoints; i++)
                {
                    float x,y,z;
                    fread(&x,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    fread(&y,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    fread(&z,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }

                    float nx=0.,ny=0.,nz=0.;
                    if (magic == TRI_MAGICNORMAL) {
                    fread(&nx,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    fread(&ny,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    fread(&nz,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    }
                    else
                    {
                    nx = ny = nz = 0; // Will be calculated in fixns()
                    }
                    assert(_vertices.size() < _vertices.capacity());
                    _vertices.push_back(Vertex(x,y,z, nx,ny,nz));
                }

                Uint32 noOfTris;
                fread(&noOfTris, 1, sizeof(Uint32), fp);
                if (feof(fp)) { cout<<"Malformed 3D file"; }

                for(Uint32 i=0; i<noOfTris; i++)
                {
                    Uint32 idx1,idx2,idx3;
                    fread(&idx1,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    fread(&idx2,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    fread(&idx3,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                    if (idx1>=(totalPoints+noOfPoints)) { cout<<"Malformed 3D file (idx1)"; }
                    if (idx2>=(totalPoints+noOfPoints)) { cout<<"Malformed 3D file (idx2)"; }
                    if (idx3>=(totalPoints+noOfPoints)) { cout<<"Malformed 3D file (idx3)"; }


                    float r, g, b;
                    if (magic == TRI_MAGIC || magic == TRI_MAGICNORMAL)
                        {
                        fread(&r,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                        fread(&g,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                        fread(&b,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                        r*=255.; g*=255.; b*=255.;
                        }
                    else {r = g = b = 255.0;} // No colors? White, then...

                    _triangles.push_back(Triangle(
                        &_vertices[idx1], &_vertices[idx2], &_vertices[idx3],
                        unsigned(r),unsigned(g),unsigned(b)));
                }

                totalPoints += noOfPoints;
                totalTris   += noOfTris;
            } while(!feof(fp));

            fclose(fp);
            if (magic != TRI_MAGICNORMAL)
            fixns();
        }
        else if (!strcmp(dt, "ra2"))
        {
            // http://ompf.org/forum/viewtopic.php?t=64
            FILE *fp = fopen(filename, "rb");
            if (!fp)
            cout<< (string("File '") + string(filename) + string("' not found!")).c_str();

            Uint32 totalPoints = 0, totalTriangles = 0;

            fseek(fp, 0, SEEK_END);
            totalTriangles = ftell(fp)/36;
            totalPoints = 3*totalTriangles;
            fseek(fp, 0, SEEK_SET);

            _vertices.reserve(totalPoints);
            _triangles.reserve(totalTriangles);

            // Now load them inside the std::vectors...
            for(Uint32 i=0; i<totalPoints; i++)
            {
                float x,y,z;
                fread(&y,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                fread(&z,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }
                fread(&x,1,4,fp); if (feof(fp)) { cout<<"Malformed 3D file"; }

                float nx=0.,ny=0.,nz=0.;
                assert(_vertices.size() < _vertices.capacity());
                _vertices.push_back(Vertex(x,y,z, nx,ny,nz));
            }

            for(Uint32 i=0; i<totalTriangles; i++)
            {
            Uint32 idx1,idx2,idx3;
            idx1 = 3*i + 0;
            if (getenv("RA2")!=NULL)
            {
                idx2 = 3*i + 2;
                idx3 = 3*i + 1;
            }
            else
            {
                idx2 = 3*i + 1;
                idx3 = 3*i + 2;
            }
            if (idx1>=(totalPoints)) { cout<<"Malformed 3D file (idx1)"; }
            if (idx2>=(totalPoints)) { cout<<"Malformed 3D file (idx2)"; }
            if (idx3>=(totalPoints)) { cout<<"Malformed 3D file (idx3)"; }

            //uchar4 color;
            float r,g,b;
            r = g = b = 255.0; // No colors? White, then..

            _triangles.push_back(
                Triangle(&_vertices[idx1], &_vertices[idx2], &_vertices[idx3],
                unsigned(r),unsigned(g),unsigned(b)));
            }
            fclose(fp);
            fixns();
        }
        else if (!strcmp(dt, "3ds") || !strcmp(dt, "3DS"))
        {
            int i = 0;
            Lib3dsFile* p3DS = lib3ds_file_load(filename);
            if (!p3DS)
            cout<<"Lib3DS couldn't load this .3ds file";
            lib3ds_file_eval(p3DS, 0);

            Lib3dsDword currentTotalTris = 0;
            Lib3dsMesh *pMesh = p3DS->meshes;
            if (!pMesh)
            cout<<"This .3ds file has no meshes";
            std::map<Lib3dsMesh*, Lib3dsVector*> normals;
            while(pMesh)
            {
                currentTotalTris   += pMesh->faces;
                Lib3dsVector *pMeshNormals = new Lib3dsVector[3*pMesh->faces];
                lib3ds_mesh_calculate_normals(pMesh, pMeshNormals);
                normals[pMesh] = pMeshNormals;
                pMesh = pMesh->next;
            }
            _vertices.reserve(3*currentTotalTris);
            _triangles.reserve(currentTotalTris);
            std::map<string, MaterialColors> colors;
            Lib3dsMaterial *pMaterial = p3DS->materials;
            while(pMaterial)
            {
                colors.insert(
                    std::map<string, MaterialColors>::value_type(
                    string(pMaterial->name),
                    MaterialColors(pMaterial->ambient, pMaterial->diffuse,  pMaterial->specular, pMaterial->two_sided!=0)));
                pMaterial = pMaterial->next;
            }
            pMesh = p3DS->meshes;
            int currentTotalPoints = 0;
            while(pMesh)
            {
            if (!pMesh->pointL) { pMesh=pMesh->next; continue; }
            if (!pMesh->faceL)  { pMesh=pMesh->next; continue; }
            for(i=0; i<(int)pMesh->faces; i++)
            {
                std::map<string, MaterialColors>::iterator pMat = colors.find(string(pMesh->faceL[i].material));
                unsigned r,g,b;
                if (pMat != colors.end())
                {
                    r = pMat->second.getByteDiffuse<Red>();
                    g = pMat->second.getByteDiffuse<Green>();
                    b = pMat->second.getByteDiffuse<Blue>();
                } else {r = g = b = 255; }
                for (int k=0; k<3; k++)
                {
                int pointIdx = pMesh->faceL[i].points[k];
                Float3 nr(
                    normals[pMesh][3*i+k][0],
                    normals[pMesh][3*i+k][1],
                    normals[pMesh][3*i+k][2]);
                 assert(_vertices.size() < _vertices.capacity());
                _vertices.push_back(
                    Vertex(
                           pMesh->pointL[pointIdx].pos[0],
                           pMesh->pointL[pointIdx].pos[1],
                           pMesh->pointL[pointIdx].pos[2],
                    nr.X,
                    nr.Y,
                    nr.Z));
                }
                assert(_triangles.size() < _triangles.capacity());
                _triangles.push_back(
                Triangle(
                    &_vertices[currentTotalPoints + 3*i],
                    &_vertices[currentTotalPoints + 3*i + 1],
                    &_vertices[currentTotalPoints + 3*i + 2],
                    255,0,0,
                    pMat->second._twoSided,
                    true,
                    Float3(pMesh->faceL[i].normal[0],
                        pMesh->faceL[i].normal[1],
                        pMesh->faceL[i].normal[2]) ));
            }//for loop
            currentTotalPoints += 3*pMesh->faces;
            pMesh = pMesh->next;
            }// while ends here
        }
        else if (!strcmp(dt, "PLY") || !strcmp(dt, "ply"))
        {
            // Only shadevis generated objects, not full blown parser!
            std::ifstream file(filename, std::ios::in);
            if (!file) {cout<<(string("Missing ")+string(filename)).c_str();}
            string line;
            unsigned totalVertices, totalTriangles, lineNo=0;
            bool inside = false;
            while(getline(file, line))
            {
            lineNo++;
            if (!inside)
            {
                if (line.substr(0, 14) == "element vertex")
                {
                    std::istringstream str(line);
                    string word1;
                    str >> word1;
                    str >> word1;
                    str >> totalVertices;
                    _vertices.reserve(totalVertices);
                }
                else if (line.substr(0, 12) == "element face")
                {
                    std::istringstream str(line);
                    string word1;
                    str >> word1;
                    str >> word1;
                    str >> totalTriangles;
                }
                else if (line.substr(0, 10) == "end_header")
                inside = true;
            }
            else {
                if (totalVertices)
                {
                    totalVertices--;
                    float x,y,z;
                    x = 0.0f; y = 0.0f; z = 0.0f;
                    unsigned ambientOcclusionCoeff;
                    std::istringstream str(line);
                    str >> x >> y >> z >> ambientOcclusionCoeff;
                    _vertices.push_back(Vertex(x,y,z,0.0f,0.0f,0.0f,ambientOcclusionCoeff));
                }
                else if (totalTriangles)
                {
                    totalTriangles--;
                    //uchar4 Color_totalTriangles;
                    unsigned dummy,r,g,b;
                    unsigned idx1, idx2, idx3;
                    std::istringstream str(line);
                    if (str >> dummy >> idx1 >> idx2 >> idx3)
                    {
                        if (str >> r >> g >> b)
                        ;
                        else {r =  g =  b = 255;}
                        _triangles.push_back(
                        Triangle( &_vertices[idx1], &_vertices[idx2], &_vertices[idx3],
                            r,g,b));
                    }
                }
            }
            }
            fixns();
        }
        else
            cout<<"Unknown extension (only .tri .3ds or .ply accepted)";
    }
    else
    cout<<"No extension in filename (only tri .3ds or .ply accepted)";

    std::cout << "Vertexes: " << _vertices.size();
    std::cout << " Triangles: " << _triangles.size() << std::endl;

    // Center scene at world's center

    Float3 minp(FLT_MAX,FLT_MAX,FLT_MAX), maxp(-FLT_MAX,-FLT_MAX,-FLT_MAX);

    for(unsigned i=0; i<_triangles.size(); i++)
    {
        minp.assignSmaller(*_triangles[i]._vertexA);
        minp.assignSmaller(*_triangles[i]._vertexB);
        minp.assignSmaller(*_triangles[i]._vertexC);
        maxp.assignBigger(*_triangles[i]._vertexA);
        maxp.assignBigger(*_triangles[i]._vertexB);
        maxp.assignBigger(*_triangles[i]._vertexC);
    }

    Float3 origCenter = Float3( (maxp.X+ minp.X)/2,
                                (maxp.Y+minp.Y)/2,
                                (maxp.Z+minp.Z)/2);

    b->minBound = minp;
    b->maxBound = maxp;

    minp -= origCenter;
    maxp -= origCenter;

    b->center = origCenter;

    float dX = (b->maxBound.X-b->minBound.X);
    float dY = (b->maxBound.Y-b->minBound.Y);
    float dZ = (b->maxBound.Z-b->minBound.Z);
    b->dia = sqrt(dX*dX+dY*dY+dZ*dZ);

    cam->At = b->center;
    if ((dX > dY) && (dX > dZ))
    {
        cam->Up.X = 1.0f; cam->Up.Y = cam->Up.Z = 0.0f;
        if (dY > dZ)
        {
            cam->From.Z = b->center.Z+b->dia;
            cam->From.X = b->center.X;cam->From.Y = b->center.Y;
        }
        else
        {
            cam->From.Y = b->center.Y+b->dia;
            cam->From.X = b->center.X;cam->From.Z = b->center.Z;
        }
    }
    else if((dY > dX) && (dY > dZ))
    {
            cam->Up.Y = 1.0f; cam->Up.X = cam->Up.Z = 0.0f;
        if (dX > dZ)
        {
            cam->From.Z = b->center.Z+b->dia;
            cam->From.Y = b->center.Y;cam->From.X = b->center.X;
        }
        else
        {
            cam->From.X = b->center.X+b->dia;
            cam->From.Y = b->center.Y;cam->From.Z = b->center.Z;
        }
    }
    else
    {
        cam->Up.Z = 1.0f; cam->Up.X = cam->Up.Y = 0.0f;
        if (dX > dY)
        {
            cam->From.Y = b->center.Y+b->dia;
            cam->From.X = b->center.X;cam->From.Z = b->center.Z;
        }
        else
        {
            cam->From.X = b->center.X-b->dia;
            cam->From.Y = b->center.Y;cam->From.Z = b->center.Z;
        }
    }
    // Scale scene so max(abs x,y,z coordinates) = MaxCoordAfterRescale
    //typedef float coord;

    float maxi = 0.0f;
    maxi = max(maxi, (float) fabs(minp.X));
    maxi = max(maxi, (float) fabs(minp.Y));
    maxi = max(maxi, (float) fabs(minp.Z));
    maxi = max(maxi, (float) fabs(maxp.X));
    maxi = max(maxi, (float) fabs(maxp.Y));
    maxi = max(maxi, (float) fabs(maxp.Z));

    for(unsigned i=0; i<_vertices.size(); i++)
    {
        _vertices[i] -= origCenter;
        _vertices[i] *= MaxCoordAfterRescale/maxi;
    }
    for(unsigned i=0; i<_triangles.size(); i++)
    {
        _triangles[i].centroid -= origCenter;
        _triangles[i].centroid *= MaxCoordAfterRescale/maxi;
    }
    // Update triangle bounding boxes (used by BVH)
    for(unsigned i=0; i<_triangles.size(); i++)
    {
        _triangles[i]._bottom.assignSmaller(*_triangles[i]._vertexA);
        _triangles[i]._bottom.assignSmaller(*_triangles[i]._vertexB);
        _triangles[i]._bottom.assignSmaller(*_triangles[i]._vertexC);
        _triangles[i]._top.assignBigger(*_triangles[i]._vertexA);
        _triangles[i]._top.assignBigger(*_triangles[i]._vertexB);
        _triangles[i]._top.assignBigger(*_triangles[i]._vertexC);
    }
    // Pre-compute triangle intersection data (used by raytracer)

    for(unsigned i=0; i<_triangles.size(); i++)
    {
        Triangle& triangle = _triangles[i];

        // Algorithm for triangle intersection is taken from Roman Kuchkuda's paper.
        // edge vectors
        Float3 vc1=*triangle._vertexB; vc1-=*triangle._vertexA;
        Float3 vc2=*triangle._vertexC; vc2-=*triangle._vertexB;
        Float3 vc3=*triangle._vertexA; vc3-=*triangle._vertexC;

        // plane of triangle
        triangle.n = Cross(vc1, vc2);
        Float3 alt1 = Cross(vc2, vc3);
        if (alt1.length() > triangle.n.length()) triangle.n = alt1;
        Float3 alt2 = Cross(vc3, vc1);
        if (alt2.length() > triangle.n.length()) triangle.n = alt2;
        triangle.n.Normalize();
        triangle._d = Dot(triangle.n, *triangle._vertexA);

        // edge planes
        triangle._e1 = Cross(triangle.n, vc1);
        triangle._e1.Normalize();
        triangle._d1 = Dot(triangle._e1, *triangle._vertexA);
        triangle._e2 = Cross(triangle.n, vc2);
        triangle._e2.Normalize();
        triangle._d2 = Dot(triangle._e2, *triangle._vertexB);
        triangle._e3 = Cross(triangle.n, vc3);
        triangle._e3.Normalize();
        triangle._d3 = Dot(triangle._e3, *triangle._vertexC);
    }
}

bool RayIntersectsBox(const Float3& originInWorldSpace, const Float3& rayInWorldSpace, CacheFriendlyBVHNode *pBox)
{
    // set Tnear = - infinity, Tfar = infinity
    //
    // For each pair of planes P associated with X, Y, and Z do:
    //     (example using X planes)
    //     if direction Xd = 0 then the ray is parallel to the X planes, so
    //         if origin Xo is not between the slabs ( Xo < Xl or Xo > Xh) then
    //             return false
    //     else, if the ray is not parallel to the plane then
    //     begin
    //         compute the intersection distance of the planes
    //         T1 = (Xl - Xo) / Xd
    //         T2 = (Xh - Xo) / Xd
    //         If T1 > T2 swap (T1, T2) /* since T1 intersection with near plane */
    //         If T1 > Tnear set Tnear =T1 /* want largest Tnear */
    //         If T2 < Tfar set Tfar="T2" /* want smallest Tfar */
    //         If Tnear > Tfar box is missed so
    //             return false
    //         If Tfar < 0 box is behind ray
    //             return false
    //     end
    // end of for loop
    //
    // If Box survived all above tests, return true with intersection point Tnear and exit point Tfar.

    CacheFriendlyBVHNode& box = *pBox;
    float Tnear, Tfar;
    Tnear = -FLT_MAX;
    Tfar = FLT_MAX;

    #define CHECK_NEAR_AND_FAR_INTERSECTION(c)                                          \
    if (rayInWorldSpace.## c == 0.){                                    \
        if (originInWorldSpace.##c < box._bottom.##c) return false;                 \
        if (originInWorldSpace.##c > box._top.##c) return false;                    \
    } else{                                            \
        float T1 = (box._bottom.##c - originInWorldSpace.##c)/rayInWorldSpace.##c;          \
        float T2 = (box._top.##c    - originInWorldSpace.##c)/rayInWorldSpace.##c;          \
        if (T1>T2) {float tmp=T1; T1=T2; T2=tmp; }                          \
        if (T1 > Tnear) Tnear = T1;                                 \
        if (T2 < Tfar)  Tfar = T2;                                  \
        if (Tnear > Tfar)                                       \
            return false;                                       \
        if (Tfar < 0.)                                          \
            return false;                                       \
        }

    CHECK_NEAR_AND_FAR_INTERSECTION(X)
    CHECK_NEAR_AND_FAR_INTERSECTION(Y)
    CHECK_NEAR_AND_FAR_INTERSECTION(Z)

    return true;
}

////////////////Raytracer////////////
class RaytraceScanline
{
    // Since this class contains only references and has no virtual methods, it (hopefully)
    // doesn't exist in runtime; it is optimized away when RaytraceHorizontalSegment is called.
   // const Triangle& scene;
    const Scene& scene;
    const Camera& eye;
    Screen& canvas;
    int& y;
    bool& antialias;
public: RaytraceScanline(const Scene& scene,
                    //const Triangle& scene,
                    const Camera& e, Screen& c,
                    int& scanline,bool withAntialiasing)
                    :
                    scene(scene), eye(e),
                    canvas(c), y(scanline),
                    antialias(withAntialiasing)
                    {}

template <bool stopAtfirstRayHit, bool doCulling>
bool BVH_IntersectTriangles(
     //Inputs
    const Float3& origin, const Float3& ray,
    const Triangle *avoidSelf,
     //outputs
    const Triangle*& pBestTri,


    // both inputs and outputs!

     //for normal rays:
      //pointHitInWorldSpace (output)
      //    kXX (outputs) perpendicular distances of intersection point from the 3 triangle edges
      //   (used for PhongNormal calculations)

     //for shadow rays:
      //pointHitInWorldSpace (input) provides the light position
    Float3& pointHitInWorldSpace,
    float& kAB, float& kBC, float& kCA
)const{
     //in the loop below, maintain the closest triangle and the point where we hit it:
        pBestTri = NULL;
        float bestTriDist;

        // light position passed-in pointHitInWorldSpace (only in shadow mode - i.e. stopAtfirstRayHit=true)
        Float3& lightPos = pointHitInWorldSpace;
        bool stopAtfirstRayHit = true;
         //Compile-time work (stopAtfirstRayHit is template param)
        if (stopAtfirstRayHit)
            // In shadow ray mode, start from light distance
            bestTriDist = distancesq(origin, lightPos);
        else
            // In normal mode, start from infinity
            bestTriDist = FLT_MAX;

        CacheFriendlyBVHNode* stack[BVH_STACK_SIZE];
        //BVHNode* stack[BVH_STACK_SIZE];
        //bvh box index
        int stackIdx = 0;
        stack[stackIdx++] = scene._pCFBVH;
        while(stackIdx)
        {
            CacheFriendlyBVHNode *pCurrent = stack[stackIdx-1];
            stackIdx--;
            //if (!pCurrent->IsLeaf()) {
            if (!(pCurrent->u.leaf._count & 0x80000000))
            {
                if (RayIntersectsBox(origin,ray, pCurrent))
                {
                    stack[stackIdx++] = &scene._pCFBVH[pCurrent->u.inner._idxRight];
                    stack[stackIdx++] = &scene._pCFBVH[pCurrent->u.inner._idxLeft];
                    assert(stackIdx<=BVH_STACK_SIZE);
                }
            }
            else
            {
                //BVHLeaf *p = dynamic_cast<BVHLeaf*>(pCurrent);
                //for(std::list<const Triangle*>::iterator it=p->_triangles.begin();
                //    it != p->_triangles.end();
                //    it++)
                for(unsigned i=pCurrent->u.leaf._startIndexInTriIndexList;
                    i<pCurrent->u.leaf._startIndexInTriIndexList +
                    (pCurrent->u.leaf._count & 0x7fffffff); i++)
                   {
                        //const triangle& triangle = *(*it);
                        const Triangle& Tri = scene._triangles[scene._triIndexList[i]];

                        if (avoidSelf == &Tri)continue; // avoid self-reflections/refractions

                         //doCulling is a compile-time param, this code will be "codegenerated"
                        // at compile time only for reflection-related calls to Raytrace (see below)

                        if (doCulling && !Tri._twoSided)
                        {
                             //Check visibility of triangle via dot product
                            Float3 fromTriToOrigin = origin;
                            fromTriToOrigin -= Tri.centroid;
                            // Normally we would normalize, but since we just need the sign
                            //of the dot product (to determine if it facing us or not)...
                            fromTriToOrigin = normalize(fromTriToOrigin);
                            if (Dot(fromTriToOrigin, Tri.n)<0)continue;
                        }

                        // Use the pre-computed triangle intersection data: normal
                        float k = Dot(Tri.n, ray);
                        if (k == 0.0f)
                        continue; // this triangle is parallel to the ray, ignore it.

                        float s = (Tri._vertexA->X - Dot(Tri.n, origin))/k;

                        if(s <= 0.0f && s <=0.0f && s <=0.0f) // this triangle is "behind" the origin.
                         continue;
                        if (s <= NUDGE_FACTOR && s <= NUDGE_FACTOR && s <= NUDGE_FACTOR )
                        continue;
                        if (s <= 0.0) // this triangle is "behind" the origin.
                        continue;
                        if (s <= NUDGE_FACTOR)
                        continue;

                        Float3 hit = Scale(ray,s);
                        hit -= origin;

                        //Is the intersection of the ray with the triangle's plane INSIDE the triangle?
                        float kt1 = Dot(Tri._e1, hit) - Tri._d1; if (kt1<0.0) continue;
                        float kt2 = Dot(Tri._e2, hit) - Tri._d2; if (kt2<0.0) continue;
                        float kt3 = Dot(Tri._e3, hit) - Tri._d3; if (kt3<0.0) continue;

                         //It is, "hit" is the world space floatinate of the intersection.

                         //check whether it Was this a normal ray or a shadow ray? (template param)
                        if (stopAtfirstRayHit)
                        {
                            // Shadow ray, check whether the triangle obstructs the light
                            float dist = distancesq(lightPos, hit);
                            if (dist < bestTriDist) // distance to light (squared) passed in kAB
                                return true; // we found a triangle obstructing the light, return true
                        }
                        else
                        { // Normal ray - it this intersection closer than all the others?
                            float hitZ = distancesq(origin, hit);
                            if (hitZ < bestTriDist)
                            {
                                // maintain the closest hit
                                bestTriDist = hitZ;
                                pBestTri = &Tri;
                                pointHitInWorldSpace = hit;
                                kAB = kt1;
                                kBC = kt2;
                                kCA = kt3;
                            }
                        }
                    }
                }
        }
         //Normal ray or shadow ray? (compile-time template param)
        if (!stopAtfirstRayHit)
           //  for normal ray, return true if we pierced a triangle
              return pBestTri != NULL;
        else
            // for shadow ray, return true if we found a triangle obstructing the light.
            return false;
}
template <bool doCulling>
Pixel Raytrace(Float3 originInWorldSpace,
               Float3 rayInWorldSpace,
               const Triangle *avoidSelf, int depth)const
{
    if (depth >= MAX_RAY_DEPTH)
        return Pixel(0.,0.,0.);

    const Triangle *pBestTri = NULL;
    Float3 pointHitInWorldSpace;
    float kAB= 0.f, kBC= 0.f, kCA= 0.f; // distances from the 3 edges of the triangle (from where we hit it)

    // Use the surface-area heuristic based, bounding volume hierarchy of axis-aligned bounding boxes
    // (keywords: SAH, BVH, AABB)
    if (!BVH_IntersectTriangles<false,doCulling>(
        originInWorldSpace, rayInWorldSpace, avoidSelf,
        pBestTri, pointHitInWorldSpace, kAB, kBC, kCA))
        // We pierced no triangle, return with no contribution (ambient is black)
        return Pixel(0.,0.,0.);

    // Set this to pass to recursive calls below, so that we don't get self-shadow or self-reflection
    // from this triangle...
    avoidSelf = pBestTri;

    // We'll also calculate the color contributed from this intersection
    // Start from the triangle's color
    Pixel color = pBestTri->_colorf;

    return color;
}
void RaytraceHorizontalSegment(int xStarting, int iOnePastEndingX) const
{
    for(int x=xStarting; x<iOnePastEndingX; x++)
    {
        Pixel finalColor(0,0,0);

        int pixelsTraced = 1;
        if (antialias)
        pixelsTraced = 4;

        while(pixelsTraced--)
        {
            // We will shoot a ray in camera space (from Eye to the screen point, so in camera
            // space, from (0,0,0) to this:
            float xx = (float)x;
            float yy = (float)y;

            if (antialias)
            {
                // nudge in a cross pattern around the uchar4 center
                xx += 0.25f - .5f*(pixelsTraced&1);
                yy += 0.25f - .5f*((pixelsTraced&2)>>1);
            }
            float lx = float((HEIGHT/2)-yy)/SCREEN_DIST;
            float ly = float(xx-(WIDTH/2))/SCREEN_DIST;
            float lz = 1.0;
            Float3 rayInCameraSpace(lx,ly,lz);
            rayInCameraSpace.Normalize();

            // We will need the origin in world space
            Float3 originInWorldSpace = eye.From;

            // We have a rayInCameraSpace, and we want to use the BVH, which was constructed
            // in World space, so we convert the ray in World space
            Float3 rayInWorldSpace = Scale(eye.From, rayInCameraSpace.X);
            rayInWorldSpace += Scale(eye.At, rayInCameraSpace.Y);
            rayInWorldSpace += Scale(eye.Up, rayInCameraSpace.Z);
            // in theory, this should not be required
            rayInWorldSpace.Normalize();

            // Primary ray, we want backface culling: <true>
            finalColor += Raytrace<true>(originInWorldSpace, rayInWorldSpace, NULL, 0);
        }
       // if (antialias)
        finalColor /= 4.;
        if (finalColor._r>255.0f) finalColor._r=255.0f;
        if (finalColor._g>255.0f) finalColor._g=255.0f;
        if (finalColor._b>255.0f) finalColor._b=255.0f;
        canvas.DrawPixel(y,x, SDL_MapRGB(
        canvas._surface->format, (Uint8)finalColor._r, (Uint8)finalColor._g, (Uint8)finalColor._b));
    }
    }
};

bool Scene :: renderRaytracer(Camera& eye, Screen& canvas, bool antialias)
{
    bool needToUpdateTitleBar = !_pSceneBVH; // see below
    // Update the BVH and its cache-friendly version
     extern const char *g_filename;
    UpdateBoundingVolumeHierarchy(g_filename);

    // Main loop: for each uchar4...
    for(int y=0; y<HEIGHT; y++)
    {

    // For both OpenMP and single-threaded, call the RaytraceHorizontalSegment member
    // of the RaytraceScanline, requesting drawing of ALL the scanline.
    // For OpenMP, the appropriate pragma inside RaytraceHorizontalSegment will make it execute via SMP...
        RaytraceScanline(*this, eye, canvas, y, antialias).
        RaytraceHorizontalSegment(0, WIDTH);

    }
    canvas.ShowScreen(true);
    return true;
}

bool g_benchmark = false;
const char *g_filename = NULL;

 int main(int argc,char *argv[])
{
    char *fname = argv[1];
    g_filename = fname;
    bool doBenchmark = false;
    g_benchmark =doBenchmark;

    ComputeCameraframe();
    try
    {
        Scene scene;
        scene.load(fname, &b, &Cam);
            if (g_benchmark )
            {
                // When benchmarking, we dont want the first frame to "suffer" the BVH creation
                puts("Creating BVH... please wait...");
                scene.UpdateBoundingVolumeHierarchy(fname);
            }

        Screen canvas(scene);
        canvas.ShowScreen();
        scene.renderRaytracer(Cam,canvas,false);
    }
    catch(const string& s) {cerr << s.c_str();}
    return 0;
 }