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

Revision 2606, 23.4 KB checked in by bittner, 17 years ago (diff)

merge on nemo

Line 
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"
11#include "Exporter.h"
12#include "Environment.h"
13#include "RayCaster.h"
14
15#ifdef GTP_INTERNAL
16#include "ArchModeler2MLRT.hxx"
17#endif
18
19namespace GtpVisibilityPreprocessor {
20
21
22
23#define MUTATION_USE_CDF 0
24
25  //#define SIL_TERMINATION_MUTATION_PROB 0.8f
26
27#define EVALUATE_MUTATION_STATS 1
28
29//#define Q_SEARCH_STEPS 3
30
31// VH - commened out !!!!! 8/1/2008
32#define SORT_RAY_ENTRIES 1
33
34// use avg ray contribution as importance
35// if 0 the importance is evaluated from the succ of mutations
36#define USE_AVG_CONTRIBUTION 1
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
48void
49MutationBasedDistribution::Update(VssRayContainer &vssRays)
50{
51  //  for (int i=0; i < mRays.size(); i++)
52  //    cout<<mRays[i].mMutations<<" ";
53  //  cout<<endl;
54  cerr<<"Mutation update..."<<endl;
55  cerr<<"rays = "<<mRays.size()<<endl;
56  if ((int)mRays.size()) {
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<<" ... "<<
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;
70  }
71  int contributingRays = 0;
72
73  int mutationRays = 0;
74  int dummyNcMutations = 0;
75  int dummyCMutations = 0;
76
77  int reverseCandidates = 0;
78
79#if 0
80  sort(mRays.begin(), mRays.end());
81  // reset the start of the buffer
82  mBufferStart = 0;
83#endif
84
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
91                  ) {
92                  mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
93#if EVALUATE_MUTATION_STATS
94                  mutationRays++;
95                 
96                  Intersectable *newObject = vssRays[i]->mTerminationObject;
97               
98                  Intersectable *oldObject = mRays[vssRays[i]->mGeneratorId].mRay->mTerminationObject;
99                  // the ray generated a contribution although it hits the same object
100                // mark this using a counter
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 {
112                // unref the old ray and add the new ray into the mutation buffer
113                *mRays[mBufferStart].mRay = *vssRays[i];
114                mRays[mBufferStart].mMutations = 0;
115                mRays[mBufferStart].mUnsuccessfulMutations = 0;
116                mRays[mBufferStart].ResetReverseMutation();
117                //              mRays[mBufferStart] = RayEntry(newRay);
118                mBufferStart++;
119                if (mBufferStart >= mMaxRays)
120                  mBufferStart = 0;
121          }
122        } else {
123          // the ray did not generate any contribution
124          if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
125                  vssRays[i]->mGeneratorId != -1
126                  ) {
127                // check whether not to store a new backward mutation candidate
128                // this happens if the new ray is occluded significantly closer than
129                // the generator ray
130                VssRay *oldRay = mRays[vssRays[i]->mGeneratorId].mRay;
131                VssRay *newRay = vssRays[i];
132
133                Intersectable *oldObject = oldRay->mTerminationObject;
134
135                 
136                // only allow one reverse mutation per generator ray
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                 
142                  if (newDist < oldDist - oldObject->GetBox().Radius()*mReverseSamplesDistance) {
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               
155                Intersectable *newObject = vssRays[i]->mTerminationObject;
156
157                if (oldObject == newObject)
158                  dummyNcMutations++;
159#endif
160          }
161        }
162  }
163 
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;
170        cout<<"Reverse candidates:"<<100.0f*reverseCandidates/(float)mutationRays<<endl;
171  }
172 
173  float pContributingRays = contributingRays/(float)vssRays.size();
174 
175  cout<<"Ratio of contributing rays:"<<pContributingRays<<endl;
176 
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)
188          importance = 1.0f/(pContributingRays + 1e-5);
189        // set this values for last contributingRays
190        int index = mBufferStart - 1;
191       
192        for (int i=0; i < contributingRays; i++, index--) {
193          if (index < 0)
194                index = (int)mRays.size()-1;
195          mRays[index].mImportance = importance;
196        }
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;
204  mLastIndex = (int)mRays.size();
205  cout<<"Mutation candidates sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
206#endif
207                                                                                                                                                                 
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<<" "<<
222        GetEntry((int)mRays.size()-1).mImportance<<endl;
223
224  cout<<"Sampling factor = "<<
225        GetEntry(0).GetSamplingFactor()<<" "<<
226        GetEntry((int)mRays.size()-1).GetSamplingFactor()<<endl;
227
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
304bool
305MutationBasedDistribution::ComputeReverseMutation(
306                                                                                                  const VssRay &oldRay,
307                                                                                                  const VssRay &newRay,
308                                                                                                  Vector3 &origin,
309                                                                                                  Vector3 &termination
310                                                                                                  )
311{
312#ifdef GTP_INTERNAL
313  static int counter = 0;
314
315 
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)
325
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
334  Intersectable *occluder = newRay.mTerminationObject;
335
336
337  AxisAlignedBox3 box = occluder->GetBox();
338  // consider slightly larger neighborhood of the occluder
339  // in search for unocclude reverse ray
340  box.Scale(2.0f);
341  //box.Scale(200.0f);
342
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];
349  static Vector3 origs[packetSize];
350  // now find the silhouette along the line
351  int i;
352  float left = 0.0f;
353  float right = 1.0f;
354
355  // cast rays to find silhouette ray
356  for (int j=0; j < mSilhouetteSearchSteps; j++) {
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 );
361          origs[i] = termination;
362          // mlrtaStoreRayASEye4(&termination.x, &dirs[i].x, i);
363        }
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(),
369                                                                                          origs,
370                                                                                          dirs,
371                                                                                          hit_triangles,
372                                                                                          dist);       
373       
374        for (i=0; i < packetSize; i++) {
375          //cout<<hit_triangles[i]<<endl;
376          if (hit_triangles[i] == -1) {
377                // break on first passing ray
378            break;
379          }
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
387
388  float t = right;
389  if (right==1.0f)
390        return false;
391
392
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         
417          Intersectable *occludee =
418                oldRay.mTerminationObject;
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
435  return true;
436#else
437  cerr << "warning: reverse mutation not supported!" << endl;
438  return false;
439#endif
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
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{
456#ifdef GTP_INTERNAL
457
458  const int packetSize = 4;
459  static int hit_triangles[packetSize];
460  static float dist[packetSize];
461  static Vector3 dirs[packetSize];
462  static Vector3 shifts[packetSize];
463  static Vector3 origs[packetSize];
464  // mutate the
465  float alpha = RandomValue(0.0f, Real(2.0*M_PI));
466  //float alpha = vr2.x*2.0f*M_PI;
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
473  int i;
474  float left = 0.0f;
475  float right = radius;
476 
477  AxisAlignedBox3 _box = box;
478  _box.Scale(1.1f);
479  // cast rays to find silhouette ray
480  for (int j=0; j < mSilhouetteSearchSteps; j++) {
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 );
485          origs[i] = origin;
486          //mlrtaStoreRayASEye4(&origin.x, &dirs[i].x, i);
487        }
488       
489        // mlrtaTraverseGroupASEye4(&box.Min().x, &box.Max().x, hit_triangles, dist);
490        assert(preprocessor->mRayCaster);
491        preprocessor->mRayCaster->CastRaysPacket4(_box.Min(),
492                                                                                          _box.Max(),
493                                                                                          origs,
494                                                                                          dirs,
495                                                                                          hit_triangles,
496                                                                                          dist);       
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;
509  }
510
511  Vector3 shift;
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
516        //cout<<"W"<<i<<endl;
517        //      return (RandomValue(1.0f, 1.5f)*radius)*line;
518        shift = right*line;
519  } else {
520        //  cout<<i<<endl;
521        shift = shifts[i];
522  }
523 
524
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 
580#else
581  //cerr << "warning: shiluette mutation not supported" << endl;
582  return Vector3(0, 0, 0);
583#endif
584 
585}
586
587
588bool
589MutationBasedDistribution::GenerateSample(SimpleRay &sray)
590{
591
592  if (mRays.size() == 0) {
593        float rr[5];
594        // use direction based distribution until we have some mutation candidates
595        Vector3 origin, direction;
596
597        for (int i=0; i < 5; i++)
598          rr[i] = RandomValue(0,1);
599       
600        mPreprocessor.mViewCellsManager->GetViewPoint(origin,
601                                                                                                  Vector3(rr[0], rr[1], rr[2]));
602       
603       
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;
609       
610        return true;
611  }
612
613  int index;
614
615#if !MUTATION_USE_CDF
616#if SORT_RAY_ENTRIES
617  // RAYS are sorted -> find mitation candidate from the tail of the buffer
618  index = mLastIndex - 1;
619  if (index < 0 || index >= mRays.size()-1) {
620        index = (int)mRays.size() - 1;
621  } else
622        if (
623                mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
624          // make another round to equalize the oversampling factor
625
626          //      cout<<"R2"<<endl;
627          //      cout<<mLastIndex<<endl;
628          //      cout<<index<<endl;
629          index = (int)mRays.size() - 1;
630        }
631#else
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  }
650#endif
651#else
652
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
657  int l=0, r=(int)mRays.size()-1;
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
675  // WE HAVE THE INDEX HERE
676
677  mLastIndex = index;
678  //  Debug<<index<<" "<<mRays[index].GetSamplingFactor()<<endl;
679 
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),
684                                         MUTATION_BASED_DISTRIBUTION,
685                                         1.0f);
686       
687       
688        sray.mGeneratorId = index;
689        mRays[index].ResetReverseMutation();
690        mRays[index].mMutations++;
691        mRays[index].mUnsuccessfulMutations++;
692       
693        return true;
694  }
695
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
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);
718 
719  // mutate the origin
720  Vector3 d = ray->GetDir();
721 
722  float objectRadius = box.Radius();
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();
734
735 
736  // use probabilitistic approach to decide for the type of mutation
737  float a = RandomValue(0.0f,1.0f);
738  bool bidirectional = true;
739 
740  if (mUseSilhouetteSamples && a < mSilhouetteProb) {
741        termination += ComputeSilhouetteTerminationMutation(*ray,
742                                                                                                                origin,
743                                                                                                                box,
744                                                                                                                U, V,
745                                                                                                                2.0f*objectRadius);
746        bidirectional = false;
747  } else {
748        mRays[index].mHalton.GetNext(4, rr);
749
750        // fuzzy random mutation
751        origin += ComputeOriginMutation(*ray, U, V,
752                                                                        Vector2(rr[0], rr[1]),
753                                                                        mMutationRadiusOrigin*objectRadius);
754       
755        termination += ComputeTerminationMutation(*ray, U, V,
756                                                                                          Vector2(rr[2], rr[3]),
757                                                                                          mMutationRadiusTermination*objectRadius);
758  }
759
760  Vector3 direction = termination - origin;
761 
762  if (Magnitude(direction) < Limits::Small)
763        return false;
764 
765  // shift the origin a little bit
766#if 1
767  origin += direction*0.5f;
768#else
769  // $$JB - 14.1. 2008 increase shift to test HavranRayCaster
770  origin = 0.5f*(origin + termination);
771#endif
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;
780  sray.SetBidirectional(bidirectional);
781 
782  return true;
783}
784
785bool
786MutationBasedDistribution::GenerateMutation(const int index, SimpleRay &sray)
787{
788  VssRay *ray = mRays[index].mRay;
789
790  Intersectable *object = ray->mTerminationObject;
791 
792  AxisAlignedBox3 box = object->GetBox();
793 
794  if (GenerateMutationCandidate(index, sray, object, box)) {
795    mRays[index].mMutations++;
796        mRays[index].mUnsuccessfulMutations++;
797        return true;
798  }
799  return false;
800}
801 
802
803
804MutationBasedDistribution::MutationBasedDistribution(Preprocessor &preprocessor
805                                                                                                         ) :
806  SamplingStrategy(preprocessor)
807{
808  mType = MUTATION_BASED_DISTRIBUTION;
809  mBufferStart = 0;
810  mLastIndex = 0;
811
812  Environment::GetSingleton()->GetIntValue("Mutation.bufferSize",
813                                                                                   mMaxRays);
814 
815  mRays.reserve(mMaxRays);
816 
817  Environment::GetSingleton()->GetFloatValue("Mutation.radiusOrigin",
818                                                                                         mMutationRadiusOrigin);
819
820  Environment::GetSingleton()->GetFloatValue("Mutation.radiusTermination",
821                                                                                         mMutationRadiusTermination);
822 
823  bool useEnhancedFeatures;
824
825#ifdef GTP_INTERNAL
826  int rayCaster;
827  Environment::GetSingleton()->GetIntValue("Preprocessor.rayCastMethod",
828                                                                                   rayCaster);
829  useEnhancedFeatures = (rayCaster !=
830                                                 RayCaster::INTERNAL_RAYCASTER);
831#else
832  useEnhancedFeatures = false;
833#endif
834 
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  }
851 
852  Environment::GetSingleton()->GetBoolValue("Mutation.useSilhouetteSamples",
853                                                                                        mUseSilhouetteSamples);
854
855  Environment::GetSingleton()->GetIntValue("Mutation.silhouetteSearchSteps",
856                                                                                   mSilhouetteSearchSteps);
857 
858 
859  Environment::GetSingleton()->GetBoolValue("Mutation.usePassImportance",
860                                                                                        mUsePassImportance);
861 
862
863  Environment::GetSingleton()->GetBoolValue("Mutation.useUnsuccCountImportance",
864                                                                                        mUseUnsuccCountImportance);
865
866 
867 
868}
869
870
871
872}
Note: See TracBrowser for help on using the repository browser.