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

Revision 2686, 23.4 KB checked in by mattausch, 16 years ago (diff)

fixed several problems

RevLine 
[1991]1#include "SamplingStrategy.h"
2#include "Ray.h"
3#include "Intersectable.h"
4#include "Preprocessor.h"
5#include "ViewCellsManager.h"
6#include "AxisAlignedBox3.h"
7#include "RssTree.h"
8#include "Vector2.h"
9#include "RndGauss.h"
10#include "Mutation.h"
[2001]11#include "Exporter.h"
[2060]12#include "Environment.h"
[2076]13#include "RayCaster.h"
[1991]14
15#ifdef GTP_INTERNAL
16#include "ArchModeler2MLRT.hxx"
17#endif
18
19namespace GtpVisibilityPreprocessor {
20
[2599]21
22
[1991]23#define MUTATION_USE_CDF 0
[1997]24
[2060]25  //#define SIL_TERMINATION_MUTATION_PROB 0.8f
[1997]26
[1991]27#define EVALUATE_MUTATION_STATS 1
28
[2590]29//#define Q_SEARCH_STEPS 3
[2001]30
[2590]31// VH - commened out !!!!! 8/1/2008
[2583]32#define SORT_RAY_ENTRIES 1
[2011]33
34// use avg ray contribution as importance
35// if 0 the importance is evaluated from the succ of mutations
[2035]36#define USE_AVG_CONTRIBUTION 1
[2011]37
38MutationBasedDistribution::RayEntry &
39MutationBasedDistribution::GetEntry(const int index)
40{
41#if SORT_RAY_ENTRIES
42  return mRays[index];
43#else
44  return mRays[(mBufferStart+index)%mRays.size()];
45#endif
46}
47
[1991]48void
49MutationBasedDistribution::Update(VssRayContainer &vssRays)
50{
51  //  for (int i=0; i < mRays.size(); i++)
52  //    cout<<mRays[i].mMutations<<" ";
53  //  cout<<endl;
[2353]54  cerr<<"Mutation update..."<<endl;
[1991]55  cerr<<"rays = "<<mRays.size()<<endl;
[2123]56  if ((int)mRays.size()) {
[1991]57        cerr<<"Oversampling factors = "<<
58          GetEntry(0).mMutations<<" "<<
59          GetEntry(1).mMutations<<" "<<
60          GetEntry(2).mMutations<<" "<<
61          GetEntry(3).mMutations<<" "<<
62          GetEntry(4).mMutations<<" "<<
63          GetEntry(5).mMutations<<" ... "<<
[2123]64          GetEntry((int)mRays.size()-6).mMutations<<" "<<
65          GetEntry((int)mRays.size()-5).mMutations<<" "<<
66          GetEntry((int)mRays.size()-4).mMutations<<" "<<
67          GetEntry((int)mRays.size()-3).mMutations<<" "<<
68          GetEntry((int)mRays.size()-2).mMutations<<" "<<
69          GetEntry((int)mRays.size()-1).mMutations<<endl;
[1991]70  }
71  int contributingRays = 0;
72
73  int mutationRays = 0;
74  int dummyNcMutations = 0;
75  int dummyCMutations = 0;
[2001]76
77  int reverseCandidates = 0;
[2011]78
79#if 0
80  sort(mRays.begin(), mRays.end());
81  // reset the start of the buffer
82  mBufferStart = 0;
83#endif
84
[2606]85  for (int i=0; i < vssRays.size(); i++)
86        if (vssRays[i]->mTerminationObject) {
87          if (vssRays[i]->mPvsContribution) {
88                // reset the counter of unsuccsseful mutation for a generating ray (if it exists)
89                if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
90                        vssRays[i]->mGeneratorId != -1
[1991]91                  ) {
[2606]92                  mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
[1991]93#if EVALUATE_MUTATION_STATS
[2606]94                  mutationRays++;
95                 
96                  Intersectable *newObject = vssRays[i]->mTerminationObject;
[1991]97               
[2606]98                  Intersectable *oldObject = mRays[vssRays[i]->mGeneratorId].mRay->mTerminationObject;
99                  // the ray generated a contribution although it hits the same object
[2592]100                // mark this using a counter
[1991]101                if (oldObject == newObject)
102                  dummyCMutations++;
103#endif
104          }
105          contributingRays++;
106          if (mRays.size() < mMaxRays) {
107                VssRay *newRay = new VssRay(*vssRays[i]);
108                // add this ray
109                newRay->Ref();
110                mRays.push_back(RayEntry(newRay));
111          } else {
[2592]112                // unref the old ray and add the new ray into the mutation buffer
[1991]113                *mRays[mBufferStart].mRay = *vssRays[i];
114                mRays[mBufferStart].mMutations = 0;
[2011]115                mRays[mBufferStart].mUnsuccessfulMutations = 0;
[2001]116                mRays[mBufferStart].ResetReverseMutation();
[1991]117                //              mRays[mBufferStart] = RayEntry(newRay);
118                mBufferStart++;
119                if (mBufferStart >= mMaxRays)
120                  mBufferStart = 0;
121          }
122        } else {
[2592]123          // the ray did not generate any contribution
[1991]124          if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
125                  vssRays[i]->mGeneratorId != -1
126                  ) {
[2001]127                // check whether not to store a new backward mutation candidate
[2592]128                // this happens if the new ray is occluded significantly closer than
129                // the generator ray
[2001]130                VssRay *oldRay = mRays[vssRays[i]->mGeneratorId].mRay;
131                VssRay *newRay = vssRays[i];
[1991]132
[2021]133                Intersectable *oldObject = oldRay->mTerminationObject;
[1991]134
[2599]135                 
[2592]136                // only allow one reverse mutation per generator ray
[2001]137                if (!mRays[newRay->mGeneratorId].HasReverseMutation()) {
138                  if (DotProd(oldRay->GetDir(), newRay->GetDir()) > 0.0f) {
139                  float oldDist = Magnitude(oldRay->mTermination - newRay->mOrigin);
140                  float newDist = Magnitude(newRay->mTermination - newRay->mOrigin);
141                 
[2060]142                  if (newDist < oldDist - oldObject->GetBox().Radius()*mReverseSamplesDistance) {
[2001]143                        Vector3 origin, termination;
144                        if (ComputeReverseMutation(*oldRay, *newRay, origin, termination)) {
145                          mRays[newRay->mGeneratorId].SetReverseMutation(origin, termination);
146                        }
147                        reverseCandidates++;
148                        //mReverseCandidates
149                  }
150                  }
151                }
152#if EVALUATE_MUTATION_STATS
153                mutationRays++;
154               
[2021]155                Intersectable *newObject = vssRays[i]->mTerminationObject;
[1991]156
157                if (oldObject == newObject)
158                  dummyNcMutations++;
[2001]159#endif
[1991]160          }
161        }
162  }
[2001]163 
[1991]164  if (mutationRays) {
165        cout<<"Mutated rays:"<<mutationRays<<endl;
166        cout<<"Dummy mutations ratio:"<<100.0f*(dummyCMutations + dummyNcMutations)/
167          (float)mutationRays<<"%"<<endl;
168        cout<<"Dummy NC mutations ratio:"<<100.0f*dummyNcMutations/(float)mutationRays<<"%"<<endl;
169        cout<<"Dummy C mutations ratio:"<<100.0f*dummyCMutations/(float)mutationRays<<"%"<<endl;
[2592]170        cout<<"Reverse candidates:"<<100.0f*reverseCandidates/(float)mutationRays<<endl;
[1991]171  }
[2011]172 
[1991]173  float pContributingRays = contributingRays/(float)vssRays.size();
[2011]174 
[2592]175  cout<<"Ratio of contributing rays:"<<pContributingRays<<endl;
[2011]176 
[2060]177  if (mUseUnsuccCountImportance) {
178        // use unsucc mutation samples as feedback on importance
179        for (int i=0; i < mRays.size(); i++) {
180          const float  minImportance = 0.1f;
181          const int minImportanceSamples = 20;
182          mRays[i].mImportance = minImportance +
183                (1-minImportance)*exp(-3.0f*mRays[i].mUnsuccessfulMutations/minImportanceSamples);
184        }
185  } else {
186        float importance = 1.0f;
187        if (mUsePassImportance)
[2686]188          importance = 1.0f/(pContributingRays + 1e-5f);
[2060]189        // set this values for last contributingRays
190        int index = mBufferStart - 1;
[2011]191       
[2060]192        for (int i=0; i < contributingRays; i++, index--) {
193          if (index < 0)
[2124]194                index = (int)mRays.size()-1;
[2060]195          mRays[index].mImportance = importance;
196        }
[2011]197  }
198 
199#if SORT_RAY_ENTRIES
200  long t1 = GetTime();
201  sort(mRays.begin(), mRays.end());
202  // reset the start of the buffer
203  mBufferStart = 0;
[2124]204  mLastIndex = (int)mRays.size();
[2011]205  cout<<"Mutation candidates sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
206#endif
[2590]207                                                                                                                                                                 
[1991]208#if MUTATION_USE_CDF
209  // compute cdf
210  mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mMutations+1);
211  for (int i=1; i < mRays.size(); i++)
212        mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mMutations+1);
213 
214  float scale = 1.0f/mRays[i-1].mCdf;
215  for (i=0; i < mRays.size(); i++) {
216        mRays[i].mCdf *= scale;
217  }
218#endif
219 
220  cout<<"Importance = "<<
221        GetEntry(0).mImportance<<" "<<
[2124]222        GetEntry((int)mRays.size()-1).mImportance<<endl;
[2011]223
224  cout<<"Sampling factor = "<<
225        GetEntry(0).GetSamplingFactor()<<" "<<
[2124]226        GetEntry((int)mRays.size()-1).GetSamplingFactor()<<endl;
[2011]227
[1991]228  cerr<<"Mutation update done."<<endl;
229}
230
231
232Vector3
233MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray,
234                                                                                                 const Vector3 &U,
235                                                                                                 const Vector3 &V,
236                                                                                                 const Vector2 vr2,
237                                                                                                 const float radius
238                                                                                                 )
239{
240#if 0
241  return (U*(r[0] - 0.5f) + V*(r[1] - 0.5f))*(2*radius);
242#endif
243
244
245  // Output random variable
246  Vector2 gaussvec2;
247 
248  // Here we apply transform to gaussian, so 2D bivariate
249  // normal distribution
250  //  float sigma = ComputeSigmaFromRadius(radius);
251  float sigma = radius;
252  GaussianOn2D(vr2,
253                           sigma, // input
254                           gaussvec2);  // output
255
256       
257    // Here we tranform the point correctly to 3D space using base
258    // vectors of the 3D space defined by the direction
259    Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
260       
261        //      cout<<shift<<endl;
262        return shift;
263}
264
265Vector3
266MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray,
267                                                                                                          const Vector3 &U,
268                                                                                                          const Vector3 &V,
269                                                                                                          const Vector2 vr2,
270                                                                                                          const float radius
271                                                                                                          )
272{
273#if 0
274  return (U*(vr2.xx - 0.5f) + V*(vr2.yy - 0.5f))*(4.0f*radius);
275#endif
276  Vector2 gaussvec2;
277#if 1
278  float sigma = radius;
279  GaussianOn2D(vr2,
280                           sigma, // input
281                           gaussvec2);  // output
282  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
283  //    cout<<shift<<endl;
284  return shift;
285#endif
286#if 0
287  // Here we estimate standard deviation (sigma) from radius
288  float sigma = 1.1f*ComputeSigmaFromRadius(radius);
289  Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1));
290  PolarGaussianOnDisk(vr3,
291                                          sigma,
292                                          radius, // input
293                                          gaussvec2); // output
294 
295  // Here we tranform the point correctly to 3D space using base
296  // vectors of the 3D space defined by the direction
297  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
298 
299  //    cout<<shift<<endl;
300  return shift;
301#endif
302}
303
[2001]304bool
305MutationBasedDistribution::ComputeReverseMutation(
306                                                                                                  const VssRay &oldRay,
307                                                                                                  const VssRay &newRay,
308                                                                                                  Vector3 &origin,
309                                                                                                  Vector3 &termination
310                                                                                                  )
311{
[2063]312#ifdef GTP_INTERNAL
[2606]313  static int counter = 0;
[2063]314
[2606]315 
[2001]316  // first reconstruct the termination point
317  Vector3 oldDir = Normalize(oldRay.GetDir());
318  Plane3 oldPlane(oldDir, oldRay.mTermination);
319 
320  termination = oldPlane.FindIntersection(newRay.mOrigin,
321                                                                                  newRay.mTermination);
322 
323  // now find the new origin of the ray by casting ray backward from the termination and termining
324  // silhouette point with respect to the occluding object (object containing the newRay termination)
[1991]325
[2001]326  Plane3 newPlane(oldDir, newRay.mTermination);
327
328  Vector3 oldPivot = newPlane.FindIntersection(oldRay.mOrigin,
329                                                                                           oldRay.mTermination);
330
331  Vector3 newPivot = newRay.mTermination;
332  Vector3 line = 2.0f*(oldPivot - newPivot);
333
[2021]334  Intersectable *occluder = newRay.mTerminationObject;
[2606]335
336
[2001]337  AxisAlignedBox3 box = occluder->GetBox();
[2592]338  // consider slightly larger neighborhood of the occluder
339  // in search for unocclude reverse ray
[2001]340  box.Scale(2.0f);
[2599]341  //box.Scale(200.0f);
[2606]342
[2001]343 
344  const int packetSize = 4;
345  static int hit_triangles[packetSize];
346  static float dist[packetSize];
347  static Vector3 dirs[packetSize];
348  static Vector3 shifts[packetSize];
[2582]349  static Vector3 origs[packetSize];
[2001]350  // now find the silhouette along the line
351  int i;
352  float left = 0.0f;
353  float right = 1.0f;
[2063]354
[2001]355  // cast rays to find silhouette ray
[2060]356  for (int j=0; j < mSilhouetteSearchSteps; j++) {
[2001]357        for (i=0; i < packetSize; i++) {
358          float r = left + (i+1)*(right-left)/(packetSize+1);
359          shifts[i] = r*line;
360          dirs[i] = Normalize(newPivot + shifts[i] - termination );
[2582]361          origs[i] = termination;
362          // mlrtaStoreRayASEye4(&termination.x, &dirs[i].x, i);
[2001]363        }
[2582]364
365        // mlrtaTraverseGroupASEye4(&box.Min().x, &box.Max().x, hit_triangles, dist);
366        assert(preprocessor->mRayCaster);
367        preprocessor->mRayCaster->CastRaysPacket4(box.Min(),
368                                                                                          box.Max(),
[2583]369                                                                                          origs,
370                                                                                          dirs,
[2582]371                                                                                          hit_triangles,
372                                                                                          dist);       
[2001]373       
374        for (i=0; i < packetSize; i++) {
[2583]375          //cout<<hit_triangles[i]<<endl;
[2001]376          if (hit_triangles[i] == -1) {
[2583]377                // break on first passing ray
[2582]378            break;
[2599]379          }
[2001]380        }
381        float rr = left + (i+1)*(right-left)/(packetSize+1);
382        float rl = left + i*(right-left)/(packetSize+1);
383        left = rl;
384        right = rr;
385  }
386
[2599]387
[2001]388  float t = right;
389  if (right==1.0f)
390        return false;
[2599]391
392
[2001]393  if (i == packetSize)
394        origin = newPivot + right*line;
395  else
396        origin = newPivot + shifts[i];
397
398  if (0) {
399       
400        static VssRayContainer rRays;
401        static int counter = 0;
402        char filename[256];
403
404        if (counter < 50) {
405          sprintf(filename, "reverse_rays_%03d.x3d", counter++);
406         
407          VssRay tRay(origin, termination, NULL, NULL);
408          rRays.push_back((VssRay *)&oldRay);
409          rRays.push_back((VssRay *)&newRay);
410          rRays.push_back(&tRay);
411         
412          Exporter *exporter = NULL;
413          exporter = Exporter::GetExporter(filename);
414         
415          exporter->SetFilled();
416         
[2021]417          Intersectable *occludee =
418                oldRay.mTerminationObject;
[2001]419         
420          exporter->SetForcedMaterial(RgbColor(0,0,1));
421          exporter->ExportIntersectable(occluder);
422          exporter->SetForcedMaterial(RgbColor(0,1,0));
423          exporter->ExportIntersectable(occludee);
424          exporter->ResetForcedMaterial();
425         
426          exporter->SetWireframe();
427         
428         
429          exporter->ExportRays(rRays, RgbColor(1, 0, 0));
430          delete exporter;
431          rRays.clear();
432        }
433  }
434
[2592]435  return true;
[2063]436#else
[2574]437  cerr << "warning: reverse mutation not supported!" << endl;
[2592]438  return false;
[2063]439#endif
[2001]440 
441 
442  // now the origin and termination is swapped compred to the generator ray
443  // swap(origin, termination);???
444  // -> perhaps not neccessary as the reverse mutation wil only be used once!
445}
446
[1997]447Vector3
448MutationBasedDistribution::ComputeSilhouetteTerminationMutation(const VssRay &ray,
449                                                                                                                                const Vector3 &origin,
450                                                                                                                                const AxisAlignedBox3 &box,
451                                                                                                                                const Vector3 &U,
452                                                                                                                                const Vector3 &V,
453                                                                                                                                const float radius
454                                                                                                                                )
455{
[2063]456#ifdef GTP_INTERNAL
457
[2001]458  const int packetSize = 4;
[1997]459  static int hit_triangles[packetSize];
460  static float dist[packetSize];
461  static Vector3 dirs[packetSize];
[2001]462  static Vector3 shifts[packetSize];
[2582]463  static Vector3 origs[packetSize];
[1997]464  // mutate the
[2071]465  float alpha = RandomValue(0.0f, Real(2.0*M_PI));
[2001]466  //float alpha = vr2.x*2.0f*M_PI;
[1997]467 
468  // direction along which we will mutate the ray
469  Vector3 line = sin(alpha)*U + cos(alpha)*V;
470 
471  //  cout<<line<<endl;
472  // create 16 rays along the selected dir
[2001]473  int i;
474  float left = 0.0f;
475  float right = radius;
[2599]476 
477  AxisAlignedBox3 _box = box;
478  _box.Scale(1.1f);
[1997]479  // cast rays to find silhouette ray
[2060]480  for (int j=0; j < mSilhouetteSearchSteps; j++) {
[2001]481        for (i=0; i < packetSize; i++) {
482          float r = left + (i+1)*(right-left)/(packetSize+1);
483          shifts[i] = r*line;
484          dirs[i] = Normalize(ray.mTermination + shifts[i] - origin );
[2582]485          origs[i] = origin;
486          //mlrtaStoreRayASEye4(&origin.x, &dirs[i].x, i);
[1997]487        }
[2001]488       
[2582]489        // mlrtaTraverseGroupASEye4(&box.Min().x, &box.Max().x, hit_triangles, dist);
490        assert(preprocessor->mRayCaster);
[2599]491        preprocessor->mRayCaster->CastRaysPacket4(_box.Min(),
492                                                                                          _box.Max(),
493                                                                                          origs,
494                                                                                          dirs,
495                                                                                          hit_triangles,
496                                                                                          dist);       
[2001]497       
498        for (i=0; i < packetSize; i++) {
499          if (hit_triangles[i] == -1) {
500                //        if (hit_triangles[i] == -1 || !box.IsInside(origin + dist[i]*dirs[i])) {
501                // break on first passing ray
502                break;
503          }
504        }
505        float rr = left + (i+1)*(right-left)/(packetSize+1);
506        float rl = left + i*(right-left)/(packetSize+1);
507        left = rl;
508        right = rr;
[1997]509  }
[2599]510
511  Vector3 shift;
[1997]512 
513  if (i == packetSize) {
514        //      cerr<<"Warning: hit the same box here should never happen!"<<endl;
515        // shift the ray even a bit more
[2001]516        //cout<<"W"<<i<<endl;
517        //      return (RandomValue(1.0f, 1.5f)*radius)*line;
[2599]518        shift = right*line;
519  } else {
520        //  cout<<i<<endl;
521        shift = shifts[i];
[1997]522  }
[2001]523 
[2063]524
[2599]525  if (0) {
526       
527        static VssRayContainer rRays;
528        static int counter = 0;
529        char filename[256];
530
531        if (counter < 50) {
532          sprintf(filename, "sil_rays_%03d.x3d", counter++);
533         
534          VssRay tRays[10];
535          rRays.push_back((VssRay *)&ray);
536          for (int k=0; k < packetSize; k++)
537                if (k!=i) {
538                  tRays[k] = VssRay(origin, ray.mTermination + shifts[k], NULL, NULL);
539                  rRays.push_back(&tRays[k]);
540                }
541         
542          Exporter *exporter = NULL;
543          exporter = Exporter::GetExporter(filename);
544         
545          exporter->SetFilled();
546         
547          Intersectable *occluder =
548                ray.mTerminationObject;
549
550          // cout<<occluder->Type()<<endl;
551         
552          exporter->SetForcedMaterial(RgbColor(0,0,1));
553          exporter->ExportIntersectable(occluder);
554
555          exporter->SetWireframe();
556         
557          exporter->SetForcedMaterial(RgbColor(0,1,0));
558          exporter->ExportBox(occluder->GetBox());
559
560          exporter->SetForcedMaterial(RgbColor(0,1,1));
561          exporter->ExportBox(_box);
562
563          exporter->ResetForcedMaterial();
564
565          exporter->ExportRays(rRays, RgbColor(1, 0, 0));
566
567          rRays.clear();
568          tRays[0] = VssRay(origin, ray.mTermination+shift, NULL, NULL);
569          rRays.push_back((VssRay *)&tRays[0]);
570          exporter->ExportRays(rRays, RgbColor(1, 1, 0));
571
572          delete exporter;
573          rRays.clear();
574        }
575  }
576
577  return shift;
578
579 
[2063]580#else
[2592]581  //cerr << "warning: shiluette mutation not supported" << endl;
582  return Vector3(0, 0, 0);
[2063]583#endif
584 
[1997]585}
586
587
[1991]588bool
589MutationBasedDistribution::GenerateSample(SimpleRay &sray)
590{
591
592  if (mRays.size() == 0) {
593        float rr[5];
[2592]594        // use direction based distribution until we have some mutation candidates
[1991]595        Vector3 origin, direction;
[2076]596
597        for (int i=0; i < 5; i++)
598          rr[i] = RandomValue(0,1);
[1991]599       
600        mPreprocessor.mViewCellsManager->GetViewPoint(origin,
601                                                                                                  Vector3(rr[0], rr[1], rr[2]));
602       
[2592]603       
[1991]604        direction = UniformRandomVector(rr[3], rr[4]);
605       
606        const float pdf = 1.0f;
607        sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
608        sray.mGeneratorId = -1;
[2592]609       
[1991]610        return true;
[2592]611  }
[1991]612
613  int index;
614
615#if !MUTATION_USE_CDF
[2011]616#if SORT_RAY_ENTRIES
[2599]617  // RAYS are sorted -> find mitation candidate from the tail of the buffer
[2011]618  index = mLastIndex - 1;
619  if (index < 0 || index >= mRays.size()-1) {
[2124]620        index = (int)mRays.size() - 1;
[2011]621  } else
622        if (
623                mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
[2592]624          // make another round to equalize the oversampling factor
[2011]625
626          //      cout<<"R2"<<endl;
627          //      cout<<mLastIndex<<endl;
628          //      cout<<index<<endl;
[2124]629          index = (int)mRays.size() - 1;
[2011]630        }
631#else
[1991]632  // get tail of the buffer
633  index = (mLastIndex+1)%mRays.size();
634  if (mRays[index].GetSamplingFactor() >
635          mRays[mLastIndex].GetSamplingFactor()) {
636        // search back for index where this is valid
637        index = (mLastIndex - 1 + mRays.size())%mRays.size();
638        for (int i=0; i < mRays.size(); i++) {
639         
640          //      if (mRays[index].mMutations > mRays[mLastIndex].mMutations)
641          //            break;
642          if (mRays[index].GetSamplingFactor() >
643                  mRays[mLastIndex].GetSamplingFactor() )
644                break;
645          index = (index - 1 + mRays.size())%mRays.size();
646        }
647        // go one step back
648        index = (index+1)%mRays.size();
649  }
[2011]650#endif
[1991]651#else
[2599]652
[1991]653  static HaltonSequence iHalton;
654  iHalton.GetNext(1, rr);
655  //rr[0] = RandomValue(0,1);
656  // use binary search to find index with this cdf
[2210]657  int l=0, r=(int)mRays.size()-1;
[1991]658  while(l<r) {
659        int i = (l+r)/2;
660        if (rr[0] < mRays[i].mCdf )
661          r = i;
662        else
663          l = i+1;
664  }
665  index = l;
666  //  if (rr[0] >= mRays[r].mCdf)
667  //    index = r;
668  //  else
669  //    index = l;
670       
671       
672#endif
673  //  cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl;
674
[2599]675  // WE HAVE THE INDEX HERE
676
[1991]677  mLastIndex = index;
[2011]678  //  Debug<<index<<" "<<mRays[index].GetSamplingFactor()<<endl;
679 
[2001]680  if (mRays[index].HasReverseMutation()) {
681        //cout<<"R "<<mRays[index].mutatedOrigin<<" "<<mRays[index].mutatedTermination<<endl;
682        sray = SimpleRay(mRays[index].mutatedOrigin,
683                                         Normalize(mRays[index].mutatedTermination - mRays[index].mutatedOrigin),
[2105]684                                         MUTATION_BASED_DISTRIBUTION,
685                                         1.0f);
686       
687       
[2001]688        sray.mGeneratorId = index;
689        mRays[index].ResetReverseMutation();
690        mRays[index].mMutations++;
[2002]691        mRays[index].mUnsuccessfulMutations++;
[2574]692       
[2001]693        return true;
694  }
[2574]695
[1991]696  return GenerateMutation(index, sray);
697}
698
699
700
701
702
703bool
704MutationBasedDistribution::GenerateMutationCandidate(const int index,
705                                                                                                         SimpleRay &sray,
706                                                                                                         Intersectable *object,
707                                                                                                         const AxisAlignedBox3 &box
708                                                                                                         )
709{
710  float rr[4];
711
712  VssRay *ray = mRays[index].mRay;
713
[2001]714  //  rr[0] = RandomValue(0.0f,0.99999f);
715  //  rr[1] = RandomValue(0.0f,0.99999f);
716  //  rr[2] = RandomValue(0.0f,0.99999f);
717  //  rr[3] = RandomValue(0.0f,0.99999f);
[1991]718 
719  // mutate the origin
720  Vector3 d = ray->GetDir();
721 
[2022]722  float objectRadius = box.Radius();
[1991]723  //  cout<<objectRadius<<endl;
724  if (objectRadius < Limits::Small)
725        return false;
726 
727  // Compute right handed coordinate system from direction
728  Vector3 U, V;
729  Vector3 nd = Normalize(d);
730  nd.RightHandedBase(U, V);
731 
732  Vector3 origin = ray->mOrigin;
733  Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center();
[1995]734
[1991]735 
[2022]736  // use probabilitistic approach to decide for the type of mutation
737  float a = RandomValue(0.0f,1.0f);
[2105]738  bool bidirectional = true;
739 
[2060]740  if (mUseSilhouetteSamples && a < mSilhouetteProb) {
[2022]741        termination += ComputeSilhouetteTerminationMutation(*ray,
742                                                                                                                origin,
743                                                                                                                box,
744                                                                                                                U, V,
745                                                                                                                2.0f*objectRadius);
[2105]746        bidirectional = false;
[2022]747  } else {
748        mRays[index].mHalton.GetNext(4, rr);
[2076]749
[2022]750        // fuzzy random mutation
751        origin += ComputeOriginMutation(*ray, U, V,
752                                                                        Vector2(rr[0], rr[1]),
[2060]753                                                                        mMutationRadiusOrigin*objectRadius);
[2022]754       
755        termination += ComputeTerminationMutation(*ray, U, V,
756                                                                                          Vector2(rr[2], rr[3]),
[2060]757                                                                                          mMutationRadiusTermination*objectRadius);
[2022]758  }
759
[1991]760  Vector3 direction = termination - origin;
761 
762  if (Magnitude(direction) < Limits::Small)
763        return false;
764 
765  // shift the origin a little bit
[2592]766#if 1
[1991]767  origin += direction*0.5f;
[2592]768#else
769  // $$JB - 14.1. 2008 increase shift to test HavranRayCaster
770  origin = 0.5f*(origin + termination);
771#endif
[1991]772 
773  direction.Normalize();
774 
775  // $$ jb the pdf is yet not correct for all sampling methods!
776  const float pdf = 1.0f;
777 
778  sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
779  sray.mGeneratorId = index;
[2105]780  sray.SetBidirectional(bidirectional);
781 
[2048]782  return true;
[1991]783}
784
785bool
786MutationBasedDistribution::GenerateMutation(const int index, SimpleRay &sray)
787{
788  VssRay *ray = mRays[index].mRay;
789
[2021]790  Intersectable *object = ray->mTerminationObject;
[1991]791 
792  AxisAlignedBox3 box = object->GetBox();
793 
794  if (GenerateMutationCandidate(index, sray, object, box)) {
795    mRays[index].mMutations++;
[2002]796        mRays[index].mUnsuccessfulMutations++;
[1991]797        return true;
798  }
799  return false;
800}
[2060]801 
[1991]802
[2060]803
804MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor
805                                                                                                         ) :
806  SamplingStrategy(preprocessor)
[1991]807{
[2060]808  mType = MUTATION_BASED_DISTRIBUTION;
809  mBufferStart = 0;
810  mLastIndex = 0;
811
812  Environment::GetSingleton()->GetIntValue("Mutation.bufferSize",
813                                                                                   mMaxRays);
[1991]814 
[2060]815  mRays.reserve(mMaxRays);
[1991]816 
[2060]817  Environment::GetSingleton()->GetFloatValue("Mutation.radiusOrigin",
818                                                                                         mMutationRadiusOrigin);
[1991]819
[2060]820  Environment::GetSingleton()->GetFloatValue("Mutation.radiusTermination",
821                                                                                         mMutationRadiusTermination);
[1991]822 
[2076]823  bool useEnhancedFeatures;
824
825#ifdef GTP_INTERNAL
826  int rayCaster;
[2172]827  Environment::GetSingleton()->GetIntValue("Preprocessor.rayCastMethod",
[2076]828                                                                                   rayCaster);
[2583]829  useEnhancedFeatures = (rayCaster !=
830                                                 RayCaster::INTERNAL_RAYCASTER);
[2076]831#else
832  useEnhancedFeatures = false;
833#endif
[1991]834 
[2076]835  if (useEnhancedFeatures) {
836        Environment::GetSingleton()->GetBoolValue("Mutation.useReverseSamples",
837                                                                                          mUseReverseSamples);
838       
839        Environment::GetSingleton()->GetFloatValue("Mutation.reverseSamplesDistance",
840
841                                                                                           mReverseSamplesDistance);
842
843        Environment::GetSingleton()->GetFloatValue("Mutation.silhouetteProb",
844                                                                                           mSilhouetteProb);
845
846  } else {
847        mUseReverseSamples = false;
848        mReverseSamplesDistance = 1e20f;
849        mSilhouetteProb = 0.0f;
850  }
[2060]851 
852  Environment::GetSingleton()->GetBoolValue("Mutation.useSilhouetteSamples",
853                                                                                        mUseSilhouetteSamples);
[1991]854
[2060]855  Environment::GetSingleton()->GetIntValue("Mutation.silhouetteSearchSteps",
856                                                                                   mSilhouetteSearchSteps);
[1991]857 
[2076]858 
[2060]859  Environment::GetSingleton()->GetBoolValue("Mutation.usePassImportance",
860                                                                                        mUsePassImportance);
[1991]861 
862
[2060]863  Environment::GetSingleton()->GetBoolValue("Mutation.useUnsuccCountImportance",
864                                                                                        mUseUnsuccCountImportance);
[1991]865
866 
[2060]867 
[1991]868}
869
870
[2011]871
[1991]872}
Note: See TracBrowser for help on using the repository browser.