source: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp @ 2635

#include "RayCaster.h"
#include "VssRay.h"
#include "Ray.h"
#include "Preprocessor.h"
#include "ViewCellsManager.h"

#include <cassert>

namespace GtpVisibilityPreprocessor {


#define DEBUG_PROCESS_RAY 0

#define DEBUG_RAYCAST 0

#define EXACT_BOX_CLIPPING 0

RayCaster::RayCaster(const Preprocessor &preprocessor):
  mVssRayPool(), mPreprocessor(preprocessor)
{
}


RayCaster::~RayCaster()
{
}


VssRay *RayCaster::CastRay(const SimpleRay &simpleRay,
                                                   const AxisAlignedBox3 &box,
                                                   const bool castDoubleRay) 
{
        VssRayContainer rays;
        CastRay(simpleRay, rays, box, castDoubleRay, true);
    
        if (!rays.empty())
                return rays.back();
        else
                return NULL;
}


bool
RayCaster::ClipToViewSpaceBox(const Vector3 &origin,
                                                          const Vector3 &termination,
                                                          Vector3 &clippedOrigin,
                                                          Vector3 &clippedTermination)
{
  
  Ray ray(origin, termination - origin, Ray::LINE_SEGMENT);     
  ray.Precompute();
  
  float tmin, tmax;
  if ((!mPreprocessor.mViewCellsManager->
           GetViewSpaceBox().ComputeMinMaxT(ray, &tmin, &tmax)) ||
          tmin>=tmax
          )
        return false;

  if (tmin >= 1.0f || tmax <= 0.0f)
        return false;
  
  if (tmin > 0.0f)
        clippedOrigin = ray.Extrap(tmin);
  else
        clippedOrigin = origin;
  
  if (tmax < 1.0f)
        clippedTermination = ray.Extrap(tmax);
  else
        clippedTermination = termination;
  
  return true;
}

/** Checks if ray is valid
        (e.g., not in empty view space or outside the view space)
*/
bool
RayCaster::ValidateRay(const Vector3 &origin,
                                           const Vector3 &direction,
                                           const AxisAlignedBox3 &box,
                                           Intersection &hit)
{
        if (!hit.mObject)  
        {
                // compute intersection with the scene bounding box
#if EXACT_BOX_CLIPPING
                static Ray ray;
                mPreprocessor.SetupRay(ray, origin, direction);

                float tmin, tmax;
                if (box.ComputeMinMaxT(ray, &tmin, &tmax) && (tmin < tmax)) 
                {
                        hit.mPoint = ray.Extrap(tmax);
                }
                else 
                {
                        // cout<<" invalid hp "<<tmin<<" "<<tmax<<endl;
                        // cout<<" box "<<box<<endl;
                        // cout<<" origin "<<origin<<endl;
                        // cout<<" direction "<<direction<<endl;
                        // cout<< "inv dir"<<ray.GetInvDir()<<endl;
                        return false;
                }
#else
                hit.mPoint = origin + direction * Magnitude(box.Diagonal());
#endif
        }
        else
        {
          if (mPreprocessor.mDetectEmptyViewSpace)  {
                if (DotProd(hit.mNormal, direction) >= -Limits::Small)  {       
                  hit.mObject = NULL;
                  return false;
                }
          }
        }
        
        return true;
}

void
RayCaster::SortRays(SimpleRayContainer &rays)
{
  AxisAlignedBox3 box =
        mPreprocessor.mViewCellsManager->GetViewSpaceBox();
  
  float b[12]={
        box.Min().x,
        box.Min().y,
        box.Min().z,
        -1,
        -1,
        -1,
        box.Max().x,
        box.Max().y,
        box.Max().z,
        1,
        1,
        1
  };

#if 0
  static vector<SimpleRay *> pointerArray;

  if (pointerArray.size()!=rays.size()) {
        // realloc the pointerarray
        pointerArray.resize(rays.size());
  }

  // init pointer array
  SimpleRay *p = &pointerArray[0];
  for (i=0; i < rays.size(); i++, p++)
        pointerArray[i] = p;
#endif
  
  _SortRays(rays,
                        0,
                        (int)rays.size()-1,
                        0,
                        b
                        );
}
                                        
void
RayCaster::_SortRays(SimpleRayContainer &rays,
                                         const int l,
                                         const int r,
                                         const int depth,
                                         float box[12])
{
  // pick-up a pivot
  int axis;
  
  float maxDiff = -1.0f;
  // get the largest axis
  int offset = 0;
  int i;

  //const int batchsize = 16384;
  const int batchsize = 8192;
  //const int batchsize = 128;

  //if (r - l < 16*batchsize)
  //    offset = 3;
        
  //  if (depth%2==0)
  //    offset = 3;
  
  for (i=offset; i < offset + 3; i++) {
        float diff = box[i + 6] - box[i];
        if (diff > maxDiff) {
          maxDiff = diff;
          axis = i;
        }
  }

  //  cout<<depth<<" "<<axis<<" "<<l<<" "<<r<<endl;
  
  i=l;
  int j=r;

  float x = (box[axis] + box[axis+6])*0.5f;
  //  float x = rays[(l+r)/2].GetParam(axis);
  do {
        while(i<j && rays[i].GetParam(axis) < x)
          i++;
        while(i<j && x < rays[j].GetParam(axis))
          j--;
        
        if (i <= j) {
          swap(rays[i], rays[j]);
          i++;
          j--;
        }
  } while (i<=j);

  
  if (l + batchsize < j ) {
        // set new max
        float save = box[axis+6];
        box[axis+6] = x;
        _SortRays(rays, l, j, depth+1, box);
        box[axis+6] = save;
  } else {
        //      for (int k=0; k < 6; k++)
        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
  }
  
  if (i + batchsize < r) {
        // set new min
        box[axis] = x;
    _SortRays(rays, i, r, depth+1, box);
  } else {
        //      for (int k=0; k < 6; k++)
        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
  }
        
}
  

void
RayCaster::SortRays2(SimpleRayContainer &rays)
{
  AxisAlignedBox3 box =
        mPreprocessor.mViewCellsManager->GetViewSpaceBox();

  const float sizeBox = Magnitude(box.Diagonal());
  // This is some size of the 
  const float sizeDir = 0.2f * sizeBox;
  
  float b[12]={
    box.Min().x,
    box.Min().y,
    box.Min().z,
    -sizeDir,
    -sizeDir,
    -sizeDir,
    box.Max().x,
    box.Max().y,
    box.Max().z,
    sizeDir,
    sizeDir,
    sizeDir
  };

#if 0
  static vector<SimpleRay *> pointerArray;

  if (pointerArray.size()!=rays.size()) {
        // realloc the pointerarray
        pointerArray.resize(rays.size());
  }

  // init pointer array
  SimpleRay *p = &pointerArray[0];
  for (i=0; i < rays.size(); i++, p++)
        pointerArray[i] = p;
#endif
  
  _SortRays2(rays, 0, (int)rays.size()-1, 0, b);

  return;
}
                                        
void
RayCaster::_SortRays2(SimpleRayContainer &rays,
                      const int l,
                      const int r,
                      const int depth,
                      float box[12])
{
  // pick-up a pivot
  int axis;
  
  float maxDiff = -1.0f;
  // get the largest axis
  int offset = 0;
  int i;

  //const int batchsize = 16384;
  //const int batchsize = 8192;
  const int batchsize = 128;

  //if (r - l < 16*batchsize)
  //    offset = 3;
        
  //  if (depth%2==0)
  //    offset = 3;
  
  for (i=offset; i < offset + 6; i++) {
    float diff = box[i + 6] - box[i];
    assert(diff >= 0.f);
    if (diff > maxDiff) {
      // Find the maximum
      maxDiff = diff;
      axis = i;
    }
  }

  //  cout<<depth<<" "<<axis<<" "<<l<<" "<<r<<endl;
  
  i=l;
  int j=r;

  float x = (box[axis] + box[axis+6])*0.5f;
  //  float x = rays[(l+r)/2].GetParam(axis);
  do {
        while(i<j && rays[i].GetParam(axis) < x)
          i++;
        while(i<j && x < rays[j].GetParam(axis))
          j--;
        
        if (i <= j) {
          swap(rays[i], rays[j]);
          i++;
          j--;
        }
  } while (i<=j);

  
  if (l + batchsize < j ) {
        // set new max
        float save = box[axis+6];
        box[axis+6] = x;
        _SortRays2(rays, l, j, depth+1, box);
        box[axis+6] = save;
  } else {
        //      for (int k=0; k < 6; k++)
        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
  }
  
  if (i + batchsize < r) {
    // set new min
    box[axis] = x;
    _SortRays2(rays, i, r, depth+1, box);
  } else {
        //      for (int k=0; k < 6; k++)
        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
  }
        
}
  
  
VssRay *RayCaster::RequestRay(const Vector3 &origin, 
                                                          const Vector3 &termination, 
                                                          Intersectable *originObject, 
                                                          Intersectable *terminationObject, 
                                                          const int pass, 
                                                          const float pdf)
{
 #if DEBUG_RAYCAST
        Debug<<"PR2a"<<flush;
#endif
        
        // old method: always allocate
        if (0) return new VssRay(origin, termination, originObject, terminationObject, pass, pdf);

        VssRay *vssRay = mVssRayPool.Alloc();

#if DEBUG_RAYCAST
        Debug<<"PR2b"<<flush;
#endif
        
        *vssRay = VssRay(origin, termination, originObject, terminationObject, pass, pdf);

#if DEBUG_RAYCAST
        Debug<<"PR2c"<<flush;
#endif

        return vssRay;
}


int
RayCaster::ProcessRay(const SimpleRay &simpleRay,
                                          Intersection &hitA,
                                          Intersection &hitB,
                                          VssRayContainer &vssRays,
                                          const AxisAlignedBox3 &box,
                                          const bool castDoubleRay,
                                          const bool pruneInvalidRays)
{
  int hits = 0;


#if DEBUG_RAYCAST
  static int id=0;
  Debug<<"PRA "<<id++<<endl<<flush;
#endif
        
        if (pruneInvalidRays) 
        {
          if (!hitA.mObject && !hitB.mObject) {
#if DEBUG_PROCESS_RAY
                cout<<"I1 ";
#endif
                return 0;
          }
        }
        
        // regardless of the pruneInvalidRays setting reject 
        // rays whic degenerate to a point
        if (EpsilonEqualV3(hitA.mPoint, hitB.mPoint, Limits::Small)) {
#if DEBUG_PROCESS_RAY
                cout<<"I2 ";
#endif
          return 0;
        }
        
        const bool validA = ValidateRay(simpleRay.mOrigin, simpleRay.mDirection, box, hitA);
        const bool validB = //castDoubleRay && 
          ValidateRay(simpleRay.mOrigin, -simpleRay.mDirection, box, hitB);
        
        
#if DEBUG_RAYCAST
        Debug<<"PR1"<<flush;
#endif
        
        // reset both contributions
        if (!validA || !validB) {
          if (0 || pruneInvalidRays)
                return 0;

#if DEBUG_PROCESS_RAY
                cout<<"I2 ";
#endif

          // reset both contributions of this ray
          hitA.mObject = NULL;
          hitB.mObject = NULL;
        }
        
        // 8.11. 2007 JB
        // degenerate rays checked by geometrical constraint...
        //      !pruneInvalidRays || (hitA.mObject != hitB.mObject);

        
#if DEBUG_RAYCAST
        Debug<<"PR2"<<flush;
#endif

        // VH - I do not know what this is for, commented out temporarily
        // !!!!!!!!!!!!!! !!!!!!!!!!!! !!!!!!!!!!
#if 1
        const bool validSample = true;
        if (validSample) {
          Vector3 clipA, clipB;
          if (!ClipToViewSpaceBox(hitA.mPoint,
                                                          hitB.mPoint,
                                                          clipA,
                                                          clipB)) {
#if DEBUG_PROCESS_RAY
                cout<<"I3 ";
#endif

                return 0;
          }

          if (!pruneInvalidRays || hitA.mObject) {

                  VssRay *vssRay = 
                          RequestRay(!castDoubleRay ? simpleRay.mOrigin : clipB,
                                                 hitA.mPoint,
                                                 hitB.mObject,
                                                 hitA.mObject,
                                                 mPreprocessor.mPass,
                                                 1.0f //simpleRay.mPdf
                                                 );

                if (validA)
                  vssRay->mFlags |= VssRay::Valid;
                
                vssRay->mDistribution = simpleRay.mDistribution;
                vssRay->mGeneratorId = simpleRay.mGeneratorId;

                vssRays.push_back(vssRay);
                ++ hits;
                //cout << "vssray 1: " << *vssRay << " " << vssRay->mTermination - vssRay->mOrigin << endl;
        }

#if DEBUG_RAYCAST
        Debug<<"PR3"<<flush;
#endif

        if (castDoubleRay && (!pruneInvalidRays || hitB.mObject))
        {
                VssRay *vssRay = RequestRay(
                                            clipA,
                                            hitB.mPoint,
                                            hitA.mObject,
                                            hitB.mObject,
                                            mPreprocessor.mPass,
                                            1.0f //simpleRay.mPdf
                                            );

                if (validB)
                        vssRay->mFlags |= VssRay::Valid;

                vssRay->mDistribution = simpleRay.mDistribution;
                vssRay->mGeneratorId = simpleRay.mGeneratorId;
                vssRays.push_back(vssRay);
                ++ hits;
                //cout << "vssray 2: " << *vssRay << endl;
        }
#if DEBUG_RAYCAST
          Debug<<"PR4"<<flush;
#endif
        } // validSample
#else
        // Just pass if intersected or not
        hits = (hitA.mObject != 0) ? 1 : 0;
        intersect = hitA;
#endif  
        
        return hits;
}


void
RayCaster::CastRays(
                    SimpleRayContainer &rays,
                    VssRayContainer &vssRays,
                    const AxisAlignedBox3 &sbox,
                    const bool castDoubleRay,
                    const bool pruneInvalidRays )
{
  SimpleRayContainer::const_iterator rit, rit_end = rays.end();
  
  for (rit = rays.begin(); rit != rit_end; ++ rit) {
    
    CastRay(
            *rit,                               
            vssRays,
            sbox,
            castDoubleRay,
            pruneInvalidRays);
        }
}
  

}