source: GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp @ 2686

#include "SamplingStrategy.h"
#include "Ray.h"
#include "Intersectable.h"
#include "Preprocessor.h"
#include "ViewCellsManager.h"
#include "AxisAlignedBox3.h"
#include "RssTree.h"
#include "Vector2.h"
#include "RndGauss.h"
#include "Mutation.h"
#include "Exporter.h"
#include "Environment.h"
#include "RayCaster.h"

#ifdef GTP_INTERNAL
#include "ArchModeler2MLRT.hxx"
#endif

namespace GtpVisibilityPreprocessor {



#define MUTATION_USE_CDF 0

  //#define SIL_TERMINATION_MUTATION_PROB 0.8f

#define EVALUATE_MUTATION_STATS 1

//#define Q_SEARCH_STEPS 3

// VH - commened out !!!!! 8/1/2008
#define SORT_RAY_ENTRIES 1

// use avg ray contribution as importance
// if 0 the importance is evaluated from the succ of mutations
#define USE_AVG_CONTRIBUTION 1

MutationBasedDistribution::RayEntry &
MutationBasedDistribution::GetEntry(const int index)
{
#if SORT_RAY_ENTRIES
  return mRays[index];
#else
  return mRays[(mBufferStart+index)%mRays.size()];
#endif
}

void
MutationBasedDistribution::Update(VssRayContainer &vssRays)
{
  //  for (int i=0; i < mRays.size(); i++)
  //    cout<<mRays[i].mMutations<<" ";
  //  cout<<endl;
  cerr<<"Mutation update..."<<endl;
  cerr<<"rays = "<<mRays.size()<<endl;
  if ((int)mRays.size()) {
        cerr<<"Oversampling factors = "<<
          GetEntry(0).mMutations<<" "<<
          GetEntry(1).mMutations<<" "<<
          GetEntry(2).mMutations<<" "<<
          GetEntry(3).mMutations<<" "<<
          GetEntry(4).mMutations<<" "<<
          GetEntry(5).mMutations<<" ... "<<
          GetEntry((int)mRays.size()-6).mMutations<<" "<<
          GetEntry((int)mRays.size()-5).mMutations<<" "<<
          GetEntry((int)mRays.size()-4).mMutations<<" "<<
          GetEntry((int)mRays.size()-3).mMutations<<" "<<
          GetEntry((int)mRays.size()-2).mMutations<<" "<<
          GetEntry((int)mRays.size()-1).mMutations<<endl;
  }
  int contributingRays = 0;

  int mutationRays = 0;
  int dummyNcMutations = 0;
  int dummyCMutations = 0;

  int reverseCandidates = 0;

#if 0
  sort(mRays.begin(), mRays.end());
  // reset the start of the buffer
  mBufferStart = 0;
#endif

  for (int i=0; i < vssRays.size(); i++) 
        if (vssRays[i]->mTerminationObject) {
          if (vssRays[i]->mPvsContribution) {
                // reset the counter of unsuccsseful mutation for a generating ray (if it exists)
                if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
                        vssRays[i]->mGeneratorId != -1
                  ) {
                  mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
#if EVALUATE_MUTATION_STATS 
                  mutationRays++;
                  
                  Intersectable *newObject = vssRays[i]->mTerminationObject;
                
                  Intersectable *oldObject = mRays[vssRays[i]->mGeneratorId].mRay->mTerminationObject;
                  // the ray generated a contribution although it hits the same object
                // mark this using a counter
                if (oldObject == newObject)
                  dummyCMutations++;
#endif
          }
          contributingRays++;
          if (mRays.size() < mMaxRays) {
                VssRay *newRay = new VssRay(*vssRays[i]);
                // add this ray
                newRay->Ref();
                mRays.push_back(RayEntry(newRay));
          } else {
                // unref the old ray and add the new ray into the mutation buffer
                *mRays[mBufferStart].mRay = *vssRays[i];
                mRays[mBufferStart].mMutations = 0;
                mRays[mBufferStart].mUnsuccessfulMutations = 0;
                mRays[mBufferStart].ResetReverseMutation();
                //              mRays[mBufferStart] = RayEntry(newRay);
                mBufferStart++;
                if (mBufferStart >= mMaxRays)
                  mBufferStart = 0;
          }
        } else {
          // the ray did not generate any contribution
          if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
                  vssRays[i]->mGeneratorId != -1
                  ) {
                // check whether not to store a new backward mutation candidate
                // this happens if the new ray is occluded significantly closer than
                // the generator ray
                VssRay *oldRay = mRays[vssRays[i]->mGeneratorId].mRay;
                VssRay *newRay = vssRays[i];

                Intersectable *oldObject = oldRay->mTerminationObject;

                  
                // only allow one reverse mutation per generator ray
                if (!mRays[newRay->mGeneratorId].HasReverseMutation()) {
                  if (DotProd(oldRay->GetDir(), newRay->GetDir()) > 0.0f) {
                  float oldDist = Magnitude(oldRay->mTermination - newRay->mOrigin);
                  float newDist = Magnitude(newRay->mTermination - newRay->mOrigin);
                  
                  if (newDist < oldDist - oldObject->GetBox().Radius()*mReverseSamplesDistance) {
                        Vector3 origin, termination;
                        if (ComputeReverseMutation(*oldRay, *newRay, origin, termination)) {
                          mRays[newRay->mGeneratorId].SetReverseMutation(origin, termination);
                        }
                        reverseCandidates++;
                        //mReverseCandidates
                  }
                  }
                }
#if EVALUATE_MUTATION_STATS 
                mutationRays++;
                
                Intersectable *newObject = vssRays[i]->mTerminationObject;

                if (oldObject == newObject)
                  dummyNcMutations++;
#endif
          }
        }
  }
  
  if (mutationRays) {
        cout<<"Mutated rays:"<<mutationRays<<endl;
        cout<<"Dummy mutations ratio:"<<100.0f*(dummyCMutations + dummyNcMutations)/
          (float)mutationRays<<"%"<<endl;
        cout<<"Dummy NC mutations ratio:"<<100.0f*dummyNcMutations/(float)mutationRays<<"%"<<endl;
        cout<<"Dummy C mutations ratio:"<<100.0f*dummyCMutations/(float)mutationRays<<"%"<<endl;
        cout<<"Reverse candidates:"<<100.0f*reverseCandidates/(float)mutationRays<<endl;
  }
  
  float pContributingRays = contributingRays/(float)vssRays.size();
  
  cout<<"Ratio of contributing rays:"<<pContributingRays<<endl;
  
  if (mUseUnsuccCountImportance) {
        // use unsucc mutation samples as feedback on importance
        for (int i=0; i < mRays.size(); i++) {
          const float  minImportance = 0.1f;
          const int minImportanceSamples = 20;
          mRays[i].mImportance = minImportance +
                (1-minImportance)*exp(-3.0f*mRays[i].mUnsuccessfulMutations/minImportanceSamples);
        }
  } else {
        float importance = 1.0f;
        if (mUsePassImportance) 
          importance = 1.0f/(pContributingRays + 1e-5f);
        // set this values for last contributingRays
        int index = mBufferStart - 1;
        
        for (int i=0; i < contributingRays; i++, index--) {
          if (index < 0)
                index = (int)mRays.size()-1;
          mRays[index].mImportance = importance;
        }
  }
  
#if SORT_RAY_ENTRIES
  long t1 = GetTime();
  sort(mRays.begin(), mRays.end());
  // reset the start of the buffer
  mBufferStart = 0;
  mLastIndex = (int)mRays.size();
  cout<<"Mutation candidates sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
#endif
                                                                                                                                                                  
#if MUTATION_USE_CDF
  // compute cdf
  mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mMutations+1);
  for (int i=1; i < mRays.size(); i++) 
        mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mMutations+1);
  
  float scale = 1.0f/mRays[i-1].mCdf;
  for (i=0; i < mRays.size(); i++) {
        mRays[i].mCdf *= scale;
  }
#endif
  
  cout<<"Importance = "<<
        GetEntry(0).mImportance<<" "<<
        GetEntry((int)mRays.size()-1).mImportance<<endl;

  cout<<"Sampling factor = "<<
        GetEntry(0).GetSamplingFactor()<<" "<<
        GetEntry((int)mRays.size()-1).GetSamplingFactor()<<endl;

  cerr<<"Mutation update done."<<endl;
}


Vector3
MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray,
                                                                                                 const Vector3 &U,
                                                                                                 const Vector3 &V,
                                                                                                 const Vector2 vr2,
                                                                                                 const float radius
                                                                                                 )
{
#if 0
  return (U*(r[0] - 0.5f) + V*(r[1] - 0.5f))*(2*radius);
#endif


  // Output random variable
  Vector2 gaussvec2; 
  
  // Here we apply transform to gaussian, so 2D bivariate
  // normal distribution
  //  float sigma = ComputeSigmaFromRadius(radius);
  float sigma = radius;
  GaussianOn2D(vr2,
                           sigma, // input
                           gaussvec2);  // output

        
    // Here we tranform the point correctly to 3D space using base
    // vectors of the 3D space defined by the direction
    Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
        
        //      cout<<shift<<endl;
        return shift;
}

Vector3
MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray,
                                                                                                          const Vector3 &U,
                                                                                                          const Vector3 &V,
                                                                                                          const Vector2 vr2,
                                                                                                          const float radius
                                                                                                          )
{
#if 0
  return (U*(vr2.xx - 0.5f) + V*(vr2.yy - 0.5f))*(4.0f*radius);
#endif
  Vector2 gaussvec2; 
#if 1
  float sigma = radius;
  GaussianOn2D(vr2,
                           sigma, // input
                           gaussvec2);  // output
  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
  //    cout<<shift<<endl;
  return shift;
#endif
#if 0
  // Here we estimate standard deviation (sigma) from radius
  float sigma = 1.1f*ComputeSigmaFromRadius(radius);
  Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1));
  PolarGaussianOnDisk(vr3,
                                          sigma,
                                          radius, // input
                                          gaussvec2); // output
  
  // Here we tranform the point correctly to 3D space using base
  // vectors of the 3D space defined by the direction
  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
  
  //    cout<<shift<<endl;
  return shift;
#endif
}

bool
MutationBasedDistribution::ComputeReverseMutation(
                                                                                                  const VssRay &oldRay,
                                                                                                  const VssRay &newRay,
                                                                                                  Vector3 &origin,
                                                                                                  Vector3 &termination
                                                                                                  )
{
#ifdef GTP_INTERNAL
  static int counter = 0;

  
  // first reconstruct the termination point
  Vector3 oldDir = Normalize(oldRay.GetDir());
  Plane3 oldPlane(oldDir, oldRay.mTermination);
  
  termination = oldPlane.FindIntersection(newRay.mOrigin,
                                                                                  newRay.mTermination);
  
  // now find the new origin of the ray by casting ray backward from the termination and termining
  // silhouette point with respect to the occluding object (object containing the newRay termination)

  Plane3 newPlane(oldDir, newRay.mTermination);

  Vector3 oldPivot = newPlane.FindIntersection(oldRay.mOrigin,
                                                                                           oldRay.mTermination);

  Vector3 newPivot = newRay.mTermination;
  Vector3 line = 2.0f*(oldPivot - newPivot);

  Intersectable *occluder = newRay.mTerminationObject;


  AxisAlignedBox3 box = occluder->GetBox();
  // consider slightly larger neighborhood of the occluder
  // in search for unocclude reverse ray
  box.Scale(2.0f);
  //box.Scale(200.0f);

  
  const int packetSize = 4;
  static int hit_triangles[packetSize];
  static float dist[packetSize];
  static Vector3 dirs[packetSize];
  static Vector3 shifts[packetSize];
  static Vector3 origs[packetSize];
  // now find the silhouette along the line
  int i;
  float left = 0.0f;
  float right = 1.0f;

  // cast rays to find silhouette ray
  for (int j=0; j < mSilhouetteSearchSteps; j++) {
        for (i=0; i < packetSize; i++) {
          float r = left + (i+1)*(right-left)/(packetSize+1);
          shifts[i] = r*line;
          dirs[i] = Normalize(newPivot + shifts[i] - termination );
          origs[i] = termination;
          // mlrtaStoreRayASEye4(&termination.x, &dirs[i].x, i);
        }

        // mlrtaTraverseGroupASEye4(&box.Min().x, &box.Max().x, hit_triangles, dist);
        assert(preprocessor->mRayCaster);
        preprocessor->mRayCaster->CastRaysPacket4(box.Min(),
                                                                                          box.Max(),
                                                                                          origs,
                                                                                          dirs,
                                                                                          hit_triangles,
                                                                                          dist);       
        
        for (i=0; i < packetSize; i++) {
          //cout<<hit_triangles[i]<<endl;
          if (hit_triangles[i] == -1) {
                // break on first passing ray
            break;
          } 
        }
        float rr = left + (i+1)*(right-left)/(packetSize+1);
        float rl = left + i*(right-left)/(packetSize+1);
        left = rl;
        right = rr;
  }


  float t = right;
  if (right==1.0f)
        return false;


  if (i == packetSize)
        origin = newPivot + right*line;
  else
        origin = newPivot + shifts[i];

  if (0) {
        
        static VssRayContainer rRays;
        static int counter = 0;
        char filename[256];

        if (counter < 50) {
          sprintf(filename, "reverse_rays_%03d.x3d", counter++);
          
          VssRay tRay(origin, termination, NULL, NULL);
          rRays.push_back((VssRay *)&oldRay);
          rRays.push_back((VssRay *)&newRay);
          rRays.push_back(&tRay);
          
          Exporter *exporter = NULL;
          exporter = Exporter::GetExporter(filename);
          
          exporter->SetFilled();
          
          Intersectable *occludee = 
                oldRay.mTerminationObject;
          
          exporter->SetForcedMaterial(RgbColor(0,0,1));
          exporter->ExportIntersectable(occluder);
          exporter->SetForcedMaterial(RgbColor(0,1,0));
          exporter->ExportIntersectable(occludee);
          exporter->ResetForcedMaterial();
          
          exporter->SetWireframe();
          
          
          exporter->ExportRays(rRays, RgbColor(1, 0, 0));
          delete exporter;
          rRays.clear();
        }
  }

  return true;
#else
  cerr << "warning: reverse mutation not supported!" << endl;
  return false;
#endif
  
  
  // now the origin and termination is swapped compred to the generator ray
  // swap(origin, termination);???
  // -> perhaps not neccessary as the reverse mutation wil only be used once!
}

Vector3
MutationBasedDistribution::ComputeSilhouetteTerminationMutation(const VssRay &ray,
                                                                                                                                const Vector3 &origin,
                                                                                                                                const AxisAlignedBox3 &box,
                                                                                                                                const Vector3 &U,
                                                                                                                                const Vector3 &V,
                                                                                                                                const float radius
                                                                                                                                )
{
#ifdef GTP_INTERNAL

  const int packetSize = 4;
  static int hit_triangles[packetSize];
  static float dist[packetSize];
  static Vector3 dirs[packetSize];
  static Vector3 shifts[packetSize];
  static Vector3 origs[packetSize];
  // mutate the 
  float alpha = RandomValue(0.0f, Real(2.0*M_PI));
  //float alpha = vr2.x*2.0f*M_PI;
  
  // direction along which we will mutate the ray
  Vector3 line = sin(alpha)*U + cos(alpha)*V;
  
  //  cout<<line<<endl;
  // create 16 rays along the selected dir
  int i;
  float left = 0.0f;
  float right = radius;
  
  AxisAlignedBox3 _box = box;
  _box.Scale(1.1f);
  // cast rays to find silhouette ray
  for (int j=0; j < mSilhouetteSearchSteps; j++) {
        for (i=0; i < packetSize; i++) {
          float r = left + (i+1)*(right-left)/(packetSize+1);
          shifts[i] = r*line;
          dirs[i] = Normalize(ray.mTermination + shifts[i] - origin );
          origs[i] = origin;
          //mlrtaStoreRayASEye4(&origin.x, &dirs[i].x, i);
        }
        
        // mlrtaTraverseGroupASEye4(&box.Min().x, &box.Max().x, hit_triangles, dist);
        assert(preprocessor->mRayCaster);
        preprocessor->mRayCaster->CastRaysPacket4(_box.Min(),
                                                                                          _box.Max(),
                                                                                          origs,
                                                                                          dirs,
                                                                                          hit_triangles,
                                                                                          dist);       
        
        for (i=0; i < packetSize; i++) {
          if (hit_triangles[i] == -1) {
                //        if (hit_triangles[i] == -1 || !box.IsInside(origin + dist[i]*dirs[i])) {
                // break on first passing ray
                break;
          }
        }
        float rr = left + (i+1)*(right-left)/(packetSize+1);
        float rl = left + i*(right-left)/(packetSize+1);
        left = rl;
        right = rr;
  }

  Vector3 shift;
  
  if (i == packetSize) {
        //      cerr<<"Warning: hit the same box here should never happen!"<<endl;
        // shift the ray even a bit more
        //cout<<"W"<<i<<endl;
        //      return (RandomValue(1.0f, 1.5f)*radius)*line;
        shift = right*line;
  } else {
        //  cout<<i<<endl;
        shift = shifts[i];
  }
  

  if (0) {
        
        static VssRayContainer rRays;
        static int counter = 0;
        char filename[256];

        if (counter < 50) {
          sprintf(filename, "sil_rays_%03d.x3d", counter++);
          
          VssRay tRays[10];
          rRays.push_back((VssRay *)&ray);
          for (int k=0; k < packetSize; k++)
                if (k!=i) {
                  tRays[k] = VssRay(origin, ray.mTermination + shifts[k], NULL, NULL);
                  rRays.push_back(&tRays[k]);
                }
          
          Exporter *exporter = NULL;
          exporter = Exporter::GetExporter(filename);
          
          exporter->SetFilled();
          
          Intersectable *occluder = 
                ray.mTerminationObject;

          // cout<<occluder->Type()<<endl;
          
          exporter->SetForcedMaterial(RgbColor(0,0,1));
          exporter->ExportIntersectable(occluder);

          exporter->SetWireframe();
          
          exporter->SetForcedMaterial(RgbColor(0,1,0));
          exporter->ExportBox(occluder->GetBox());

          exporter->SetForcedMaterial(RgbColor(0,1,1));
          exporter->ExportBox(_box);

          exporter->ResetForcedMaterial();

          exporter->ExportRays(rRays, RgbColor(1, 0, 0));

          rRays.clear();
          tRays[0] = VssRay(origin, ray.mTermination+shift, NULL, NULL);
          rRays.push_back((VssRay *)&tRays[0]);
          exporter->ExportRays(rRays, RgbColor(1, 1, 0));

          delete exporter;
          rRays.clear();
        }
  }

  return shift;

  
#else
  //cerr << "warning: shiluette mutation not supported" << endl;
  return Vector3(0, 0, 0);
#endif
  
}


bool
MutationBasedDistribution::GenerateSample(SimpleRay &sray)
{

  if (mRays.size() == 0) {
        float rr[5];
        // use direction based distribution until we have some mutation candidates
        Vector3 origin, direction;

        for (int i=0; i < 5; i++)
          rr[i] = RandomValue(0,1);
        
        mPreprocessor.mViewCellsManager->GetViewPoint(origin,
                                                                                                  Vector3(rr[0], rr[1], rr[2]));
        
        
        direction = UniformRandomVector(rr[3], rr[4]);
        
        const float pdf = 1.0f;
        sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
        sray.mGeneratorId = -1;
        
        return true;
  }

  int index;

#if !MUTATION_USE_CDF
#if SORT_RAY_ENTRIES
  // RAYS are sorted -> find mitation candidate from the tail of the buffer
  index = mLastIndex - 1;
  if (index < 0 || index >= mRays.size()-1) {
        index = (int)mRays.size() - 1;
  } else
        if (
                mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
          // make another round to equalize the oversampling factor

          //      cout<<"R2"<<endl;
          //      cout<<mLastIndex<<endl;
          //      cout<<index<<endl;
          index = (int)mRays.size() - 1;
        }
#else
  // get tail of the buffer
  index = (mLastIndex+1)%mRays.size();
  if (mRays[index].GetSamplingFactor() >
          mRays[mLastIndex].GetSamplingFactor()) {
        // search back for index where this is valid
        index = (mLastIndex - 1 + mRays.size())%mRays.size();
        for (int i=0; i < mRays.size(); i++) {
          
          //      if (mRays[index].mMutations > mRays[mLastIndex].mMutations)
          //            break;
          if (mRays[index].GetSamplingFactor() >
                  mRays[mLastIndex].GetSamplingFactor() )
                break;
          index = (index - 1 + mRays.size())%mRays.size();
        }
        // go one step back
        index = (index+1)%mRays.size();
  }
#endif
#else

  static HaltonSequence iHalton;
  iHalton.GetNext(1, rr);
  //rr[0] = RandomValue(0,1);
  // use binary search to find index with this cdf
  int l=0, r=(int)mRays.size()-1;
  while(l<r) {
        int i = (l+r)/2;
        if (rr[0] < mRays[i].mCdf )
          r = i;
        else
          l = i+1;
  }
  index = l;
  //  if (rr[0] >= mRays[r].mCdf)
  //    index = r;
  //  else
  //    index = l;
        
        
#endif
  //  cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl;

  // WE HAVE THE INDEX HERE

  mLastIndex = index;
  //  Debug<<index<<" "<<mRays[index].GetSamplingFactor()<<endl;
  
  if (mRays[index].HasReverseMutation()) {
        //cout<<"R "<<mRays[index].mutatedOrigin<<" "<<mRays[index].mutatedTermination<<endl;
        sray = SimpleRay(mRays[index].mutatedOrigin,
                                         Normalize(mRays[index].mutatedTermination - mRays[index].mutatedOrigin),
                                         MUTATION_BASED_DISTRIBUTION,
                                         1.0f);
        
        
        sray.mGeneratorId = index;
        mRays[index].ResetReverseMutation();
        mRays[index].mMutations++;
        mRays[index].mUnsuccessfulMutations++;
        
        return true;
  }

  return GenerateMutation(index, sray);
}





bool
MutationBasedDistribution::GenerateMutationCandidate(const int index,
                                                                                                         SimpleRay &sray,
                                                                                                         Intersectable *object,
                                                                                                         const AxisAlignedBox3 &box
                                                                                                         )
{
  float rr[4];

  VssRay *ray = mRays[index].mRay;

  //  rr[0] = RandomValue(0.0f,0.99999f);
  //  rr[1] = RandomValue(0.0f,0.99999f);
  //  rr[2] = RandomValue(0.0f,0.99999f);
  //  rr[3] = RandomValue(0.0f,0.99999f);
  
  // mutate the origin
  Vector3 d = ray->GetDir();
  
  float objectRadius = box.Radius();
  //  cout<<objectRadius<<endl;
  if (objectRadius < Limits::Small)
        return false;
  
  // Compute right handed coordinate system from direction
  Vector3 U, V;
  Vector3 nd = Normalize(d);
  nd.RightHandedBase(U, V);
  
  Vector3 origin = ray->mOrigin;
  Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center();

  
  // use probabilitistic approach to decide for the type of mutation
  float a = RandomValue(0.0f,1.0f);
  bool bidirectional = true;
  
  if (mUseSilhouetteSamples && a < mSilhouetteProb) {
        termination += ComputeSilhouetteTerminationMutation(*ray,
                                                                                                                origin,
                                                                                                                box,
                                                                                                                U, V,
                                                                                                                2.0f*objectRadius);
        bidirectional = false;
  } else {
        mRays[index].mHalton.GetNext(4, rr);

        // fuzzy random mutation
        origin += ComputeOriginMutation(*ray, U, V,
                                                                        Vector2(rr[0], rr[1]),
                                                                        mMutationRadiusOrigin*objectRadius);
        
        termination += ComputeTerminationMutation(*ray, U, V,
                                                                                          Vector2(rr[2], rr[3]),
                                                                                          mMutationRadiusTermination*objectRadius);
  }

  Vector3 direction = termination - origin;
  
  if (Magnitude(direction) < Limits::Small)
        return false;
  
  // shift the origin a little bit
#if 1
  origin += direction*0.5f;
#else
  // $$JB - 14.1. 2008 increase shift to test HavranRayCaster
  origin = 0.5f*(origin + termination);
#endif
  
  direction.Normalize();
  
  // $$ jb the pdf is yet not correct for all sampling methods!
  const float pdf = 1.0f;
  
  sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
  sray.mGeneratorId = index;
  sray.SetBidirectional(bidirectional);
  
  return true;
}

bool
MutationBasedDistribution::GenerateMutation(const int index, SimpleRay &sray)
{
  VssRay *ray = mRays[index].mRay;

  Intersectable *object = ray->mTerminationObject;
  
  AxisAlignedBox3 box = object->GetBox();
  
  if (GenerateMutationCandidate(index, sray, object, box)) {
    mRays[index].mMutations++;
        mRays[index].mUnsuccessfulMutations++;
        return true;
  }
  return false;
}
  


MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor
                                                                                                         ) :
  SamplingStrategy(preprocessor)
{
  mType = MUTATION_BASED_DISTRIBUTION;
  mBufferStart = 0;
  mLastIndex = 0;

  Environment::GetSingleton()->GetIntValue("Mutation.bufferSize",
                                                                                   mMaxRays);
  
  mRays.reserve(mMaxRays);
  
  Environment::GetSingleton()->GetFloatValue("Mutation.radiusOrigin",
                                                                                         mMutationRadiusOrigin);

  Environment::GetSingleton()->GetFloatValue("Mutation.radiusTermination",
                                                                                         mMutationRadiusTermination);
  
  bool useEnhancedFeatures;

#ifdef GTP_INTERNAL
  int rayCaster;
  Environment::GetSingleton()->GetIntValue("Preprocessor.rayCastMethod",
                                                                                   rayCaster);
  useEnhancedFeatures = (rayCaster !=
                                                 RayCaster::INTERNAL_RAYCASTER);
#else
  useEnhancedFeatures = false;
#endif
  
  if (useEnhancedFeatures) {
        Environment::GetSingleton()->GetBoolValue("Mutation.useReverseSamples",
                                                                                          mUseReverseSamples);
        
        Environment::GetSingleton()->GetFloatValue("Mutation.reverseSamplesDistance",

                                                                                           mReverseSamplesDistance);

        Environment::GetSingleton()->GetFloatValue("Mutation.silhouetteProb",
                                                                                           mSilhouetteProb);

  } else {
        mUseReverseSamples = false;
        mReverseSamplesDistance = 1e20f;
        mSilhouetteProb = 0.0f;
  }
  
  Environment::GetSingleton()->GetBoolValue("Mutation.useSilhouetteSamples",
                                                                                        mUseSilhouetteSamples);

  Environment::GetSingleton()->GetIntValue("Mutation.silhouetteSearchSteps",
                                                                                   mSilhouetteSearchSteps);
  
  
  Environment::GetSingleton()->GetBoolValue("Mutation.usePassImportance",
                                                                                        mUsePassImportance);
  

  Environment::GetSingleton()->GetBoolValue("Mutation.useUnsuccCountImportance",
                                                                                        mUseUnsuccCountImportance);

  
  
}



}