source: GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp @ 1942

Revision 1942, 15.4 KB checked in by bittner, 18 years ago (diff)

tmp commit

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
9namespace GtpVisibilityPreprocessor {
10
11
12//HaltonSequence SamplingStrategy::sHalton;
13
14HaltonSequence ObjectBasedDistribution::sHalton;
15HaltonSequence MixtureDistribution::sHalton;
16HaltonSequence GlobalLinesDistribution::sHalton;
17HaltonSequence SpatialBoxBasedDistribution::sHalton;
18HaltonSequence ObjectDirectionBasedDistribution::sHalton;
19
20
21SamplingStrategy::SamplingStrategy(Preprocessor &preprocessor):
22  mPreprocessor(preprocessor), mRatio(1.0f),  mTotalRays(0), mTotalContribution(0.0f)
23{
24 
25}
26
27
28SamplingStrategy::~SamplingStrategy()
29{
30}
31
32int
33SamplingStrategy::GenerateSamples(const int number,
34                                                                  SimpleRayContainer &rays)
35{
36        SimpleRay ray;
37        int samples = 0;
38        int i = 0;
39        const int maxTries = 20;
40        // tmp changed matt. Q: should one rejected sample
41        // terminate the whole method?
42        if (0)
43        {
44                for (; i < number; i++) {
45                        if (!GenerateSample(ray))
46                                return i;
47                        rays.push_back(ray);
48                }
49        }
50        else
51        {
52                for (; i < number; i++)
53                {
54                        int j = 0;
55                        bool sampleGenerated = false;
56
57                        for (j = 0; !sampleGenerated && (j < maxTries); ++ j)
58                        {
59                                sampleGenerated = GenerateSample(ray);
60
61                                if (sampleGenerated)
62                                {               
63                                        ++ samples;
64                                        rays.push_back(ray);
65                                }
66                        }
67                }
68        }
69       
70        return samples;
71}
72
73
74/*********************************************************************/
75/*            Individual sampling strategies implementation          */
76/*********************************************************************/
77
78
79bool ObjectBasedDistribution::GenerateSample(SimpleRay &ray)
80{
81  Vector3 origin, direction;
82 
83  float r[5];
84  sHalton.GetNext(5, r);
85 
86  mPreprocessor.mViewCellsManager->GetViewPoint(origin,
87                                                                                                Vector3(r[2],r[3],r[4]));
88 
89
90  Vector3 point, normal;
91 
92  r[0] *= (float)mPreprocessor.mObjects.size() - 1;
93  const int i = (int)r[0];
94 
95  Intersectable *object = mPreprocessor.mObjects[i];
96 
97  // take the remainder as a parameter over objects surface
98  r[0] -= (float)i;
99 
100  object->GetRandomSurfacePoint(r[0], r[1], point, normal);
101 
102  direction = point - origin;
103 
104  const float c = Magnitude(direction);
105 
106  if (c <= Limits::Small)
107        return false;
108 
109  // $$ jb the pdf is yet not correct for all sampling methods!
110  const float pdf = 1.0f;
111 
112  direction *= 1.0f / c;
113  ray = SimpleRay(origin, direction, OBJECT_BASED_DISTRIBUTION, pdf);
114 
115  return true;
116}
117
118
119bool
120ObjectDirectionBasedDistribution::GenerateSample(SimpleRay &ray)
121{       
122  Vector3 origin, direction;
123
124
125  float r[4];
126  sHalton.GetNext(4, r);
127
128  r[0] *= (float)mPreprocessor.mObjects.size() - 1;
129  const int i = (int)r[0];
130 
131  Intersectable *object = mPreprocessor.mObjects[i];
132 
133  // take the remainder as a parameter over objects surface
134  r[0] -= (float)i;
135
136  Vector3 normal;
137
138  object->GetRandomSurfacePoint(r[0], r[1], origin, normal);
139 
140  direction = Normalize(CosineRandomVector(r[2], r[3], normal));
141 
142  origin += 1e-2f*direction;
143 
144  // $$ jb the pdf is yet not correct for all sampling methods!
145  const float pdf = 1.0f;
146 
147  ray = SimpleRay(origin, direction, OBJECT_DIRECTION_BASED_DISTRIBUTION, pdf);
148 
149  return true;
150}
151
152
153bool DirectionBasedDistribution::GenerateSample(SimpleRay &ray)
154{       
155        Vector3 origin, direction;
156        mPreprocessor.mViewCellsManager->GetViewPoint(origin);
157         
158        direction = UniformRandomVector();
159        const float c = Magnitude(direction);
160
161        if (c <= Limits::Small)
162                return false;
163
164        const float pdf = 1.0f;
165
166        direction *= 1.0f / c;
167        ray = SimpleRay(origin, direction, DIRECTION_BASED_DISTRIBUTION, pdf);
168
169        return true;
170}
171
172
173bool DirectionBoxBasedDistribution::GenerateSample(SimpleRay &ray)
174{
175        Vector3 origin, direction;
176        mPreprocessor.mViewCellsManager->GetViewPoint(origin);
177
178        const float alpha = RandomValue(0.0f, 2.0f * (float)M_PI);
179        const float beta = RandomValue((float)-M_PI * 0.5f, (float)M_PI * 0.5f);
180       
181        direction = VssRay::GetDirection(alpha, beta);
182       
183        const float c = Magnitude(direction);
184
185        if (c <= Limits::Small)
186                return false;
187
188        const float pdf = 1.0f;
189
190        direction *= 1.0f / c;
191        ray = SimpleRay(origin, direction, DIRECTION_BOX_BASED_DISTRIBUTION, pdf);
192
193        return true;
194}
195
196
197bool SpatialBoxBasedDistribution::GenerateSample(SimpleRay &ray)
198{
199  Vector3 origin, direction;
200
201  float r[6];
202  sHalton.GetNext(6, r);
203  mPreprocessor.mViewCellsManager->GetViewPoint(origin, Vector3(r[0],r[1],r[2]));
204 
205  direction = mPreprocessor.mKdTree->GetBox().GetRandomPoint(Vector3(r[3],
206                                                                                                                                         r[4],
207                                                                                                                                         r[5])
208                                                                                                                                         ) - origin;
209  //cout << "z";
210  const float c = Magnitude(direction);
211 
212  if (c <= Limits::Small)
213        return false;
214 
215  const float pdf = 1.0f;
216 
217  direction *= 1.0f / c;
218  ray = SimpleRay(origin, direction, SPATIAL_BOX_BASED_DISTRIBUTION, pdf);
219 
220  return true;
221}
222
223
224bool ReverseObjectBasedDistribution::GenerateSample(SimpleRay &ray)
225{
226        Vector3 origin, direction;
227
228        mPreprocessor.mViewCellsManager->GetViewPoint(origin);
229
230        Vector3 point;
231        Vector3 normal;
232        //cout << "y";
233        const int i = (int)RandomValue(0, (float)mPreprocessor.mObjects.size() - 0.5f);
234
235        Intersectable *object = mPreprocessor.mObjects[i];
236
237        object->GetRandomSurfacePoint(point, normal);
238       
239        direction = origin - point;
240       
241        // $$ jb the pdf is yet not correct for all sampling methods!
242        const float c = Magnitude(direction);
243       
244        if ((c <= Limits::Small) || (DotProd(direction, normal) < 0))
245        {
246                return false;
247        }
248
249        // $$ jb the pdf is yet not correct for all sampling methods!
250        const float pdf = 1.0f;
251        //cout << "p: " << point << " ";
252        direction *= 1.0f / c;
253        // a little offset
254        point += direction * 0.001f;
255
256        ray = SimpleRay(point, direction, REVERSE_OBJECT_BASED_DISTRIBUTION, pdf);
257       
258        return true;
259}
260
261
262bool ViewCellBorderBasedDistribution::GenerateSample(SimpleRay &ray)
263{
264        Vector3 origin, direction;
265
266        ViewCellContainer &viewCells = mPreprocessor.mViewCellsManager->GetViewCells();
267
268        Vector3 point;
269        Vector3 normal, normal2;
270       
271        const int vcIdx = (int)RandomValue(0, (float)viewCells.size() - 0.5f);
272        const int objIdx = (int)RandomValue(0, (float)mPreprocessor.mObjects.size() - 0.5f);
273
274        Intersectable *object = mPreprocessor.mObjects[objIdx];
275        ViewCell *viewCell = viewCells[vcIdx];
276
277        //cout << "vc: " << vcIdx << endl;
278        //cout << "obj: " << objIdx << endl;
279
280        object->GetRandomSurfacePoint(point, normal);
281        viewCell->GetRandomEdgePoint(origin, normal2);
282
283        direction = point - origin;
284
285        // $$ jb the pdf is yet not correct for all sampling methods!
286        const float c = Magnitude(direction);
287
288        if ((c <= Limits::Small) /*|| (DotProd(direction, normal) < 0)*/)
289        {
290                return false;
291        }
292
293        // $$ jb the pdf is yet not correct for all sampling methods!
294        const float pdf = 1.0f;
295        //cout << "p: " << point << " ";
296        direction *= 1.0f / c;
297        ray = SimpleRay(origin, direction, VIEWCELL_BORDER_BASED_DISTRIBUTION, pdf);
298
299        //cout << "ray: " << ray.mOrigin << " " << ray.mDirection << endl;
300
301        return true;
302}
303
304
305#if 0
306bool ObjectsInteriorDistribution::GenerateSample(SimpleRay &ray)
307{
308        Vector3 origin, direction;
309
310        // get random object
311        const int i = RandomValue(0, mPreprocessor.mObjects.size() - 1);
312
313        const Intersectable *obj = mPreprocessor.mObjects[i];
314
315        // note: if we load the polygons as meshes,
316        // asymtotically every second sample is lost!
317        origin = obj->GetBox().GetRandomPoint();
318
319        // uniformly distributed direction
320        direction = UniformRandomVector();
321
322        const float c = Magnitude(direction);
323
324        if (c <= Limits::Small)
325                return false;
326
327        const float pdf = 1.0f;
328
329        direction *= 1.0f / c;
330        ray = SimpleRay(origin, direction, pdf);
331
332        return true;
333}
334
335#endif
336
337
338bool ReverseViewSpaceBorderBasedDistribution::GenerateSample(SimpleRay &ray)
339{
340        Vector3 origin, direction;
341
342        origin = mPreprocessor.mViewCellsManager->GetViewSpaceBox().GetRandomSurfacePoint();
343
344        Vector3 point;
345        Vector3 normal;
346        //cout << "y";
347        const int i = (int)RandomValue(0, (float)mPreprocessor.mObjects.size() - 0.5f);
348
349        Intersectable *object = mPreprocessor.mObjects[i];
350
351        object->GetRandomSurfacePoint(point, normal);
352       
353        direction = origin - point;
354       
355        // $$ jb the pdf is yet not correct for all sampling methods!
356        const float c = Magnitude(direction);
357       
358        if ((c <= Limits::Small) || (DotProd(direction, normal) < 0))
359        {
360                return false;
361        }
362
363        // $$ jb the pdf is yet not correct for all sampling methods!
364        const float pdf = 1.0f;
365        //cout << "p: " << point << " ";
366        direction *= 1.0f / c;
367        // a little offset
368        point += direction * 0.001f;
369
370        ray = SimpleRay(point, direction, REVERSE_VIEWSPACE_BORDER_BASED_DISTRIBUTION, pdf);
371       
372        return true;
373}
374
375
376bool ViewSpaceBorderBasedDistribution::GenerateSample(SimpleRay &ray)
377{
378        Vector3 origin, direction;
379
380        origin = mPreprocessor.mViewCellsManager->GetViewSpaceBox().GetRandomSurfacePoint();
381
382        Vector3 point;
383        Vector3 normal;
384        //cout << "w";
385        const int i = (int)RandomValue(0, (float)mPreprocessor.mObjects.size() - 0.5f);
386
387        Intersectable *object = mPreprocessor.mObjects[i];
388
389        object->GetRandomSurfacePoint(point, normal);
390        direction = point - origin;
391
392        // $$ jb the pdf is yet not correct for all sampling methods!
393        const float c = Magnitude(direction);
394
395        if (c <= Limits::Small)
396                return false;
397
398        // $$ jb the pdf is yet not correct for all sampling methods!
399        const float pdf = 1.0f;
400       
401        direction *= 1.0f / c;
402
403        // a little offset
404        origin += direction * 0.001f;
405
406        ray = SimpleRay(origin, direction, VIEWSPACE_BORDER_BASED_DISTRIBUTION, pdf);
407
408        return true;
409}
410
411
412
413
414bool
415GlobalLinesDistribution::GenerateSample(SimpleRay &ray)
416{
417  Vector3 origin, termination, direction;
418
419  float radius = 0.5f*Magnitude(mPreprocessor.mViewCellsManager->GetViewSpaceBox().Size());
420  Vector3 center = mPreprocessor.mViewCellsManager->GetViewSpaceBox().Center();
421
422  const int tries = 1000;
423  int i;
424  for (i=0; i < tries; i++) {
425        float r[4];
426        sHalton.GetNext(4, r);
427       
428        origin = center + (radius*UniformRandomVector(r[0], r[1]));
429        termination = center + (radius*UniformRandomVector(r[2], r[3]));
430       
431        direction = termination - origin;
432       
433       
434        // $$ jb the pdf is yet not correct for all sampling methods!
435        const float c = Magnitude(direction);
436        if (c <= Limits::Small)
437          return false;
438       
439        direction *= 1.0f / c;
440
441        // check if the ray intersects the view space box
442        static Ray ray;
443        ray.Init(origin, direction, Ray::LOCAL_RAY);   
444       
445        float tmin, tmax;
446        if (mPreprocessor.mViewCellsManager->
447                GetViewSpaceBox().ComputeMinMaxT(ray, &tmin, &tmax) && (tmin < tmax))
448          break;
449  }
450 
451  if (i!=tries) {
452        // $$ jb the pdf is yet not correct for all sampling methods!
453        const float pdf = 1.0f;
454       
455       
456        ray = SimpleRay(origin, direction, GLOBAL_LINES_DISTRIBUTION, pdf);
457        ray.mType = Ray::GLOBAL_RAY;
458        return true;
459  }
460 
461  return false;
462}
463
464
465  // has to called before first usage
466void
467MixtureDistribution::Init()
468{
469  for (int i=0; i < mDistributions.size(); i++) {
470        // small non-zero value
471        mDistributions[i]->mRays = 1;
472        // unit contribution per ray
473        if (1 || mDistributions[i]->mType != RSS_BASED_DISTRIBUTION)
474          mDistributions[i]->mContribution = 1.0f;
475        else
476          mDistributions[i]->mContribution = 0.0f;
477  }
478  UpdateRatios();
479}
480
481void
482MixtureDistribution::Reset()
483{
484  for (int i=0; i < mDistributions.size(); i++) {
485        // small non-zero value
486        mDistributions[i]->mTotalRays = 0;
487        // unit contribution per ray
488        mDistributions[i]->mTotalContribution = 0.0f;
489  }
490  UpdateRatios();
491}
492
493// Generate a new sample according to a mixture distribution
494bool
495MixtureDistribution::GenerateSample(SimpleRay &ray)
496{
497  float r;
498  sHalton.GetNext(1, &r);
499  int i;
500  // pickup a distribution
501  for (i=0; i < mDistributions.size()-1; i++)
502        if (r < mDistributions[i]->mRatio)
503          break;
504
505  return mDistributions[i]->GenerateSample(ray);
506}
507
508  // add contributions of the sample to the strategies
509void
510MixtureDistribution::ComputeContributions(VssRayContainer &vssRays)
511{
512  int i;
513 
514  VssRayContainer::iterator it = vssRays.begin();
515
516  for (i=0; i < mDistributions.size(); i++) {
517        mDistributions[i]->mContribution = 0;
518        mDistributions[i]->mRays = 0;
519  }
520
521  for(; it != vssRays.end(); ++it) {
522        VssRay *ray = *it;
523        for (i=0; i < mDistributions.size()-1; i++) {
524          if (mDistributions[i]->mType == ray->mDistribution)
525                break;
526        }
527 
528        float contribution =
529          mPreprocessor.mViewCellsManager->ComputeSampleContribution(*ray,
530                                                                                                                                 true,
531                                                                                                                                 false);
532
533        mDistributions[i]->mContribution += contribution;
534        mDistributions[i]->mRays ++;
535       
536        mDistributions[i]->mTotalContribution += contribution;
537        mDistributions[i]->mTotalRays ++;
538  }
539
540 
541  UpdateRatios();
542}
543
544void
545MixtureDistribution::UpdateDistributions(VssRayContainer &vssRays)
546{
547  // now update the distributions with all the rays
548  for (int i=0; i < mDistributions.size(); i++) {
549        mDistributions[i]->Update(vssRays);
550  }
551}
552 
553void
554MixtureDistribution::UpdateRatios()
555{
556  // now compute importance (ratio) of all distributions
557  float sum = 0.0f;
558  int i;
559  for (i=0; i < mDistributions.size(); i++) {
560        cout<<i<<": c="<<mDistributions[i]->mContribution<<
561          " rays="<<mDistributions[i]->mRays<<endl;
562        float importance = 0.0f;
563        if (mDistributions[i]->mRays != 0) {
564          //      importance = pow(mDistributions[i]->mContribution/mDistributions[i]->mRays, 3);
565          importance = mDistributions[i]->mContribution/mDistributions[i]->mRays;
566        }
567        mDistributions[i]->mRatio = importance;
568        sum += importance;
569  }
570
571  const float minratio = 0.02f;
572  float threshold = minratio*sum;
573  for (i=0; i < mDistributions.size(); i++) {
574        if (mDistributions[i]->mRatio < threshold) {
575          sum += threshold - mDistributions[i]->mRatio;
576          mDistributions[i]->mRatio = threshold;
577        }
578  }
579 
580  mDistributions[0]->mRatio /= sum;
581       
582  for (i=1; i < mDistributions.size(); i++) {
583        float r = mDistributions[i]->mRatio / sum;
584        mDistributions[i]->mRatio = mDistributions[i-1]->mRatio + r;
585  }
586 
587  cout<<"ratios: ";
588  float last = 0.0f;
589  for (i=0; i < mDistributions.size(); i++) {
590        cout<<mDistributions[i]->mRatio-last<<" ";
591        last = mDistributions[i]->mRatio;
592  }
593  cout<<endl;
594}
595
596
597
598bool
599MixtureDistribution::Construct(char *str)
600{
601  char *curr = str;
602
603  while (1) {
604        char *e = strchr(curr,'+');
605        if (e!=NULL) {
606          *e=0;
607        }
608       
609        if (strcmp(curr, "rss")==0) {
610          mDistributions.push_back(new RssBasedDistribution(mPreprocessor));
611        } else
612          if (strcmp(curr, "object")==0) {
613                mDistributions.push_back(new ObjectBasedDistribution(mPreprocessor));
614          } else
615                if (strcmp(curr, "spatial")==0) {
616                  mDistributions.push_back(new SpatialBoxBasedDistribution(mPreprocessor));
617                } else
618                  if (strcmp(curr, "global")==0) {
619                        mDistributions.push_back(new GlobalLinesDistribution(mPreprocessor));
620                  } else
621                        if (strcmp(curr, "direction")==0) {
622                          mDistributions.push_back(new DirectionBasedDistribution(mPreprocessor));
623                        } else
624                          if (strcmp(curr, "object_direction")==0) {
625                                mDistributions.push_back(new ObjectDirectionBasedDistribution(mPreprocessor));
626                          } else
627                                if (strcmp(curr, "reverse_object")==0) {
628                                  mDistributions.push_back(new ReverseObjectBasedDistribution(mPreprocessor));
629                                } else
630                                  if (strcmp(curr, "reverse_viewspace_border")==0) {
631                                        mDistributions.push_back(new ReverseViewSpaceBorderBasedDistribution(mPreprocessor));
632                                  }
633       
634        if (e==NULL)
635          break;
636        curr = e+1;
637  }
638
639  Init();
640  return true;
641}
642
643}
644
645
Note: See TracBrowser for help on using the repository browser.