• fullscreen
  • Camera.pde
  • bvh.pde
  • bvhBuilder.pde
  • cache.pde
  • objLoader.pde
  • sampling.pde
  • traverse.pde
  • workQueue.pde
  • // camera
    
    class PolarCamera
    {
      float  m_camD;
      float   m_CamAngleX = -0.4f;
      float   m_CamAngleY = -0.5f;
      PVector  m_camDir;
      float    m_camDelta;
      PVector m_camPos;
      PolarCamera()
      {
        update(0,0);
      }
      String Settings()
      {
        return "new PolarCamera( new PVector(" + m_camPos.x +", "
              +  m_camPos.y +", " + m_camPos.z + "), "
              + m_CamAngleX + ", " +  m_CamAngleY + ") ";
      }
      PolarCamera( float d, float cax, float cay, float cd )
      {
        m_camD= d;
        m_CamAngleX = cax;
        m_CamAngleY = cay;
         m_camDelta = cd;
        update( 0,0);
      }
      void update(float dy, float dz)
      {
         m_CamAngleX += dy;
         m_CamAngleY += -dz;
        rotate(m_CamAngleX, m_CamAngleY);              
      }
      void forward(float v )
      {
        m_camD +=v*0.25* m_camDelta;
         rotate(m_CamAngleX, m_CamAngleY);              
      }
       void rotate(float ang, float angy )
       {
         m_camDir = new PVector();
         m_camDir.x = sin(ang)*cos(angy);
         m_camDir.z = cos(ang)*cos(angy);
         m_camDir.y = sin(angy);
         m_camPos = new PVector(0.0,0.0f,0.0);
         m_camPos.x = -m_camDir.x*m_camD;
         m_camPos.z = -m_camDir.z*m_camD;     
         m_camPos.y =-m_camDir.y*m_camD +.25f;
         m_camDir.normalize();
         
       }
      
      void side(float v )
      {
        m_camD +=(v*0.25* m_camDelta);
       rotate(m_CamAngleX, m_CamAngleY);              
      }
    }
    PVector lerp( PVector a, PVector b, float t) {
      return new PVector( lerp(a.x,b.x,t),lerp(a.y,b.y,t),lerp(a.z,b.z,t));
    }
    class Ray
    {
      public PVector origin;
      public PVector direction;
      Ray() {
        origin =new PVector();
        direction = new PVector();
      }
      Ray( PVector o, PVector d)
      {
        origin=o;
        direction=d;
      }
      void Lerp( Ray or, float t)
      {
        origin = lerp( origin, or.origin,t);
        direction = lerp( direction, or.direction,t);
        direction.normalize();
      }
    }
    
    class Camera
    {
      PVector position;
      PVector lower_left_corner;
      PVector uvw_u;
      PVector uvw_v;
      PVector uvw_w;
    
      public Camera(PVector eye, PVector up, PVector gaze,
      float resX, float resY, float d)
      {
        position = eye;
        uvw_v = up;
        uvw_w = gaze;
        uvw_u = uvw_v.cross(uvw_w);
        uvw_u.normalize();
        
        uvw_v = uvw_w.cross(uvw_u );
        uvw_v.normalize();
        
        uvw_w.normalize();
        
        float sX = -resX/2;
        float sY = -resY/2;
        float eX = resX/2;
        float eY = resY/2;
        
        lower_left_corner = new PVector();
        PVector t1= new PVector();
        PVector t2= new PVector();
        t1.set( uvw_w);
        t1.mult(d);
        t2.set(uvw_v);
        t2.mult(sY);
        lower_left_corner.set( uvw_u);
        lower_left_corner.mult(sX);
        lower_left_corner.add(t1);
        lower_left_corner.add(t2);
            
        //uvw_u.mult(resX);
        //uvw_v.mult( resY);
      }
    
      void generateRay( Ray ray, float screenCoordX, float screenCoordY )
      {
        ray.direction.set( lower_left_corner );
        ray.direction.add( PVector.mult(uvw_v,screenCoordY));
        ray.direction.add(PVector.mult(uvw_u,screenCoordX));
        ray.direction.normalize();
        ray.origin.set(position);
        return;
      }
      Ray generateRay( float screenCoordX, float screenCoordY )
      {
        PVector origin =new PVector();
        PVector up =new PVector();
        PVector acc =new PVector();
        PVector direction = new PVector();
        direction.set( lower_left_corner );
        up.set( uvw_v);
        up.mult( screenCoordY);
        acc.set(uvw_u);
        acc.mult( screenCoordX);
        direction.add(up);
        direction.add(acc);
        direction.normalize();
        origin.set(position);
        return new Ray( origin, direction);
      }
      float[] generateScreen( Ray ray)
      {
        float[] r = new float[2];
        r[0]=-1.;
        r[1]=-1.;
        
        // find intersection with screen plane
        /*float vd = ray.direction.dot(uvw_w);
        if ( abs(vd)<0.0001)
          return r;
          
        float v0 = -(ray.origin.dot(n)+d);
        float t =  v0/vd;
        if ( t < 0.)
          return r;
        */  
        // find hitpoint andreturn it
        return r;
      }
    };
    
    
    // BVH - maybe separate sample
    
    TraversalTree ttree;
    
    class RenderState
    {
      float[] pdf;  
      float[] cdf;    
      float[] samples2;
      hitResult hit;
      int[] traversalStack;
      RenderState( TraversalTree ttree)
      {
        hit = new hitResult();
        traversalStack = new int[32];
        pdf = new float[ttree.luminares.length];
        cdf = new float[ttree.luminares.length];
        samples2 = new float[2];
      }
    }
    
    // http://www.cs.caltech.edu/~njlitke/meshes/toc.html
    class Scene
    {
      String n;
      String des;
      float zflip;
      float gplane;
      float rescale;
      boolean addLight;
      int     numRequireSamples;
      float  cangX;
      int    colorSel;
      boolean useReflect;
      
      Scene( String _n, String _d, float _z, float _g) {
        n =_n;
        zflip=_z;
        gplane = _g;
        rescale = 1.f;
        addLight =true;
        numRequireSamples=64;
        des=_d;
        cangX=0.;
        colorSel=-1;
        useReflect=true;
      }
      Scene( String _n, String _d, float _z, float _g, float _r, boolean _a, int rs, float cax, boolean _ur, int _cs) {
        n =_n;
        zflip=_z;
        gplane = _g;
        rescale=_r;
        addLight =_a;
        numRequireSamples=rs;
        des=_d;
        cangX=cax;
        colorSel=_cs;
        useReflect=_ur;
      }
    };
    int g_NumSamplesRequired = 128;
    int g_CurrentScene = 0;
    
    Scene g_Scenes[]=new Scene[] {
       new Scene( "killerroo.obj",
                 "Model property of headus (metamorphosis) Pty Ltd",0.,-0.73,1.6,true,64,-0.5,true,6),
      
      /*  new Scene( "bottles.obj",
                 "Model property of headus (metamorphosis) Pty Ltd",0.,-100.,0.6,true,64,-0.5,true,6),
      
       new Scene( "killeroo-25.obj",
                 "Model property of headus (metamorphosis) Pty Ltd",0.,-0.7,1.6,true,64,-0.5,true,6),
        */new Scene( "teapot.obj",
                  "Model from The Mesh Compendium", 0.,-.5),
        new Scene( "tRex.obj",
                  "Model by Joel Anderson (joel3d.com )",0.,0.,0.8,true,128,-3.0,false,0),
       new Scene( "knot.obj",
                 "Model from The Mesh Compendium",0.,-1.,1.0,true,128,0.,true,0),
      new Scene( "spock.obj",
                 "Model courtesy of the University of Washington", 0.,-0.825 - 10.,1.,true,64,1.5,true,1),
      new Scene( "triceratops.obj",
                  "The triceratops model is the property of Viewpoint",0.,-0.7,1.6,true,64,-0.5,true,8),
      new Scene("geometry.obj",
                "Model from the Global Illumination Test Scenes V 0.9", 0.,-0.825 - 10., 0.01, false,100000,0.,false,0), 
      new Scene( "secondary.obj",
                  "Model from the Global Illumination Test Scenes V 0.9", 0.,-0.825 - 10., 0.01, false,100000,0.,false,0), 
      new Scene( "shadow.obj",
                  "Model from the Global Illumination Test Scenes V 0.9", 0.,-0.825 - 10., 0.01, false,64,0.,false,0),
    };
    WorkQueue  g_Workers;
    int g_FrameStart;
    void setup()
    {
      size(384,384);
      int numThreads = Runtime.getRuntime().availableProcessors();
      println("Using "+ numThreads+ " Threads");
      g_Workers = new WorkQueue(numThreads);
      frameRate(60);
      LoadScene( g_Scenes[g_CurrentScene%g_Scenes.length]);
      frameBuffer = new PVector[width*height];
      for (int i=0; i<width*height;i++)
        frameBuffer[i] = new PVector();
    
      rendered = createImage( width, height, RGB);
      g_FrameStart=millis();
    }
    
    PolarCamera g_camControl = new PolarCamera();
    
    int hash_rand(int a, int b, int c)
    {
      a -= b;
      a -= c;
      a ^= (c>>13);
      b -= c;
      b -= a;
      b ^= (a<<8);
      c -= a;
      c -= b;
      c ^= (b>>13);
      a -= b;
      a -= c;
      a ^= (c>>12);
      b -= c;
      b -= a;
      b ^= (a<<16);
      c -= a;
      c -= b;
      c ^= (b>>5);
      a -= b;
      a -= c;
      a ^= (c>>3);
      b -= c;
      b -= a;
      b ^= (a<<10);
      c -= a;
      c -= b;
      c ^= (b>>15);
      return c;
    }
    void resetAccumulationBuffer()
    {
      g_stationaryFrame=0;
    }
    void mouseDragged() {
    
      float dx = (float)(mouseX-pmouseX)/width;
      float dy = (float)(mouseY-pmouseY)/height;
      if (keyPressed && keyCode == SHIFT)
      {
        g_camControl.forward(dy*4.);
        g_camControl.side(dx*4.);
      }
      else if ( keyPressed && keyCode == CONTROL)
      {
        g_lw = max(g_lw+dy,0.01);
        g_lh = max(g_lh+dx,0.01);
        AddLight();
        g_applyFlush=true; // flush the cache
      }
      else if ( keyPressed && keyCode == ALT)
      {
        g_langX = g_langX+dy;
        g_langY = g_langY+dx;
        AddLight();
        g_applyFlush=true; // flush the cache
      }
      else
        g_camControl.update(dx,dy);
    
      resetAccumulationBuffer();
    }
    
    boolean g_DrawDirect = true;
    boolean g_useReflect =true;
    boolean g_useCacheOnly =false;
    int  g_colorSelection=4;
    boolean g_applyFlush = false;
    
      boolean g_Animate = false;
    // Real-Time Path Tracer
    // Loads a series of obj models ( based from ahmet.kizilay obj loader)
    // and renders then using path tracing. 
    // Path tracer uses progressive spatial caching to get real-time 
    // update and views of model.
    // Keys
    // '+' go to next model.
    // Drag mouse to view around scene
    // Drag and hold shift to zoom/ pan
    // Drag and hold control to resize area light
    // Drag and hold alt to move area light.
    // 'D' toggle direct lighting.
    // 'R' toggle reflection
    // 'C' show spatial caching 
    // 'S' change model colors
    void keyPressed()
    {
      if ( key=='+' || key =='=')
      {
        g_CurrentScene++;
        LoadScene( g_Scenes[g_CurrentScene%g_Scenes.length]);
        g_applyFlush=true; // flush the cache
      }
      if (key=='D' || key=='d')
      {
        g_DrawDirect =!g_DrawDirect;
        g_applyFlush=true; // flush the cache
      }
      if (key=='R' || key=='r')
        g_useReflect=!g_useReflect;
      if (key=='C' || key=='c')
        g_useCacheOnly=!g_useCacheOnly;
    
      if (key=='S' || key=='s') {
        g_colorSelection++;
        g_applyFlush=true; // flush the cache
      }
      if ( key=='A' || key=='a'){
        g_Animate = !g_Animate;
      }
      resetAccumulationBuffer();
    }
    PVector  GenerateDirectionVector( PVector direction, int sc, QMCSamples sampler,
    PVector nrm, int D)
    {
      PVector tg = new PVector();
      if ( abs(nrm.y)<0.99f)
      {
        tg = new PVector(0.f,1.f,0.f);
      }
      else
      {
        tg = new PVector(1.f,0.f,0.f);
      }
      PVector bi = tg.cross(nrm);
      bi.normalize();
      tg = nrm.cross(bi);
      float[] hemiSample = sampler.Get(sc, D);
      PVector hemi = to_unit_disk( hemiSample[0], hemiSample[1],1.);
      PVector hdir = new PVector();
      hdir.x = tg.x*hemi.x + nrm.x*hemi.y + bi.x*hemi.z;
      hdir.y = tg.y*hemi.x + nrm.y*hemi.y + bi.y*hemi.z;
      hdir.z = tg.z*hemi.x + nrm.z*hemi.y + bi.z*hemi.z;
      direction=hdir;
    
      return direction;
    }
    
    
    RenderCache g_RenderCache = new RenderCache();
    int frame =0;
    
    PVector  GetDiffuse( PVector pt, PVector objCol, PVector grdCol)
    {
      return (pt.y> g_Scenes[g_CurrentScene%g_Scenes.length].gplane+0.001) ?objCol : grdCol;
    }
    int CalcNoiseOctaves( float pixSize, float s, float octSize )
    {
      pixSize = pixSize*octSize ; 
      float si = 1/s;
      int numOctaves = 0;
      while( si > pixSize)
      {
        si*=2.;
        numOctaves++;
        assert(numOctaves < 30);
      }
      return numOctaves;
    }
    int g_stationaryFrame=0;
    float Fresnel( PVector nrm, PVector I)
    {
      float f = 1.-abs(nrm.dot(I));
      float f2 = f*f;
      f2 = f2*f2 * f;
      return 0.1+ 0.9*f2;
    }
    PVector Reflect( PVector nrm, PVector I) {
      float f = -2.* nrm.dot(I);
      PVector r = new PVector();
      r.set(nrm);
      r.mult(f);
      r.add(I);
      return r;
    }
    PVector CalcReflectionColor( hitResult ht,RenderState rs, Ray ray,float [] origin, float [] dv2, 
    PVector hpt, PVector nrm, PVector lgtCol )
    {
      PVector reflectCol=new PVector();
      PVector reflectDir = Reflect(  nrm,ray.direction);
      float maxT = 10.f;
    
      dv2[0]= reflectDir.x;
      dv2[1]=reflectDir.y;
      dv2[2]=reflectDir.z;
      ht = ttree.Trace(rs, origin, dv2, maxT );
      if ( ht.t < maxT )
      {       
        PVector htpt= new PVector();  
        htpt.set(PVector.add(PVector.mult( reflectDir,ht.t),ray.origin));
        float ds =1./0.02;;
        int lod = 4;
        ds *= (float)(5-lod);
        float off = ((float)lod-1.)*10.+ht.tri;
    
        int ix = (int)(htpt.x*ds);
        int iy= (int)((htpt.y+off)*ds);
        int iz= (int)(htpt.z*ds);
        int key1 = hash_rand(ix,iy,iz);
        int cacheIdx = g_RenderCache.FindIdx(key1);
        g_RenderCache.Get( cacheIdx, key1, reflectCol, frame);
      }
      else
      {
          reflectCol.set(skyCol);
      }
      float it=ht.t;
      ht= ttree.luminares[0].Intersect(ht, origin, dv2, it );
      if (ht.t<it) 
        reflectCol.set(lgtCol);
      reflectCol.mult(Fresnel(nrm,ray.direction));
      return reflectCol;
    }
    
    PVector ColourTable[]= {
        new PVector(0.15,0.15,0.15),  // sky co
      new PVector(1.,1.,1.),   // lightCol
     new PVector(1.,.8,.6),    // obj co
      new PVector(1.,0.0,0.0),  // ground col
    
     
      new PVector(0.05*2.,0.1*2.,0.2*2.),  // sky col
      new PVector(1.2,1.1,0.9),   // lightCol
      new PVector(1.,.8,.6),    // obj col
      new PVector(54./255.,117./255.,136./255.),  // ground col
    
      new PVector(0.,0.,0.),  // sky col
      new PVector(1.,1.,1.),   // lightCol
      new PVector(0.7,0.7,0.7),    // obj col
      new PVector(0.7,0.7,0.7),  // ground col
    
    new PVector(0.1,0.1,0.1),  // sky col
      new PVector(1.,1.,1.),   // lightCol
      new PVector(0.7,0.7,0.7),    // obj col
      new PVector(1.0,0.1,0.05),  // ground col
      
      new PVector(0.1,0.1,0.1),  // sky col
      new PVector(1.,1.,1.),   // lightCol
      new PVector(0.2,1.,0.1),    // obj col
      new PVector(1.0,1.,1.),  // ground col
      
      new PVector(0.2,0.1,0.05),  // sunset
      new PVector(1.,0.6,0.3),   // lightCol
      new PVector(1.,.1,.1),    // obj col
      new PVector(0.5,1.,0.5),  // ground col
    
      new PVector(0.15,0.15,0.15),  // sky col
      new PVector(1.2,1.1,0.9),   // lightCol
      new PVector(1.,0.5,1.),    // obj col
      new PVector(1.,1.,0.4f),  // ground col
    
      new PVector(0.2,0.2,0.2),  // sky col
      new PVector(0.2,0.2,1.),   // lightCol
      new PVector(1.,1.,1.),    // obj col
      new PVector(0.,0.,0.),  // ground col
    
      new PVector(0.2,0.2,0.2),  // sky col
      new PVector(1.,1.,1.),   // lightCol
      new PVector(1.,0.2,0.5),    // obj col
      new PVector(0.5,1.0,0.2),  // ground col
    };
    
    
    PVector skyCol = new PVector(0.05*0.5,0.1*0.5,0.2*0.5);
    PVector[] frameBuffer = null;
    
    int ZOrder(int i)
    {
      return  (i&0x1) | 
        ((i&0x4)>>1)|
        ((i&0x10)>>2)|
        ((i&0x40)>>3)|
        ((i&0x100)>>4)|
        ((i&0x400)>>5)|
        ((i&0x1000)>>6)|
        ((i&0x4000)>>7)|
        ((i&0x10000)>>8)|
        ((i&0x40000)>>9)|
        ((i&0x100000)>>10);
    }
    PImage  rendered;
    PImage  renderedFrame0;
    PImage  renderedFrame1;
    final int NumRenderTiles = 16;
    
    synchronized int doColor(float x, float y, float z) {
      return color(x, y, z);
    }
    
    class RenderTileTraced implements Runnable
    {
      int tileNo;
      PVector lgtCol;
      PVector objCol;
      PVector grdCol;
      PVector  skyCol;
      Camera    cam;
      float[]   frameOff;
      int        m_frame;
      int        m_stationaryFrameNo;
      PImage     m_output;
      
      public void run() {
        
         QMCSamples sampler =new QMCSamples(0x87FA126A,0X1EACABF);
        RenderState rs = new RenderState(ttree);
        hitResult hr = new hitResult();  
        PVector lgtVCol = new PVector();
        lgtVCol.set(lgtCol);
        lgtVCol.mult(2.);
    
        float maxT = 1000.;
        Ray ray=new Ray();
        PVector htpt=new PVector();
    
        int tileWidth = width/4;
        int tileHieght=height/4;
        int sx = (tileNo&3)*tileWidth;
        int sy = (tileNo/4)*tileHieght;
        int ex = sx+tileWidth;
        int ey = sy+tileWidth;
        PVector col= new PVector(0.,0.,0.);
        PVector nrm=new PVector();
        PVector refcol=new PVector(0.,0.,0.);
        final float bias = 0.0001f;
        
        for(int y=sy;y<ey;y++)
          for(int x =sx;x<ex;x++)
          {
            int idx=(height-y-1)*width + x;
            
          if (  m_stationaryFrameNo ==1)
              frameBuffer[idx].set(0.,0.,0.);
            
            cam.generateRay( ray, x+frameOff[0],y+frameOff[1]);
       
            float [] origin = {
              ray.origin.x,ray.origin.y,ray.origin.z
            };
            float[] dv2 = {
              ray.direction.x,ray.direction.y,ray.direction.z
            };      
    
            hitResult ht = ttree.Trace(rs, origin, dv2, maxT );  
            // todo add some block based rreduction?
            hr.t=ht.t;
            hr= ttree.luminares[0].Intersect(hr, origin, dv2, ht.t );
            if (hr.t<ht.t) {
              frameBuffer[idx].add(lgtVCol);
              continue;
            }
            if ( ht.t>=maxT )
            {
              frameBuffer[idx].add(skyCol);
              continue;
            }
            htpt.set(PVector.add(PVector.mult( ray.direction,ht.t),ray.origin));
            // todo make use lod stuff
            float ds =1./0.02;//5;//1./ 2.;//(0.1*ht.t);
            int lod = constrain( (int)(ht.t*2.), 1, 4);
            ds *= (float)(5-lod);
            float off = ((float)lod-1.)*10. +ht.tri;
    
            int ix = (int)(htpt.x*ds);
            int iy= (int)((htpt.y+off)*ds);
            int iz= (int)(htpt.z*ds);
            int key1 = hash_rand(ix,iy,iz);
    
            int cacheIdx = g_RenderCache.FindIdx(key1);
            int sc= g_RenderCache.GetSampleCnt(cacheIdx,key1);
    
            int key2 = hash_rand(iz, ix, iy);
            nrm.set(ttree.triList[ht.tri].GetNormal(ray.direction, ht.u, ht.v));
            htpt.add( PVector.mult(nrm,bias));
    
            // could skip on higher rendercache sample values.
            if ( g_useReflect)
            {
              ray.origin = htpt;
              origin[0]= ray.origin.x;
              origin[1]=ray.origin.y;
              origin[2]=ray.origin.z;
              refcol=CalcReflectionColor( ht, rs, ray,origin, dv2, htpt, nrm, lgtVCol );
            }
            if ( sc > 64 ) //g_NumSamplesRequired)
            {
              PVector vc=new PVector();
              g_RenderCache.Get( cacheIdx, key1, vc, m_frame);
              vc.add(refcol);
              frameBuffer[idx].add(vc);
              continue;
            }
    
            sampler.Set(key1, key2);
            PVector rcol =ttree.Luminance( sampler, sc, rs, htpt, nrm,lgtCol,0);        
    
            PVector diff = new PVector();
            diff.set(GetDiffuse( htpt, objCol, grdCol ));
            col.set(0.,0.,0.);
            if ( g_DrawDirect) 
              col.add(rcol);
          
            col.mult(diff);
            int numBounces =2;
            for (int i =0; i < numBounces; i++)
            {
              ray.origin = htpt;
              origin[0]= ray.origin.x;
              origin[1]=ray.origin.y;
              origin[2]=ray.origin.z;
              ray.direction= GenerateDirectionVector(ray.direction, sc, sampler, nrm, 1+i);
    
              dv2[0]= ray.direction.x;
              dv2[1]=ray.direction.y;
              dv2[2]=ray.direction.z;
              ht = ttree.Trace(rs, origin, dv2, maxT );
              if ( ht.t < maxT )
              {         
                htpt.set(PVector.add(PVector.mult( ray.direction,ht.t),ray.origin));
                nrm.set(ttree.triList[ht.tri].GetNormal(ray.direction,ht.u,ht.v));
                htpt.add( PVector.mult(nrm,bias));
                PVector bdiff = GetDiffuse( htpt, objCol, grdCol);
                diff.mult(bdiff);
                rcol = ttree.Luminance( sampler, sc, rs, htpt, nrm, lgtCol, 3+i);
    
                rcol.mult(diff);
                col.add(rcol);
              } 
              else
              {
                col.add(PVector.mult(skyCol,diff));
                i=numBounces+1;
              }
            }  
        //  col.set(PVector.add(PVector.mult(nrm,.5),new PVector(.5,.5,.5)));
               
            col= g_RenderCache.Add( cacheIdx, key1,col, m_frame );
            if ( g_useReflect)
              col.add(refcol);
            frameBuffer[idx].add(col);
          }
          
          
        float wieght=1./(float)m_stationaryFrameNo;
         for(int y=sy;y<ey;y++)
          for(int x =sx;x<ex;x++)
          {
            int idx=(height-y-1)*width + x;
            m_output.pixels[idx]=doColor( sqrt(frameBuffer[idx].x*wieght)*255.,
            sqrt(frameBuffer[idx].y*wieght)*255.,
            sqrt(frameBuffer[idx].z*wieght)*255.);
          }
       
      }
    };
    
    class RenderTile implements Runnable
    {
      int tileNo;
      PVector lgtCol;
      PVector objCol;
      PVector grdCol;
      PVector  skyCol;
      Camera    cam;
      PImage     m_output;
      
      public void run() {
        RenderState rs = new RenderState(ttree);
        hitResult hr = new hitResult();  
        Ray ray=new Ray();
        PVector htpt=new PVector();
        PVector col=new PVector();
    
        int tileWidth = width/4;
        int tileHieght=height/4;
        int sx = (tileNo&3)*tileWidth;
        int sy = (tileNo/4)*tileHieght;
        int ex = sx+tileWidth;
        int ey = sy+tileWidth;
    
        for(int y=sy;y<ey;y++)
          for(int x =sx;x<ex;x++)
          {
            int idx=(height-y-1)*width + x;
            cam.generateRay( ray, x,y);
    
            float [] origin = {
              ray.origin.x,ray.origin.y,ray.origin.z
            };
            float[] dv2 = {
              ray.direction.x,ray.direction.y,ray.direction.z
            };      
            hitResult ht = ttree.Trace(rs, origin, dv2, 10. ); 
            hr.t=ht.t;
            hr= ttree.luminares[0].Intersect(hr, origin, dv2, ht.t );
            if (hr.t<ht.t) 
            {
              m_output.pixels[idx]=doColor(255.,255.,255.);
              continue;
            }
            if ( ht.t >=10. )
            { 
              m_output.pixels[idx]=doColor(0,0,0);
              continue;
            } 
            htpt.set(PVector.add(PVector.mult( ray.direction,ht.t),ray.origin));
            // todo make use lod stuff
            float ds =1./0.02;;
            int lod = constrain( (int)(ht.t*2.), 1, 4);
            ds *= (float)(5-lod);
            float off = ((float)lod-1.)*10. + ht.tri;
    
            int ix = (int)(htpt.x*ds);
            int iy= (int)((htpt.y+off)*ds);
            int iz= (int)(htpt.z*ds);
            int key1 = hash_rand(ix,iy,iz);
    
            int cacheIdx = g_RenderCache.FindIdx(key1);
            g_RenderCache.Get( cacheIdx, key1, col, frame);
            m_output.pixels[idx]=doColor( sqrt(col.x)*255,sqrt(col.y)*255,sqrt(col.z)*255);
          }
      }
    }
    
    void draw()
    {
      Camera cam = new Camera( g_camControl.m_camPos,
      new PVector(0.0f,1.0f,0.0f),
      g_camControl.m_camDir,width, height,width);
      QMCSamples sampler =new QMCSamples(0x87FA126A,0X1EACABF);
      g_colorSelection = g_colorSelection%(ColourTable.length/4);
    
      if (g_Workers.IsFinished())
      {
        if ( g_applyFlush){
          g_RenderCache.Flush();
          g_applyFlush=false;
        }
        if ( g_Animate){
        float ang =  (float)(millis()-g_FrameStart)*PI*0.0001f;
        g_camControl.update(ang,0.);
       resetAccumulationBuffer();
      }
        frame++;
        g_stationaryFrame++;
         float[] frameOff=sampler.Get(g_stationaryFrame,0);
        for (int i = 0; i < NumRenderTiles;i++)
        {
          if ( g_useCacheOnly)
          {
            RenderTile rt= new RenderTile();
            rt.skyCol = ColourTable[g_colorSelection*4];
            rt.lgtCol = ColourTable[g_colorSelection*4+1];
            rt.objCol = ColourTable[g_colorSelection*4+2];
            rt.grdCol = ColourTable[g_colorSelection*4+3];
            rt.tileNo=i;
            rt.cam =cam;
            rt.m_output=rendered;
            g_Workers.execute(rt);
          }
          else
          {
            RenderTileTraced rt= new RenderTileTraced();
            rt.skyCol = ColourTable[g_colorSelection*4];
            rt.lgtCol = ColourTable[g_colorSelection*4+1];
            rt.objCol = ColourTable[g_colorSelection*4+2];
            rt.grdCol = ColourTable[g_colorSelection*4+3];
            rt.m_stationaryFrameNo = g_stationaryFrame;
            rt.m_frame = frame;
            rt.tileNo=i;
            rt.cam =cam;
            rt.frameOff = frameOff;
            rt.m_output=rendered;
            g_Workers.execute(rt);
          }
        }
        background(rendered);
        fill(~0);
        text("Frame Time "+ (millis()-g_FrameStart), 20,20);
        text( g_Scenes[g_CurrentScene%g_Scenes.length].des,10,height-10);
        g_FrameStart=millis();
        // swap the frame buffer
    
        
      }
    }
    
    
    PVector minv( PVector a, PVector b)
    {
      return new PVector(min(a.x,b.x),min(a.y,b.y),min(a.z,b.z));
    }
    PVector maxv( PVector a, PVector b)
    {
      return new PVector(max(a.x,b.x),max(a.y,b.y),max(a.z,b.z));
    }
    
    Node BuildNode(PrimBound[] pbounds, AABB bound, Tri[] tris, 
    TriAccel[] fTris, int lstart, int lend, int nodeIdx) {
      PVector bmin = tris[pbounds[lstart].idx].v0;
      PVector bmax = tris[pbounds[lstart].idx].v0;
      AABB bc2 = new AABB();
      for (int i =lstart; i <lend;i++)
      {
        int ix = pbounds[i].idx;
        bmin = minv(tris[ix].v0, bmin);
        bmin = minv(tris[ix].v1, bmin);
        bmin = minv(tris[ix].v2, bmin);   
    
        bmax = maxv(tris[ix].v0, bmax);
        bmax = maxv(tris[ix].v1, bmax);
        bmax = maxv(tris[ix].v2, bmax); 
       // bc2.grow(tris[ix].v0);
        //bc2.grow(tris[ix].v1);
        //bc2.grow(tris[ix].v2);
        bc2.grow( pbounds[i].bound);
        
     //   println( "Tri " + tris[ix].v0 + " "+tris[ix].v1+" "+tris[ix].v2);
        fTris[i]=new TriAccel(tris[ix],ix);
      }
      
     // println( bound.bmin + " " + bmin + " " + bc2.bmin );
     // assert( bound.bmin.x == bmin.x);
     // assert( bound.bmin.y == bmin.y);
     //   assert( bound.bmin.z == bmin.z);
        
      //println( bound.bmax  + " " + bmax + " " + bc2.bmax );   
       //   assert( bound.bmax.x == bmax.x);
      //assert( bound.bmax.y == bmax.y);
       // assert( bound.bmax.z == bmax.z);
    
      Node n = new Node();
      n.minx = bmin.x;
      n.miny = bmin.y;
      n.minz = bmin.z;
    
      n.maxx = bmax.x;
      n.maxy = bmax.y;
      n.maxz = bmax.z;
    
      n.leafStart =lstart;
      n.leafEnd = lend;
      n.childIdx = -nodeIdx-1;
      //println("ChildNode: nodeIdx "+nodeIdx + " ChildIdx " +  n.childIdx );
      return n;
    }
    class nodeInfo
    {
      Node n;
      int nidx;
      nodeInfo( Node _n, int _nidx )
      {
        n=_n; 
        nidx=_nidx;
      }
    };
    
    void DumpTreeInfo(Node[] nodes, int idx,int level) {
      Node n = nodes[idx];
      for (int i =0; i < level;i++)
        print("  ");
    
      print( n.childIdx >=0 ? "Inner " : "Leaf ");
      print( " | Idx "+idx );
      print(" | Bx "+ n.minx +" "+n.maxx +" | By "+ n.miny +" "+n.maxy );
      if ( n.childIdx <0 )
      {
        println(" | Leaf Start " + n.leafStart +" | Leaf End " + n.leafEnd);
      }
      else
        println(" |  Child "+n.childIdx);
    
      if ( n.childIdx>=0)
      {
        DumpTreeInfo( nodes, n.childIdx, level+1);
        DumpTreeInfo( nodes, n.childIdx+1, level+1);
      }
    }
    class PrimBound
    {
      AABB bound;
      float[] centre;
      int idx;
    
      PrimBound( Tri t, int i )
      {
        bound = new AABB();
        bound.grow(t.v0);
        bound.grow(t.v1);
        bound.grow(t.v2);
    
        idx = i;
        PVector bc =  PVector.mult(PVector.add(bound.bmin, bound.bmax),0.5);
        centre=new float[3];
        centre[0]=bc.x;
        centre[1]=bc.y;
        centre[2]=bc.z;
      }
    };
    
    class BuildState
    {
      AABB[]       m_rightBounds;
      PrimBound[]  m_primBounds;
      AABB         m_totalBound;
      BuildState( Tri[] trilist)
      {
        m_primBounds = new PrimBound[trilist.length];
        m_rightBounds= new AABB[trilist.length];
        m_totalBound = new AABB();
        for (int i = 0; i < trilist.length;i++) {
          m_primBounds[i]=new PrimBound( trilist[i],i);
          m_rightBounds[i]= new AABB();
          m_totalBound.grow(m_primBounds[i].bound);
        }
        println("TB "+m_totalBound.bmin+ " " + m_totalBound.bmax);
      }
      void SortBounds( int dim, int lstart, int lend)
      {
        float shrinkFactor = 1.247330950103979f;
    
        boolean swaps=true;
        int gap = lend-lstart;
    
        while( gap > 1 || swaps == true)
        {
          assert( gap >0);
          if ( gap>1)
            gap = (int)((float)gap/shrinkFactor);
    
          swaps = false;
          for (int i =lstart; i <lend-gap;i++)
          {
            if ( m_primBounds[i].centre[dim] > m_primBounds[i+gap].centre[dim])
            {
              PrimBound p = m_primBounds[i];
              m_primBounds[i] = m_primBounds[i+gap];
              m_primBounds[i+gap] = p;
              swaps = true;
            }
          }
        }
      }
    
      final float TriangleCost = 3;
      final float NodeCost = 3;
    
      float getTriangleCost(int n) { 
        return TriangleCost*n;
      }
      float NodeSAH() { 
        return NodeCost;
      }
    };
    
    class AABB
    {
      PVector bmin;
      PVector bmax;
      AABB() {
        bmin = new PVector(10000000.f,10000000.f,10000000.f);
        bmax = new PVector(-10000000.f,-10000000.f,-10000000.f);
      }
      void grow( final PVector other) {
        bmin = minv( bmin, other); 
        bmax = maxv( bmax, other);
      }
      void grow( final AABB other) {
        bmin = minv( bmin, other.bmin); 
        bmax = maxv( bmax, other.bmax);
      }
      void set( final AABB other)
      {
        bmin.set(other.bmin);
        bmax.set(other.bmax);
      }       
      float area() {
        PVector d = PVector.sub(bmax,bmin);
        return (d.x*d.y + d.y*d.z + d.z*d.x)*2.f;
      }
    };
    
    
    nodeInfo BuildTree(BuildState bs, AABB bound, Node[] nodes, int cidx, int nodeIdx,
    Tri[] tris, TriAccel[] fTris, int lstart, int lend,int level, int maxLevel ) {
      // calc bounds
      // check if leaf
      int num = lend-lstart;
      // compute SAH split
      float splitSAH=bs.getTriangleCost(num)*bound.area();
      int splitAxis =-1;
      int  splitNumLeft=0;
      AABB splitLeftBounds = new AABB();
      AABB splitRightBounds = new AABB();
      float nodeSAH = bs.NodeSAH();
    
      for (int dim =0; dim < 3; dim++)
      {
        //sort on centres
        bs.SortBounds( dim,lstart, lend);
        // get right bounds.
        AABB rightBounds = new AABB();
        for (int i = lend-1; i >lstart; i--)
        {
          rightBounds.grow(bs.m_primBounds[i].bound);
          bs.m_rightBounds[i - 1].set( rightBounds);
          assert( bs.m_primBounds[i-1].centre[dim]<=bs.m_primBounds[i-1].centre[dim]);
        }
        // Sweep left to right and select lowest SAH.
        AABB leftBounds = new AABB();
        for (int i = lstart+1; i < lend; i++)
        {
          leftBounds.grow(bs.m_primBounds[i-1].bound);
          int numLeft = i-lstart;
          int numRight = lend-i;
          assert( (numLeft + numRight) == num);
          float sah = nodeSAH + leftBounds.area() * bs.getTriangleCost(numLeft) 
            + bs.m_rightBounds[i - 1].area() * bs.getTriangleCost(numRight);
          if (sah < splitSAH)
          {
            splitSAH = sah;
            splitAxis = dim;
            splitNumLeft = i;
            splitLeftBounds.set(leftBounds);
            splitRightBounds.set( bs.m_rightBounds[i - 1]);
          }
        }
      }
      if ( splitAxis==-1)
      {
         // println("No Axis found lstart "+lstart + " lend " +lend);
        return new nodeInfo(BuildNode( bs.m_primBounds, bound,tris, fTris, lstart, lend,cidx), nodeIdx);
      }   
     
        //sort on centres
      bs.SortBounds( splitAxis,lstart, lend); 
        
      int end = splitNumLeft;
    
      assert( end != lend);
         assert( end != lstart);
    
      Node n = new Node();
      n.minx = bound.bmin.x;
      n.miny = bound.bmin.y;
      n.minz = bound.bmin.z;
    
      n.maxx = bound.bmax.x;
      n.maxy = bound.bmax.y;
      n.maxz = bound.bmax.z;
      n.childIdx = nodeIdx;
      // println("\ninternal node : nodeIdx " +nodeIdx + " Level "+level);
    
      nodeInfo lni = BuildTree(bs,splitLeftBounds, nodes, nodeIdx, nodeIdx+2, tris, fTris,lstart,end, level+1,maxLevel);
      nodes[nodeIdx]= lni.n;
      nodeInfo rni=BuildTree(bs,splitRightBounds,nodes, nodeIdx+1, lni.nidx,tris, fTris,end,lend,level+1,maxLevel);
      nodes[nodeIdx+1]= rni.n;
      return new nodeInfo(n,  rni.nidx);
    }
    
    
    // Thread safe spatial cache
    import java.util.concurrent.atomic.AtomicIntegerArray;
    
    class RenderCache
    {
      AtomicIntegerArray keys;
      float[] col;
    
      final int csize = 4096*512;
      final int csizemask = csize-1;
      final int cmask = (csize/4)-1;
    
      RenderCache() {
        keys = new AtomicIntegerArray(csize);
        col = new float[csize*4];
    
        Flush();  
      }
      void Flush()
      {
        for (int i = 0; i < csize; i++)
          keys.set(i,0xFFFFFFFF);
    
      }
      int FindIdx( int k )
      {
        int idx = k&cmask;
        int eIdx =idx+8;
        for (;idx < eIdx;idx++)
        {
          int kv = keys.get(idx&csizemask) ;
          if ( kv ==k || kv ==0xFFFFFFFF) 
             return idx&csizemask;
        }
        return idx&csizemask;
      }
      
      void Get( int idx, int k, PVector v, float fno)
      {    
        if ( keys.get(idx) !=  k ) 
        {
         v.set(0.,0.,0.);
         return;
        }
        int idx4=idx*4;
        v.set(col[idx4+0], col[idx4+1],col[idx4+2]);
        v.mult(1./(float)col[idx4+3]);
        return;
      }
    
      PVector Add( int idx, int k, PVector v, int fno )
      {
        int kv = keys.get(idx);
        
        // only set the key once
        if ( kv!=  k )  // mm, how to do this across mm threads
        {
          if ( kv !=0xFFFFFFFF )
          {
            return v;
          }
          if (!keys.compareAndSet( idx,0xFFFFFFFF,k) )
          {
            return v;
          }
          //assert(kv==k);
          col[idx*4+0]=0.f;
          col[idx*4+1]=0.f;
          col[idx*4+2]=0.f;
          col[idx*4+3]=0.f;
        }
         int idx4 = idx*4;
        col[idx4+0]+=v.x;
        col[idx4+1]+=v.y;
        col[idx4+2]+=v.z;
        col[idx4+3]=col[idx4+3]+1.f;
    
        v.set(col[idx*4+0], col[idx*4+1],col[idx*4+2]);
        v.mult(1.f/(float)col[idx4+3]);
        return v;
      }
      int GetSampleCnt( int idx, int k )
      {
        return keys.get(idx)==k ? ((int)col[idx*4+3]) : 0;
      }
    };
    
    
    // objloader
    class Tri
    {
      PVector v0;
      PVector v1;
      PVector v2;
      
      PVector n0;
      
      PVector sn0;
      PVector sn1;
      PVector sn2;
      
      int     materialId;
      
       Tri( PVector A, PVector B, PVector C)
       {
         v0 = A; v1 = B; v2 = C;
         n0 =PVector.sub(C,A);
         n0.normalize();
         PVector n2 = PVector.sub(B,A);
         n2.normalize();
         n0=n2.cross( n0);
         n0.normalize();
         //n0.set(0.,1.,0.);
         sn0=new PVector();
         sn0.set(n0);
         sn1=new PVector();
         sn1.set(n0);
         sn2=new PVector();
         sn2.set(n0);
       }
       boolean IsDegenerate()
       {
         return n0.x==0. && n0.y==0. && n0.z==0.;
       }
       
       Tri( PVector A, PVector B, PVector C,
               PVector _n0, PVector _n1, PVector _n2)
       {
         v0 = A; v1 = B; v2 = C;
         n0 =PVector.sub(C,A);
         n0.normalize();
         PVector n2 = PVector.sub(B,A);
         n2.normalize();
         n0=n2.cross( n0);
         n0.normalize();
         
         sn0=new PVector();
         sn0.set(_n0);
         sn1=new PVector();
         sn1.set(_n1);
         sn2=new PVector();
         sn2.set(_n2);
       }
       boolean onLeft( PVector v)
       {
         return v0.x <= v.x && v0.y <= v.y  && v0.z <= v.z &&
                   v1.x <= v.x && v1.y <= v.y  && v1.z <= v.z &&
                   v2.x <= v.x && v2.y <= v.y  && v2.z <= v.z;
       }
       PVector GetNormal( PVector I, float u, float v )
       {
         float w = 1.f-u-v;
         PVector n = PVector.add(PVector.mult(sn1,u),PVector.mult(sn0,w));
         n = PVector.add(PVector.mult(sn2,v), n);
         n.normalize();
         //n.mult( I.dot(n)>0 ? -1 :1);
         return n;
       }
    };
    
    class Luminare
    {
      final PVector v0;
      final PVector e0;
      final PVector e1;
      final PVector n0;
      final float   area;
      
      final PVector v1;
      final PVector v2;
      final PVector v3;
      
      final TriAccel m_triA;
      final TriAccel m_triB;
       
      Luminare( PVector A, PVector B, PVector C,PVector D)
      {
        
        println("Luminarie " + A+B+C+D);
        v1=B;
        v2 =C;
        v3 =D;
        v0 = A; e0 = PVector.sub(B,A); e1 = PVector.sub(D,A);
        n0=e0.cross( e1);
        area = e0.mag()*e1.mag();
        n0.normalize();
        m_triA = new TriAccel(A,B,C,0);
        m_triB = new TriAccel(C,D,A,1);
      }
      PVector Sample( float x, float y)
      {
       return PVector.add(PVector.add( PVector.mult(e0,x), PVector.mult(e1,y)),v0);
      }
      hitResult Intersect(hitResult hit, float[] org, float[] dir, float maxT )
      {
        float it = hit.t;
        hit = m_triA.Intersect(hit,org,dir,it);
        if (hit.t >=it )
          hit = m_triB.Intersect(hit,org,dir,hit.t);
         return hit;
      }
      float SolidAngle( PVector pt, PVector nrm)
      {
        PVector d = PVector.sub(pt,v0);
        float r2 = max(PVector.dot(d,d),0.0001);
        return area/r2;
      }
    };
    
    /*
    this class stores the mtl information.
    I am assuming that only ambient, diffuse, specular and shininess values are specified.
    ahmet.kizilay@gmail.com
    */
    class MTLInfo{
       PVector ka;
       PVector kd;
       PVector ks;
       float ns;
       String name;
       MTLInfo(){  
       }
       void setAmbient(PVector am){
          ka =new PVector();
          ka.set(am);
       }
       void setDiffuse(PVector dif){
          kd =new PVector();
          kd.set(dif);
       }
       void setSpecular(PVector spec){
          ks =new PVector();
          ks.set(spec);
       }
       void setShininess(float shine){
          ns = shine;
       }      
       void setName(String _name){
          name = _name;
       }
    }
    
    
    class OBJReader{
      String fileName;
      Vector vertices;
      Vector normals;
      Vector tris;
      Vector textures;
      Vector mtlinfo;
      
      Vector luminares;
      boolean isLuminare=false;
      
      boolean mtlincluded = false;
      int mtlnum = -1; // mtlnum is initially -1 and it keeps the current index of the mtl.
     
     
      OBJReader(String fn, boolean mtlexists){
        fileName = fn;
        vertices = new Vector();
        normals = new Vector();
        textures = new Vector();
        mtlinfo = new Vector();
        tris = new Vector();
        luminares = new Vector();
        mtlincluded = mtlexists;
      }
     Tri[]  GenerateTriList()
     {
       Tri[] trilist = new Tri[tris.size()];
       for ( int i=0;i<trilist.length;i++)
       {
         trilist[i]=(Tri)tris.get(i);
       }
       return trilist;
     }
     
     Luminare[] GetLuminares()
     {
        Luminare[] li = new Luminare[luminares.size()];
       for ( int i=0;i<li.length;i++)
       {
         li[i]=(Luminare)luminares.get(i);
       }
       return li;
     }
     
      void readFile(){
        println("Loaded File "+fileName);
        if(mtlincluded) readMTL();
     
        try
        {
          BufferedReader input;
          input = createReader(fileName);
          try
          {
            String newLine = null;
            while((newLine = input.readLine()) != null){
             // println(newLine);
              int ind = newLine.indexOf("vn ");
              if(ind != -1){
                readNormal(newLine);
                continue;
              }
     
              ind = newLine.indexOf("v ");
              if(ind != -1){
                readVertex(newLine);
                continue;
              }
     
              ind = newLine.indexOf("f ");
              if(ind != -1) {
                readPolygon(newLine);
                continue;
              }    
     
              ind = newLine.indexOf("vt ");
              if(ind != -1){
               // readTexture(newLine);
                continue;
              }
     
              ind = newLine.indexOf("usemtl ");
              if(ind != -1){
                readMTLInfo(newLine);
                continue;
              }
            }
          }
          finally{
            input.close();
          }
       }
        catch(IOException ex){
          ex.printStackTrace();
        }
      }
     
      void readVertex(String newLine){
        String pieces[] = splitTokens(newLine, " ");
        PVector vert = new PVector(float(pieces[1]), float(pieces[2]), float(pieces[3]));
        vertices.addElement(vert);
      }
     
      void readNormal(String newLine){
        String pieces[] = splitTokens(newLine, " ");
        PVector norms = new PVector(float(pieces[1]), float(pieces[2]), float(pieces[3]));
        normals.addElement(norms);
      }
     
      void readTexture(String newLine){
        String pieces[] = splitTokens(newLine, " ");
        PVector tex = new PVector(float(pieces[1]), float(pieces[2]));
        textures.addElement(tex);
      }
      
      int GetV( int v )
      {
        if (v<=0)
        {
          println(v);
          v+=1;
        }
          
        assert( v >0);
        return v-1;
      }
     
      void readPolygon(String newLine){
        String pieces[] = splitTokens(newLine, " ");
        int npts = pieces.length-1;
        assert( pieces.length>=4 );
          
        PVector[] pts = new PVector[npts];
        boolean hasNrms = false;
        PVector[] nrms = new PVector[npts];
       // PVector[] tex = new PVector[npts];
        
        for(int i = 1; i < pieces.length; i++){
          String smallerPieces[] = splitTokens(pieces[i], "//");
          int ip=i-1;
          pts[ip]= (PVector)vertices.get( GetV(int(smallerPieces[0])));
          
          //if texture is not specified
         if(smallerPieces.length == 2){
            hasNrms =true;
            nrms[ip]=(PVector)normals.get(GetV(int(smallerPieces[1])));
          }
          //if texture is specified
          if(smallerPieces.length == 3){
            hasNrms =true;
           nrms[ip]=(PVector)normals.get(GetV(int(smallerPieces[2])));
         //   tex[ip]=(PVector)normals.get( GetV(int(smallerPieces[1])));
          }
        }
        
        if ( isLuminare)
        {
          println("Found Luminare");
          luminares.addElement( new Luminare(pts[0],pts[1],pts[2],pts[3]));
        }
        else
        {
          // now build the triangles
          for (int i=2; i < npts; i++)
          {
            Tri t;
             if ( hasNrms)
                t = new Tri( pts[0],pts[i-1],pts[i],nrms[0],nrms[i-1],nrms[i]);
             else
               t = new Tri( pts[0],pts[i-1],pts[i]);
               
            if ( !t.IsDegenerate())
               tris.addElement(t);
          }
        }
      }
     
      void readMTLInfo(String newLine){
        String pieces[] = splitTokens(newLine, " ");
       
     if ( mtlincluded == false)
     {
       isLuminare=pieces[1].equals("diffuseLuminaire1SG");
     }
        int mtlsize = mtlinfo.size();
        int i = 0;
        boolean mtlfound = false;
        while(i < mtlsize && !mtlfound){
          MTLInfo mtl = (MTLInfo) mtlinfo.elementAt(i);
          if(pieces[1].equals(mtl.name)) {
            mtlfound = true;
            mtlnum = i;
          }
          i++;
        }
      }
    
      void readMTL(){
     
        String name = fileName;
        name = name.substring(0, name.length() - 3);
        name = name + "mtl";
     
        try{
          BufferedReader input;
          input = createReader(name);
          try{
            String newLine = null;
            while((newLine = input.readLine()) != null){
              int ind = newLine.indexOf("newmtl ");
              if(ind != -1){
                MTLInfo mtl = new MTLInfo();
                String pieces[] = splitTokens(newLine, " ");
                mtl.setName(pieces[1]);
                for(int i = 0; i < 4; i++){
                  String material;
                  if((material = input.readLine()) != null){
                    String smallerPieces[] = splitTokens(material, " ");
                    if(i == 0) mtl.setAmbient(new PVector(float(smallerPieces[1]), float(smallerPieces[2]), float(smallerPieces[3])));
                    if(i == 1) mtl.setDiffuse(new PVector(float(smallerPieces[1]), float(smallerPieces[2]), float(smallerPieces[3])));
                    if(i == 2) mtl.setSpecular(new PVector(float(smallerPieces[1]), float(smallerPieces[2]), float(smallerPieces[3])));
                    if(i == 3) mtl.setShininess(float(smallerPieces[1]));
                  }
                }
                mtlinfo.addElement(mtl);
              }
            }
          }
          finally{
            input.close();
          }
        }
        catch(IOException ex){
          ex.printStackTrace();
        }
      }
    }
    
    
    //sampling
    
    // from  http://jgt.akpeters.com/papers/ShirleyChiu97/
    
    /* seedx, seedy is point on [0,1]^2.  x, y is point on radius 1 disk */
    
    PVector to_unit_disk( float seedx, float seedy, float angle) {
      float phi, r;
      float a = 2.*seedx - 1.;   /* (a,b) is now on [-1,1]^2 */
      float b = 2.*seedy - 1.;
      if (a > -b) {     /* region 1 or 2 */
        if (a > b) {  /* region 1, also |a| > |b| */
          r = a;
          phi = (PI*1./4 ) * (b/a);
        }
        else  {  /* region 2, also |b| > |a| */
        r = b;
        phi = (PI*1./4 ) * (2 - (a/b));
        }
      }
      else {        /* region 3 or 4 */
        if (a < b) {  /* region 3, also |a| >= |b|, a != 0 */
            r = -a;
            phi = (PI*1./4 ) * (4 + (b/a));
        }
        else {  /* region 4, |b| >= |a|, but a==0 and b==0 could occur. */
            r = -b;
            if (b != 0)
                phi = (PI*1./4 ) * (6 - (a/b));
              else
                  phi = 0;
        }
      }
      float x = (r * cos(phi))*angle;
      float y = (r * sin(phi))*angle;
     float z =sqrt(1.-x*x-y*y);
      return new PVector( x,z,y);  //
    }
    
     float Sobol2(int n2,  int r2) {
      long n = n2;
      long r = r2;
      for (long  v = 1 << 31L; n != 0L; n >>= 1L, v ^= v >> 1L)
             if ((n & 0x1L) != 0 )
                  r ^= v;
      r  = (r &0xFFFFFFFFL);
      return ((float)r)/  ( (float)0x100000000L); } 
      
    
    
    
    long SobolMatrices[]={
    2147483648L, 3221225472L, 1610612736L, 805306368L, 402653184L, 201326592L, 100663296L, 50331648L, 25165824L, 12582912L, 6291456L, 3145728L, 1572864L, 786432L, 393216L, 196608L, 98304L, 49152L, 24576L, 12288L, 6144L, 3072L, 1536L, 768L, 384L, 192L, 96L, 48L, 24L, 12L, 6L, 2147483648L, 1073741824L, 1610612736L, 1342177280L, 2013265920L, 1140850688L, 1711276032L, 1426063360L, 2139095040L, 1077936128L, 1616904192L, 1347420160L, 2021130240L, 1145307136L, 1717960704L, 1431633920L, 2147450880L, 1073758208L, 1610637312L, 1342197760L, 2013296640L, 1140868096L, 1711302144L, 1426085120L, 2139127680L, 1077952576L, 1616928864L, 1347440720L, 2021161080L, 1145324612L, 1717986918L, 2147483648L, 3221225472L, 2684354560L, 1342177280L, 3623878656L, 2617245696L, 1912602624L, 3372220416L, 2810183680L, 1556086784L, 3533701120L, 2572156928L, 2136473600L, 3227254784L, 2698117120L, 1352204288L, 3632234496L, 2629877760L, 1923129344L, 3377483776L, 2807617536L, 1549573120L, 3537007104L, 2576992512L, 2147428224L, 3221265600L, 2684383904L, 1342228816L, 3623921496L, 2617269340L, 1912656594L, 2147483648L, 3221225472L, 1610612736L, 4026531840L, 3087007744L, 2617245696L, 1442840576L, 855638016L, 3649044480L, 1874853888L, 3974103040L, 2909798400L, 2396520448L, 1544290304L, 906362880L, 3272540160L, 1636532224L, 4090085376L, 3135594496L, 2658349056L, 1465227264L, 868940800L, 3672690176L, 1853216512L, 4015746432L, 2948595904L, 2363490400L, 1567621360L, 920125624L, 3222012060L, 1611006038L, 2147483648L, 1073741824L, 3758096384L, 2415919104L, 3623878656L, 603979776L, 2785017856L, 1694498816L, 1166016512L, 4148166656L, 2665480192L, 3553624064L, 980942848L, 3074686976L, 2129526784L, 1138294784L, 3799351296L, 2470854656L, 3638845440L, 652169216L, 2817644544L, 1678324736L, 1175293440L, 4110576896L, 2643099520L, 3544596544L, 954491104L, 3068088464L, 2146556120L, 1074344996L, 3758984870L, 2147483648L, 1073741824L, 536870912L, 4026531840L, 1476395008L, 1677721600L, 3523215360L, 2801795072L, 964689920L, 2956984320L, 2015363072L, 2498756608L, 2321022976L, 3278110720L, 3948019712L, 401014784L, 1100578816L, 603996160L, 4060094464L, 1459679232L, 1635801088L, 3560989696L, 2854277632L, 871409408L, 3017292160L, 1931784256L, 2473752608L, 2199360752L, 3418425944L, 3882124132L, 424856402L, 2147483648L, 3221225472L, 536870912L, 268435456L, 2281701376L, 1140850688L, 1711276032L, 1996488704L, 4152360960L, 3082813440L, 3617587200L, 2815426560L, 1596456960L, 3818127360L, 961413120L, 2497904640L, 3459874816L, 596951040L, 424550400L, 2229473280L, 1178142720L, 1737786368L, 2135836160L, 4091766528L, 2981756800L, 3495179200L, 2834207136L, 1412677840L, 4002521256L, 868274260L, 2443225326L, 2147483648L, 1073741824L, 3758096384L, 1342177280L, 3355443200L, 469762048L, 301989888L, 788529152L, 1317011456L, 3695181824L, 870318080L, 1544552448L, 4063756288L, 2132017152L, 2259812352L, 3225878528L, 568623104L, 1931198464L, 3165364224L, 2741178368L, 3042899968L, 2621441024L, 3554676224L, 202376448L, 976751744L, 1662255552L, 2494693664L, 4014408432L, 1868858600L, 2942029252L, 2404098878L, 2147483648L, 1073741824L, 3758096384L, 805306368L, 2013265920L, 1409286144L, 2382364672L, 218103808L, 3397386240L, 3896508416L, 4200595456L, 2505048064L, 2832728064L, 421265408L, 2762342400L, 3579379712L, 1222475776L, 688275456L, 3702071296L, 2171949056L, 3337762816L, 605701120L, 370093568L, 1764512512L, 1018150784L, 2975596864L, 3200648928L, 1879638992L, 2550500136L, 1679444932L, 4128188230L, 2147483648L, 3221225472L, 2684354560L, 2415919104L, 939524096L, 469762048L, 3456106496L, 285212672L, 2122317824L, 4206886912L, 3177185280L, 1563426816L, 682098688L, 1742995456L, 904790016L, 2868314112L, 628850688L, 3507503104L, 3736707072L, 1792913408L, 2238896128L, 1092881408L, 3868068352L, 1992755456L, 1263764352L, 1345601984L, 2553785568L, 2350559504L, 4129382376L, 219169708L, 2964955606L, 2147483648L, 3221225472L, 3758096384L, 1342177280L, 1476395008L, 3556769792L, 436207616L, 1895825408L, 3984588800L, 3695181824L, 2455764992L, 3111124992L, 2245525504L, 1143209984L, 2720923648L, 4112580608L, 2945482752L, 2036940800L, 1709039616L, 337334272L, 4196923392L, 554519552L, 3046370816L, 139089152L, 2289008000L, 3359637824L, 1751125920L, 2557477392L, 811076440L, 1281622148L, 709757062L, 2147483648L, 1073741824L, 1610612736L, 2415919104L, 2281701376L, 3825205248L, 1442840576L, 754974720L, 3581935616L, 3955228672L, 1944059904L, 3754950656L, 3455582208L, 2264137728L, 3327787008L, 2752315392L, 909082624L, 3172352000L, 1568858112L, 264810496L, 635181056L, 4074128384L, 409134592L, 1830336768L, 3047428224L, 2078427712L, 4225295200L, 1002814416L, 2616434136L, 2881488124L, 333448286L, 2147483648L, 3221225472L, 1610612736L, 3489660928L, 671088640L, 2885681152L, 3992977408L, 3472883712L, 947912704L, 3862953984L, 2447376384L, 254803968L, 1486356480L, 913571840L, 3119120384L, 2736455680L, 3064168448L, 4182228992L, 2170822656L, 1164087296L, 661125120L, 4134800384L, 3655992832L, 1931152640L, 2661515904L, 1430768320L, 1868836576L, 2321718512L, 535296904L, 271843940L, 1208618782L, 2147483648L, 1073741824L, 536870912L, 1879048192L, 2550136832L, 738197504L, 1711276032L, 1627389952L, 1350565888L, 3921674240L, 3059744768L, 1253048320L, 110624768L, 805568512L, 3087138816L, 1543962624L, 4262035456L, 1292025856L, 914776064L, 2294681600L, 3873769472L, 2742983680L, 2968741376L, 2058398464L, 3197790592L, 1811939392L, 1174405152L, 285212784L, 3363831960L, 3317694508L, 3495952486L, 2147483648L, 3221225472L, 2684354560L, 268435456L, 134217728L, 1811939328L, 2181038080L, 1090519040L, 3816816640L, 4089446400L, 4225761280L, 2521825280L, 374865920L, 1448869888L, 3058827264L, 1180499968L, 3193077760L, 710656000L, 1012572160L, 1801441280L, 3755251712L, 2559634432L, 607856128L, 248088832L, 845300608L, 1485892800L, 2218469024L, 516524304L, 979518344L, 881913004L, 104539682L, 2147483648L, 3221225472L, 3758096384L, 4026531840L, 2818572288L, 3690987520L, 1308622848L, 1459617792L, 2172649472L, 1119879168L, 2686451712L, 1378877440L, 4208459776L, 634126336L, 1773010944L, 1023344640L, 3170009088L, 4243701760L, 1559420928L, 217477120L, 4109670400L, 3506105344L, 3136446976L, 2281607936L, 998233472L, 3309158592L, 2606555360L, 2533978352L, 1636178088L, 2999659740L, 136348238L, 2147483648L, 1073741824L, 3758096384L, 4026531840L, 3355443200L, 3690987520L, 1912602624L, 184549376L, 4236247040L, 2705326080L, 325058560L, 999292928L, 362283008L, 2895380480L, 2065039360L, 4141547520L, 1559920640L, 2953232384L, 671457280L, 738357248L, 3120850944L, 3607585792L, 2391399936L, 2856918272L, 4025460608L, 2597545024L, 117432032L, 2533498352L, 1854744392L, 1514251420L, 669810322L, 2147483648L, 3221225472L, 2684354560L, 268435456L, 4160749568L, 1006632960L, 570425344L, 1358954496L, 461373440L, 3485466624L, 706740224L, 1708130304L, 232259584L, 2312896512L, 2010251264L, 4108124160L, 3059580928L, 1444855808L, 2784993280L, 2903232512L, 2575702016L, 2395679744L, 3383671296L, 2518446336L, 100263808L, 3171626176L, 1636743840L, 2999742736L, 3953256312L, 3340184828L, 511490690L, 2147483648L, 1073741824L, 2684354560L, 2952790016L, 939524096L, 2080374784L, 2852126720L, 922746880L, 2407530496L, 843055104L, 4179623936L, 479199232L, 2048393216L, 3738435584L, 456261632L, 240058368L, 4079583232L, 2610741248L, 3448840192L, 2425401344L, 1209030656L, 635393024L, 112650752L, 1997155584L, 797643648L, 2185466944L, 3240610464L, 1620745648L, 3491449784L, 3923558460L, 2494636042L, 2147483648L, 3221225472L, 1610612736L, 1879048192L, 939524096L, 2751463424L, 2516582400L, 2600468480L, 2709520384L, 281018368L, 1222639616L, 2624585728L, 860356608L, 204210176L, 1001783296L, 2952986624L, 1476493312L, 3556884480L, 2919292928L, 1057132544L, 931289088L, 2344774656L, 3915548160L, 2360361728L, 2074682240L, 2422108608L, 150916384L, 3157209264L, 1672965848L, 1677721612L, 4127195142L, 2147483648L, 1073741824L, 3758096384L, 1342177280L, 3892314112L, 3288334336L, 973078528L, 2164260864L, 3246391296L, 541065216L, 1893728256L, 2570059776L, 1555562496L, 1724645376L, 3872260096L, 650969088L, 114130944L, 1993097216L, 4006486016L, 2999930880L, 3570321408L, 869094400L, 357420544L, 328062208L, 1705774464L, 2327818688L, 957710176L, 3966905072L, 3755533400L, 3401560516L, 3642064750L, 2147483648L, 1073741824L, 536870912L, 4026531840L, 3355443200L, 1811939328L, 3456106496L, 3942645760L, 2290089984L, 1279262720L, 1071644672L, 579862528L, 3830972416L, 2212233216L, 3577085952L, 2868969472L, 2826993664L, 3158556672L, 4158857216L, 1318170624L, 710621184L, 1759489024L, 1572086272L, 3879743744L, 2539683712L, 2664647232L, 330871136L, 3444555632L, 4285426904L, 3285975044L, 4113956866L, 2147483648L, 1073741824L, 536870912L, 3489660928L, 2013265920L, 2214592512L, 3456106496L, 3405774848L, 3296722944L, 3997171712L, 446693376L, 3161456640L, 1799880704L, 3577479168L, 1994784768L, 2923495424L, 983597056L, 1819492352L, 323608576L, 1363021824L, 3102029824L, 1698835456L, 4263744000L, 2184247552L, 166300288L, 3981238720L, 3551302944L, 2961012976L, 2294565272L, 745750532L, 860479490L, 2147483648L, 1073741824L, 3758096384L, 2415919104L, 1744830464L, 3556769792L, 570425344L, 4278190080L, 1484783616L, 3451912192L, 253755392L, 1641021440L, 1344798720L, 1208221696L, 1688076288L, 4211146752L, 4068704256L, 1992441856L, 497131520L, 2266189824L, 629704704L, 454556672L, 1662476800L, 521549568L, 3373644928L, 2752956224L, 3683999904L, 1121267728L, 2922723960L, 298057732L, 2847801358L, 2147483648L, 1073741824L, 536870912L, 1879048192L, 3087007744L, 2751463424L, 1577058304L, 620756992L, 2675965952L, 3300917248L, 245366784L, 3968860160L, 2209873920L, 1060372480L, 1968570368L, 1466236928L, 3629023232L, 4106272768L, 2527223808L, 949178368L, 1706792960L, 3191211008L, 3033859584L, 3064366848L, 1231864960L, 3700449600L, 465065056L, 3946061456L, 2471904184L, 3597051204L, 433009250L, 2147483648L, 3221225472L, 1610612736L, 805306368L, 2281701376L, 3422552064L, 3724541952L, 2734686208L, 3045064704L, 4206886912L, 2170552320L, 1099956224L, 567803904L, 295436288L, 2557870080L, 1438842880L, 2330034176L, 676380672L, 2637832192L, 1724854272L, 3872225280L, 2798488576L, 2261603840L, 2529962752L, 248489344L, 1523447360L, 3503185760L, 4191042160L, 1699639256L, 51167244L, 3828097030L, 2147483648L, 1073741824L, 2684354560L, 805306368L, 134217728L, 2885681152L, 2785017856L, 788529152L, 1115684864L, 2218786816L, 4225761280L, 3788505088L, 2423783424L, 960233472L, 2769682432L, 177274880L, 2300215296L, 1823260672L, 3331203072L, 2122297344L, 443947008L, 1879145472L, 2818812416L, 2617270528L, 2919323520L, 2197916480L, 3833713632L, 2873248240L, 3110105560L, 1704049668L, 1805355530L, 2147483648L, 3221225472L, 536870912L, 2952790016L, 2818572288L, 1006632960L, 1979711488L, 3607101440L, 662700032L, 4013948928L, 48234496L, 202375168L, 2677538816L, 2323382272L, 2180382720L, 1074987008L, 1637384192L, 3514744832L, 2020909056L, 1146351616L, 855201792L, 3844971520L, 3263970816L, 763120896L, 779572096L, 592620608L, 3178807392L, 913286384L, 3060757592L, 4134378508L, 2523676162L, 2147483648L, 3221225472L, 1610612736L, 4026531840L, 1207959552L, 603979776L, 2583691264L, 486539264L, 3061841920L, 2202009600L, 3760193536L, 2985295872L, 1754791936L, 4108058624L, 37617664L, 2704211968L, 2426896384L, 3102785536L, 1813553152L, 3198111744L, 2258454528L, 2862605312L, 882848256L, 1671346432L, 1344273024L, 3655113408L, 2629757664L, 4134351024L, 2726803608L, 805535436L, 671209702L, 2147483648L, 3221225472L, 2684354560L, 4026531840L, 3892314112L, 2751463424L, 2717908992L, 1929379840L, 729808896L, 2252341248L, 2422210560L, 3096444928L, 3165126656L, 4002676736L, 1957560320L, 4220715008L, 3744104448L, 1020772352L, 2932711424L, 2483130368L, 3925870592L, 654445568L, 1635978752L, 1363183872L, 434162048L, 1305685568L, 122313888L, 3494224720L, 1483051832L, 2912291852L, 388828682L, 2147483648L, 3221225472L, 2684354560L, 4026531840L, 3892314112L, 335544320L, 436207616L, 2566914048L, 3162505216L, 4022337536L, 2468347904L, 3698327552L, 3197632512L, 2332819456L, 546177024L, 2966355968L, 162299904L, 96256000L, 3820429312L, 1970311168L, 1243170816L, 2160980992L, 1102002688L, 3787361024L, 297031296L, 4168248000L, 3969606880L, 4128076528L, 1846410456L, 3524176588L, 1007437034L, 2147483648L, 3221225472L, 3758096384L, 4026531840L, 939524096L, 4093640704L, 1174405120L, 520093696L, 1535115264L, 2227175424L, 1063256064L, 1244659712L, 1288175616L, 4087087104L, 4192600064L, 368246784L, 3056304128L, 644857856L, 2949570560L, 3281833984L, 568629248L, 300989440L, 3358199296L, 3446539008L, 2991523200L, 1482523840L, 1157940512L, 3741409648L, 3145798680L, 1958968332L, 123923982L, }; 
    
    
    final int maxBits = 31;
    float Sobol( int n2, int D, int scramble) {
      long n = n2;
      long r = (long)scramble;
      for (int idx = D*maxBits; n !=0 ; idx++)
      {
        if ((n & 0x1L) != 0L ) 
          r ^= SobolMatrices[idx];
        n>>=1;
      }
       r  = (r &0xFFFFFFFFL);
        
      return ((float)r)/  ( (float)0x100000000L); }
    
    class QMCSamples
    {
      int ranX;
      int ranY;
      QMCSamples( int rx, int ry)
      {
        ranX = rx;
        ranY = ry;
      }
      void Next()
      {
        ranX *=16807;
        ranY *=16807;
      }
      float Getf(int sc,int D)
      {
        return  Sobol(sc, D , ranX);
      }
      void Set( int x, int y )
      {
        ranX = x | (x<<16);
        ranY = y | (y<<16);
        ranX *=16807;
        ranY *=16807;
      }
      float[] Get( int sc, int D )
       {
         float samples[] = new float[2];
        samples[0]  = Sobol(sc, D*2 , ranX);
        samples[1]  =  Sobol(sc, D*2+1 , ranY);
        return samples;
       }
       float[] GetS( float samples[], int sc, int D )
       {
        samples[0]  = Sobol(sc, D*2 , ranX);
        samples[1]  =  Sobol(sc, D*2+1 , ranY);
        return samples;
       }
    }
    
    
    class hitResult
    {
      float t = 1000000.f;
      int tri;
      float u;
      float v;
    };
    
    // from Walds thesis
    final int modulo3[] = {
      0,1,2,0,1
    };
    
    class TriAccel
    {
      // first 16 byte half cache line
      // plane:
      float n_u; //!< == normal.u / normal.k
      float n_v; //!< == normal.v / normal.k
      float n_d; //!< constant of plane equation
      int k; // projection dimension
      // second 16 byte half cache line
      // line equation for line ac
      float b_nu;
      float b_nv;
      float b_d;
      int triNum; // pad to next cache line
      // third 16 byte half cache line
      // line equation for line ab
      float c_nu;
      float c_nv;
      float c_d;
      int pad1; // pad to 48 bytes for cache alignment purposes
    
      // note could pack 3 vertex indices and material id
    
      TriAccel( Tri t, int num)
      {
        Gen(t.v0, t.v1,t.v2, num);
      }
      TriAccel( PVector A, PVector B, PVector C, int num)
      {
        Gen(A,B,C,num);
      }
      void Gen( PVector A, PVector B, PVector C, int num)
      {
        // calc normal
        PVector BT = PVector.sub(C,A);
        PVector CT = PVector.sub(B,A);
        PVector  N = CT.cross(BT);
        float[] b = {
          BT.x,BT.y,BT.z
        };
        float[] c = {
          CT.x,CT.y,CT.z
        };
    
        if ( abs(N.x) > abs(N.y))
          if (abs(N.x) > abs(N.z))
            k = 0; // X
          else
            k=2; // Z
        else
          if (abs(N.y) > abs(N.z))
            k = 1; // Y
          else
            k=2; // Z
        int u = (k+1) %3;
        int v = (k+2) % 3;
        float[] n = {
          N.x,N.y,N.z
        };
    
        float[] a = {
          A.x,A.y,A.z
        };
    
        if ( n[k] == 0.0)
        {
          println("Degenerate Tri " + A + B + C + N);
        }
        assert( abs(n[k])>0.00000f);
    
        n_u =n[u]/n[k]; //!< == normal.u / normal.k
        n_v = n[v]/n[k]; //!< == normal.v / normal.k
        n_d = a[u]*n_u + a[v]*n_v + a[k]; //!< constant of plane equation
    
        // do B
        float bdenom = (b[u]*c[v]-b[v]*c[u]);
        assert( abs(bdenom)>0.00000f);
        bdenom = 1./bdenom;
    
        b_nu = -b[v]*bdenom;
        b_nv = b[u]*bdenom;
        b_d = (b[v]*a[u]-b[u]*a[v])*bdenom;
    
        // do C
        c_nu = c[v]*bdenom;
        c_nv = -c[u]*bdenom;
        c_d = (-c[v]*a[u]+c[u]*a[v])*bdenom;
        triNum = num;
      }
     /* PVector GetNormal( PVector I)
       {
         int ku = modulo3[k+1];
        int kv = modulo3[k+2];
        float[] nv = new float[3];
        nv[ku]= n_u;
        nv[kv]= n_v;
        nv[k]=.1;
        
        // normalize nv
        
         PVector n = new PVector(uv[0],nv[1],nv[2]);
         n.set(n0);
         n.mult( I.dot(n)>0 ? -1 :1);
         return n;
      }*/
      hitResult Intersect(hitResult hit, float[] org, float[] dir, float maxT )
      {
        int ku = modulo3[k+1];
        int kv = modulo3[k+2];
        //hitResult hit = new hitResult();
        hit.t = 1000000.f;
    
        // don’t prefetch here, assume data has already been prefetched
        // start high-latency division as early as possible
        final float nd = 1./(dir[k] +  n_u * dir[ku] + n_v * dir[kv]);
        final float f = (n_d - org[k] - n_u * org[ku] - n_v * org[kv]) * nd;
    
        // check for valid distance.
        if (!(maxT > f && f > 0.0001f ))
          return hit;
        // compute hitpoint positions on uv plane
        final float hu = org[ku] + f * dir[ku];
        final float hv = org[kv] + f * dir[kv];
    
        // check first barycentric coordinate
        final float lambda = hu * b_nu + hv * b_nv + b_d;
        if (lambda < 0.0f)
          return hit;
    
        // check second barycentric coordinate
        final float mue = hu * c_nu + hv * c_nv + c_d;
        if (mue < 0.0f || lambda+mue > 1.0f) 
          return hit;
        // have a valid hitpoint here. store it.
        hit.t = f;
        hit.tri = triNum;
        hit.u = lambda;
        hit.v = mue;
        return hit;
      }
    };
    
    
    class Node
    {
      float minx;
      float maxx;
      float miny;
      float maxy;
      float minz;
      float maxz;
      int   childIdx;
      int   leafStart;
      int   leafEnd;
    };
    
    float max4( float a, float b, float c, float d) {
      return max(max(a,b),max(c,d));
    }
    float min4( float a, float b, float c, float d) {
      return min(min(a,b),min(c,d));
    }
    hitResult TraverseBVH( hitResult hit, int[] traversalStack, float[] dir, Node[] nodes, TriAccel[] Tris, boolean isShadow, float maxT, float origx,float origy,float origz, float idirx,float idiry,float idirz) {
      float tmin = 0.000001;
      hit.t = maxT;
      final int EntrypointSentinel = 0x7F00FF00;
      traversalStack[0]=EntrypointSentinel;
      // Intersect the ray against the child nodes.
      float oodx  = origx * idirx;
      float oody  = origy * idiry;
      float oodz  = origz * idirz;
    
      int stackPtr =0;
      int nodeAddr = 0;
    
      // maybe grab 16 rays here?
      while(true)
      {
        while(nodeAddr >= 0 && nodeAddr != EntrypointSentinel )
        {
          // Fetch AABBs of the two child nodes.
          Node n0 = nodes[nodeAddr];
          Node n1 = nodes[nodeAddr+1];
    
          float c0lox = n0.minx * idirx - oodx;
          float c0hix = n0.maxx * idirx - oodx;
    
          float c0loy = n0.miny * idiry - oody;
          float c0hiy = n0.maxy* idiry - oody;
    
          float c0loz = n0.minz  * idirz - oodz;
          float c0hiz = n0.maxz  * idirz - oodz;
    
          float c1loz = n1.minz  * idirz - oodz;
          float c1hiz = n1.maxz  * idirz - oodz;
    
          float c0min = max4( min(c0lox, c0hix), min(c0loy, c0hiy), min(c0loz, c0hiz), tmin);
          float c0max = min4( max(c0lox, c0hix), max(c0loy, c0hiy), max(c0loz, c0hiz), maxT);
    
          float c1lox = n1.minx  * idirx - oodx;
          float c1hix = n1.maxx * idirx - oodx;
    
          float c1loy = n1.miny  * idiry - oody;
          float c1hiy = n1.maxy * idiry - oody;
    
          float c1min = max4(min(c1lox, c1hix), min(c1loy, c1hiy), min(c1loz, c1hiz), tmin);
          float c1max = min4(max(c1lox, c1hix), max(c1loy, c1hiy), max(c1loz, c1hiz), maxT);
    
          // Decide where to go next.
          boolean traverseChild0 = (c0max >= c0min);
          boolean traverseChild1 = (c1max >= c1min);
    
          nodeAddr           = n0.childIdx;      // stored as int
          int nodeAddrChild1 = n1.childIdx;      // stored as int
    
    
          // One child was intersected => go there.
          if(traverseChild0 != traverseChild1)
          {
            if(traverseChild1)
              nodeAddr = nodeAddrChild1;
          }
          else
          {
            // Neither child was intersected => pop.
            if (!traverseChild0)
            {
              nodeAddr = traversalStack[stackPtr];
              --stackPtr;
            }
            // Both children were intersected => push the farther one.
            else
            {
              if(c1min < c0min)
              {
                int temp = nodeAddr;
                nodeAddr = nodeAddrChild1;
                nodeAddrChild1 = temp;//            swap(nodeAddr, nodeAddrChild1);
              }
              ++stackPtr;
              traversalStack[stackPtr] = nodeAddrChild1;
            }
          }
        }
    
        if ( nodeAddr == EntrypointSentinel )
          return hit;
    
        // do triangle tests
        int lidx = -nodeAddr-1;
        int triAddr  = nodes[lidx].leafStart;                  // stored as int
        int triAddr2 = nodes[lidx].leafEnd;                  // stored as int
    
    
        nodeAddr = traversalStack[stackPtr];  // wonder if doing a while nodeAddr <0 here would be nice
        --stackPtr;
    
        float[] org  = { 
          origx, origy,origz
        };
        while( triAddr < triAddr2)
        {
          TriAccel acc = Tris[triAddr++];
          int ku = modulo3[acc.k+1];
          int kv = modulo3[acc.k+2];
    
          // don’t prefetch here, assume data has already been prefetched
          // start high-latency division as early as possible
          final float nd = 1./(dir[acc.k] +  acc.n_u * dir[ku] + acc.n_v * dir[kv]);
          final float f = (acc.n_d - org[acc.k] - acc.n_u * org[ku] - acc.n_v * org[kv]) * nd;
    
          // check for valid distance.
          if ((maxT > f && f > 0.0001f ))
          {
            // compute hitpoint positions on uv plane
            final float hu = org[ku] + f * dir[ku];
            final float hv = org[kv] + f * dir[kv];
    
            // check first barycentric coordinate
            final float u = hu * acc.b_nu + hv * acc.b_nv + acc.b_d;
            if (u > 0.0f)
            {
              // check second barycentric coordinate
              final float v = hu * acc.c_nu + hv * acc.c_nv + acc.c_d;
              if (v > 0.0f && u+v < 1.0f) 
              {
                maxT = f;
    
                // have a valid hitpoint here. store it.
                hit.t = f;
                hit.tri = acc.triNum;
                hit.u = u;
                hit.v = v;
    
                if ( isShadow )
                {
                  return hit;
                }
              }
            }
          }
        }
      }
      // done tri
    }
    
    class TraversalTree {
      TriAccel[] accTris;
      Tri[]    triList;
      Node[]    bvh;
      Luminare[] luminares=null;
    
      void AddLuminare( PVector a, PVector b, PVector c, PVector d)
      {
        if (luminares == null || luminares.length==0)
            luminares = new Luminare[1];
        luminares[0]= new Luminare(a,b,c,d);
      }
      PVector Luminance( QMCSamples sampler, int sc, RenderState rs, 
                PVector pt, PVector nrm, PVector lcol, int D )
      {
        //assert (luminares.length ==1);
        //assert (numSamples ==1);
        float[] origin= {
          pt.x,pt.y,pt.z
        };
        PVector res=new PVector(0.,0.,0.);
        for (int i =0; i < luminares.length;i++)
        {
          rs.samples2=sampler.GetS(rs.samples2, sc,D);
          PVector lpt = luminares[i].Sample( rs.samples2[0], rs.samples2[1]); // could select number of shadow rays depending on importance
          PVector dt = PVector.sub(lpt,pt);
          float d = dt.mag();
          dt.mult(1./d);
    
          float[] dir= {
            dt.x,dt.y,dt.z
          };
          float lum =max(dt.dot(nrm),0);
          // see if in shadow
          if (lum > 0 && TraceShadow(rs.hit, rs.traversalStack, origin,dir, d))
          {
            // apply light could store it in a list of norml, light, pairs
            float atten=luminares[i].area/(d*d);
            lum*=atten;     
            res.add(PVector.mult(lcol,lum));
          }
        }
        return res;
      }
    
      boolean TraceShadow(hitResult hit, int[] traversalStack, float[] origin, float[] dir, float maxT )
      {
        float[] dv = { 
          dir[0] ==0. ? 10000.f :  1./dir[0], dir[1] ==0. ? 10000.f :  1./dir[1], dir[2] ==0. ? 10000.f : 1./dir[2]
        };
        hit= TraverseBVH( hit, traversalStack, dir,bvh, accTris, true,maxT, origin[0],origin[1],origin[2],
        dv[0],dv[1],dv[2]);
        return hit.t>=maxT;
      }
      hitResult Trace(RenderState rs, float[] origin, float[] dir, float maxT )
      {
        float[] dv = { 
          dir[0] ==0. ? 10000.f :  1./dir[0], dir[1] ==0. ? 10000.f :  1./dir[1], dir[2] ==0. ? 10000.f : 1./dir[2]
        };
        return TraverseBVH( rs.hit, rs.traversalStack, dir,bvh, accTris, false,maxT, origin[0],origin[1],origin[2],
        dv[0],dv[1],dv[2]);
      }
    };
    
    TraversalTree BuildBVHTree( Tri[] triList, boolean addGroundPlane, float zflip, float gplaneHieght) {
      println("Build Tree");
    
      if ( zflip != 0.)
      {
        for(int i =0; i < triList.length;i++)
        {
          Tri t = triList[i];
          triList[i] = new Tri( new PVector(t.v0.x,t.v0.z, t.v0.y),  
          new PVector(t.v1.x,t.v1.z, t.v1.y),
          new PVector(t.v2.x,t.v2.z, t.v2.y),
          t.sn0,t.sn1,t.sn2);
        }
      }  
      if ( addGroundPlane )
      {
        float s=2.;
        triList =(Tri[]) append(triList, new Tri(new PVector(-s,gplaneHieght,-s), new PVector(s,gplaneHieght,s) ,new PVector(s,gplaneHieght,-s)));
        triList =(Tri[]) append(triList, new Tri(new PVector(s,gplaneHieght,s), new PVector(-s,gplaneHieght,-s),new PVector(-s,gplaneHieght,s)));
      }
    
      TraversalTree ttree = new TraversalTree();
      ttree.accTris = new TriAccel[triList.length];
      ttree.bvh= new Node[triList.length*2];
      ttree.triList = triList;
    
      // build a tree
      // get log2 of trilist.length
      int l = triList.length;
      int ml=0;
      for (; l !=0; ml++,l>>=1);
      println("Max Level "+ml);
    
      int maxLevel = ml;
      BuildState bs = new BuildState( ttree.triList);
      nodeInfo rt = BuildTree( bs, bs.m_totalBound, ttree.bvh, 0, 0,ttree.triList, ttree.accTris,0, ttree.accTris.length,0, maxLevel );
    
      /* println("Tree --- ");
       DumpTreeInfo(ttree.bvh, 0,0);
       DumpTreeInfo(ttree.bvh, 1,0);
       println("number of nodes "+rt.nidx);
       */
      return ttree;
    }
    
    TraversalTree GenerateScene()
    {
      int amt=128*1024;
      Tri[] triList = new Tri[amt];
      randomSeed(10);
      for (int i =0; i < triList.length; i++)
      {
        PVector c = new PVector(random(-1.,1.),random(-1.,1.),random(-1.,1.));
        float r = 1./sqrt((float)amt);
        triList[i] = new Tri( new PVector(c.x + random(-r,r),c.y + random(-r,r),c.z + random(-r,r)),
        new PVector(c.x + random(-r,r),c.y + random(-r,r),c.z + random(-r,r)),
        new PVector(c.x + random(-r,r),c.y + random(-r,r),c.z + random(-r,r)));
      }
    
      return BuildBVHTree( triList, true,0.,-1.);
    }
    
    void LoadScene( Scene scene )
    {
       OBJReader ObjFile = new OBJReader( scene.n,false);
      ObjFile.readFile();
      g_camControl = new PolarCamera(4., scene.cangX,-0.3,1.);
      Tri[] triList = ObjFile.GenerateTriList();
      if ( scene.rescale != 1.)
        for (int i=0; i < triList.length;i++)
        {
          Tri t = triList[i];
          triList[i] = new Tri( PVector.mult(triList[i].v0, scene.rescale),
                          PVector.mult(triList[i].v1, scene.rescale),
                          PVector.mult(triList[i].v2, scene.rescale),
                          t.sn0,t.sn1,t.sn2);
        }
       g_NumSamplesRequired = scene.numRequireSamples;
      //ttree=GenerateScene();
      Luminare[] luminares = ObjFile.GetLuminares();
      ObjFile=null;
      Runtime r = Runtime.getRuntime(); 
       r.gc();
       
      ttree = BuildBVHTree( triList, true, scene.zflip,scene.gplane);
      ttree.luminares = luminares;
      if ( scene.rescale != 1.) {
       Luminare[] lums = ttree.luminares;
        for (int i=0; i < ttree.luminares.length;i++)
        {
          lums[i] =new Luminare( PVector.mult(lums[i].v0, scene.rescale),
                          PVector.mult(lums[i].v1, scene.rescale),
                          PVector.mult(lums[i].v2, scene.rescale),
                          PVector.mult(lums[i].v3, scene.rescale));
        }
      }
      if ( scene.addLight)
        AddLight();
    
      r.gc();
      println("At program start we have : " + r.freeMemory() + " bytes");
      
      g_useReflect=scene.useReflect;
      if ( scene.colorSel !=-1)
        g_colorSelection=scene.colorSel;
    } 
    float g_lw = 0.5;//2.;
    float g_lh=1.4;
    float g_ly=.7;
    float g_langX = 0.;
    float g_langY = 1.;
        
    void AddLight()
    {  
      PVector off = new PVector();
      off.x = sin(g_langX)*cos(g_langY);
      off.z = cos(g_langX)*cos(g_langY);
      off.y = sin(g_langY);
      ttree.AddLuminare( PVector.add( off, new PVector(-g_lw,g_ly,-g_lh)), 
       PVector.add( off, new PVector(g_lw,g_ly,-g_lh)),
      PVector.add( off, new PVector(g_lw,g_ly,g_lh)),
      PVector.add( off, new PVector(-g_lw,g_ly,g_lh))
     );
    }
    
    public class WorkQueue
    {
        private final int nThreads;
        private final PoolWorker[] threads;
        private final LinkedList queue;
     
        public WorkQueue(int nThreads)
        {
            this.nThreads = nThreads;
            queue = new LinkedList();
            threads = new PoolWorker[nThreads];
     
            for (int i=0; i<nThreads; i++) {
                threads[i] = new PoolWorker();
                threads[i].start();
            }
        }
        public void clear(){
          synchronized(queue) {
            queue.clear();
          }
        }
         int NumWorkItems()
         {
           return queue.size();
         }
        public boolean IsFinished() {
            for (int i=0; i<threads.length; i++)
              if ( threads[i].working )
                return false;
            return queue.isEmpty();
        }
        public void execute(Runnable r) {
            synchronized(queue) {
                queue.addLast(r);
                queue.notify();
            }
        }
     
        private class PoolWorker extends Thread {
            boolean working;
                 
            public void run() {
                working = true;
                Runnable r;
                 
                while (true) {
                    synchronized(queue) {
                        while (queue.isEmpty()) {
                            working = false;
                            try {                           
                                queue.wait();
                            }
                            catch (InterruptedException ignored)
                            {
                            }
                        }                   
                        working = true;
                        r = (Runnable) queue.removeFirst();
                    }
                    // If we don't catch RuntimeException,
                    // the pool could leak threads
                    try {
                        r.run();
                    }
                    catch (RuntimeException e) {
                        // You might want to log something here
                        println("TASK ERROR");
                    }
                }
            }
        }
    }
    
    

    code

    tweaks (0)

    about this sketch

    This sketch is running as Java applet, exported from Processing.

    license

    advertisement

    edward Porten

    Real-Time Path Tracing Obj Viewer

    Add to Faves Me Likey@! 33
    You must login/register to add this sketch to your favorites.

    Loads a series of obj models ( based from ahmet.kizilay obj loader)
    and renders then using path tracing.
    Path tracer uses progressive spatial caching to get real-time
    update and views of model.
    Keys
    '+' go to next model.
    Drag mouse to view around scene
    Drag and hold shift to zoom/ pan
    Drag and hold control to resize area light
    Drag and hold alt to move area light.
    'D' toggle direct lighting.
    'R' toggle reflection
    'C' show spatial caching
    'S' change model colors

    Felix Woitzel
    31 Oct 2010
    Really nice man! At first glance i thought this would be a T-Rex. ;D
    edward Porten
    2 Nov 2010
    Uploaded a new version with improved reflection, better mult-core support and a T-Rex. The t-rex needs a fast machine, so it's the third model ( press + twice)
    Felix Woitzel
    2 Nov 2010
    I've got a 2.66GHz Quadcore and it runs fine here. The dino is awesome! :D
    What about some normal interpolation.
    Regards!
    edward Porten
    2 Nov 2010
    I was thinking of that but most of the models are from the mesh compendieum and don't have normals supplied. I could add code to auto generate the normals I suppose.
    edward Porten
    2 Nov 2010
    Oh yeah, and the dino model is brilliant. It really suits ray tracing as it a lot of rays miss it.
    edward Porten
    5 Nov 2010
    Uploaded new version with normal interpolation. I imported the models into blender and exported it with smooth normals. Thanks for the idea Felix, it looks a lot better.
    great sketch! would it be possible to use it for the geometry generated in processing?
    edward Porten
    6 Dec 2010
    The function BuildBVHTree( triList, true, scene.zflip,scene.gplane) builds a bounding volume of the triangles, this is slow to do per frame and also it uses a spatial cache so that the lighting is calculated over numerous frames, so if the geometry changes it would need to be flushed, so for dynamic stuff it wouldn't work that well sadly.
    thx. i meant rather creating geometry, freezing, rendering -looking at it, going back... could be nice for all the architecture-related processing users.
    DJLinux007
    18 Feb 2011
    good job
    d.1
    27 Dec 2013
    Hello i have imported an .obj file ..is a 3d molecule and i just cant move the smoothly the molecule...do you know why?
    You need to login/register to comment.