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

Revision 2011, 21.9 KB checked in by bittner, 17 years ago (diff)

mutation strategy updated, importance depending on rsuccess

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