source: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp @ 2718

Revision 2718, 11.1 KB checked in by mattausch, 16 years ago (diff)
Line 
1#include "RayCaster.h"
2#include "VssRay.h"
3#include "Ray.h"
4#include "Preprocessor.h"
5#include "ViewCellsManager.h"
6
7#include <cassert>
8
9namespace GtpVisibilityPreprocessor {
10
11
12#define DEBUG_PROCESS_RAY 0
13
14#define DEBUG_RAYCAST 0
15
16#define EXACT_BOX_CLIPPING 0
17
18RayCaster::RayCaster(const Preprocessor &preprocessor):
19  mVssRayPool(), mPreprocessor(preprocessor)
20{
21}
22
23
24RayCaster::~RayCaster()
25{
26}
27
28
29VssRay *RayCaster::CastRay(const SimpleRay &simpleRay,
30                                                   const AxisAlignedBox3 &box,
31                                                   const bool castDoubleRay)
32{
33        VssRayContainer rays;
34        CastRay(simpleRay, rays, box, castDoubleRay, true);
35   
36        if (!rays.empty())
37                return rays.back();
38        else
39                return NULL;
40}
41
42
43bool
44RayCaster::ClipToViewSpaceBox(const Vector3 &origin,
45                                                          const Vector3 &termination,
46                                                          Vector3 &clippedOrigin,
47                                                          Vector3 &clippedTermination)
48{
49 
50  Ray ray(origin, termination - origin, Ray::LINE_SEGMENT);     
51  ray.Precompute();
52 
53  float tmin, tmax;
54  if ((!mPreprocessor.mViewCellsManager->
55           GetViewSpaceBox().ComputeMinMaxT(ray, &tmin, &tmax)) ||
56          tmin>=tmax
57          )
58        return false;
59
60  if (tmin >= 1.0f || tmax <= 0.0f)
61        return false;
62 
63  if (tmin > 0.0f)
64        clippedOrigin = ray.Extrap(tmin);
65  else
66        clippedOrigin = origin;
67 
68  if (tmax < 1.0f)
69        clippedTermination = ray.Extrap(tmax);
70  else
71        clippedTermination = termination;
72 
73  return true;
74}
75
76/** Checks if ray is valid
77        (e.g., not in empty view space or outside the view space)
78*/
79bool
80RayCaster::ValidateRay(const Vector3 &origin,
81                                           const Vector3 &direction,
82                                           const AxisAlignedBox3 &box,
83                                           Intersection &hit)
84{
85        if (!hit.mObject) 
86        {
87                // compute intersection with the scene bounding box
88#if EXACT_BOX_CLIPPING
89                static Ray ray;
90                mPreprocessor.SetupRay(ray, origin, direction);
91
92                float tmin, tmax;
93                if (box.ComputeMinMaxT(ray, &tmin, &tmax) && (tmin < tmax))
94                {
95                        hit.mPoint = ray.Extrap(tmax);
96                }
97                else
98                {
99                        // cout<<" invalid hp "<<tmin<<" "<<tmax<<endl;
100                        // cout<<" box "<<box<<endl;
101                        // cout<<" origin "<<origin<<endl;
102                        // cout<<" direction "<<direction<<endl;
103                        // cout<< "inv dir"<<ray.GetInvDir()<<endl;
104                        return false;
105                }
106#else
107                hit.mPoint = origin + direction * Magnitude(box.Diagonal());
108#endif
109        }
110        else
111        {
112          if (mPreprocessor.mDetectEmptyViewSpace)  {
113                if (DotProd(hit.mNormal, direction) >= -Limits::Small)  {       
114                  hit.mObject = NULL;
115                  return false;
116                }
117          }
118        }
119       
120        return true;
121}
122
123void
124RayCaster::SortRays(SimpleRayContainer &rays)
125{
126  AxisAlignedBox3 box =
127        mPreprocessor.mViewCellsManager->GetViewSpaceBox();
128 
129  float b[12]={
130        box.Min().x,
131        box.Min().y,
132        box.Min().z,
133        -1,
134        -1,
135        -1,
136        box.Max().x,
137        box.Max().y,
138        box.Max().z,
139        1,
140        1,
141        1
142  };
143
144#if 0
145  static vector<SimpleRay *> pointerArray;
146
147  if (pointerArray.size()!=rays.size()) {
148        // realloc the pointerarray
149        pointerArray.resize(rays.size());
150  }
151
152  // init pointer array
153  SimpleRay *p = &pointerArray[0];
154  for (i=0; i < rays.size(); i++, p++)
155        pointerArray[i] = p;
156#endif
157 
158#if 1
159  _SortRays(rays,
160                        0,
161                        (int)rays.size()-1,
162                        0,
163                        b
164                        );
165#endif
166}
167                                       
168void
169RayCaster::_SortRays(SimpleRayContainer &rays,
170                                         const int l,
171                                         const int r,
172                                         const int depth,
173                                         float *box)
174{
175  // pick-up a pivot
176  int axis = 0;
177 
178  float maxDiff = -1.0f;
179  // get the largest axis
180  int offset = 0;
181  int i;
182
183  //const int batchsize = 16384;
184  const int batchsize = 8192;
185  //const int batchsize = 128;
186
187  //if (r - l < 16*batchsize)
188  //    offset = 3;
189       
190  //  if (depth%2==0)
191  //    offset = 3;
192 
193  for (i=offset; i < offset + 3; i++) {
194        float diff = box[i + 6] - box[i];
195        if (diff > maxDiff) {
196          maxDiff = diff;
197          axis = i;
198        }
199  }
200
201  //  cout<<depth<<" "<<axis<<" "<<l<<" "<<r<<endl;
202 
203  i=l;
204  int j=r;
205
206  float x = (box[axis] + box[axis+6])*0.5f;
207  //  float x = rays[(l+r)/2].GetParam(axis);
208  do {
209        while(i<j && rays[i].GetParam(axis) < x)
210          i++;
211        while(i<j && x < rays[j].GetParam(axis))
212          j--;
213       
214        if (i <= j) {
215          swap(rays[i], rays[j]);
216          i++;
217          j--;
218        }
219  } while (i<=j);
220
221 
222  if (l + batchsize < j ) {
223        // set new max
224        float save = box[axis+6];
225        box[axis+6] = x;
226        _SortRays(rays, l, j, depth+1, box);
227        box[axis+6] = save;
228  } else {
229        //      for (int k=0; k < 6; k++)
230        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
231  }
232 
233  if (i + batchsize < r) {
234        // set new min
235        box[axis] = x;
236    _SortRays(rays, i, r, depth+1, box);
237  } else {
238        //      for (int k=0; k < 6; k++)
239        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
240  }
241       
242}
243 
244
245void
246RayCaster::SortRays2(SimpleRayContainer &rays)
247{
248  AxisAlignedBox3 box =
249        mPreprocessor.mViewCellsManager->GetViewSpaceBox();
250
251  const float sizeBox = Magnitude(box.Diagonal());
252  // This is some size of the
253  const float sizeDir = 0.2f * sizeBox;
254 
255  float b[12]={
256    box.Min().x,
257    box.Min().y,
258    box.Min().z,
259    -sizeDir,
260    -sizeDir,
261    -sizeDir,
262    box.Max().x,
263    box.Max().y,
264    box.Max().z,
265    sizeDir,
266    sizeDir,
267    sizeDir
268  };
269
270#if 0
271  static vector<SimpleRay *> pointerArray;
272
273  if (pointerArray.size()!=rays.size()) {
274        // realloc the pointerarray
275        pointerArray.resize(rays.size());
276  }
277
278  // init pointer array
279  SimpleRay *p = &pointerArray[0];
280  for (i=0; i < rays.size(); i++, p++)
281        pointerArray[i] = p;
282#endif
283 
284  _SortRays2(rays, 0, (int)rays.size()-1, 0, b);
285
286  return;
287}
288                                       
289void
290RayCaster::_SortRays2(SimpleRayContainer &rays,
291                      const int l,
292                      const int r,
293                      const int depth,
294                      float box[12])
295{
296  // pick-up a pivot
297  int axis;
298 
299  float maxDiff = -1.0f;
300  // get the largest axis
301  int offset = 0;
302  int i;
303
304  //const int batchsize = 16384;
305  //const int batchsize = 8192;
306  const int batchsize = 128;
307
308  //if (r - l < 16*batchsize)
309  //    offset = 3;
310       
311  //  if (depth%2==0)
312  //    offset = 3;
313 
314  for (i=offset; i < offset + 6; i++) {
315    float diff = box[i + 6] - box[i];
316    assert(diff >= 0.f);
317    if (diff > maxDiff) {
318      // Find the maximum
319      maxDiff = diff;
320      axis = i;
321    }
322  }
323
324  //  cout<<depth<<" "<<axis<<" "<<l<<" "<<r<<endl;
325 
326  i=l;
327  int j=r;
328
329  float x = (box[axis] + box[axis+6])*0.5f;
330  //  float x = rays[(l+r)/2].GetParam(axis);
331  do {
332        while(i<j && rays[i].GetParam(axis) < x)
333          i++;
334        while(i<j && x < rays[j].GetParam(axis))
335          j--;
336       
337        if (i <= j) {
338          swap(rays[i], rays[j]);
339          i++;
340          j--;
341        }
342  } while (i<=j);
343
344 
345  if (l + batchsize < j ) {
346        // set new max
347        float save = box[axis+6];
348        box[axis+6] = x;
349        _SortRays2(rays, l, j, depth+1, box);
350        box[axis+6] = save;
351  } else {
352        //      for (int k=0; k < 6; k++)
353        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
354  }
355 
356  if (i + batchsize < r) {
357    // set new min
358    box[axis] = x;
359    _SortRays2(rays, i, r, depth+1, box);
360  } else {
361        //      for (int k=0; k < 6; k++)
362        //        cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
363  }
364       
365}
366 
367 
368VssRay *RayCaster::RequestRay(const Vector3 &origin,
369                                                          const Vector3 &termination,
370                                                          Intersectable *originObject,
371                                                          Intersectable *terminationObject,
372                                                          const int pass,
373                                                          const float pdf)
374{
375 #if DEBUG_RAYCAST
376        Debug<<"PR2a"<<flush;
377#endif
378       
379        // old method: always allocate
380        if (0) return new VssRay(origin, termination, originObject, terminationObject, pass, pdf);
381
382        VssRay *vssRay = mVssRayPool.Alloc();
383
384#if DEBUG_RAYCAST
385        Debug<<"PR2b"<<flush;
386#endif
387       
388        *vssRay = VssRay(origin, termination, originObject, terminationObject, pass, pdf);
389
390#if DEBUG_RAYCAST
391        Debug<<"PR2c"<<flush;
392#endif
393
394        return vssRay;
395}
396
397
398int
399RayCaster::ProcessRay(const SimpleRay &simpleRay,
400                                          Intersection &hitA,
401                                          Intersection &hitB,
402                                          VssRayContainer &vssRays,
403                                          const AxisAlignedBox3 &box,
404                                          const bool castDoubleRay,
405                                          const bool pruneInvalidRays)
406{
407  int hits = 0;
408
409
410#if DEBUG_RAYCAST
411  static int id=0;
412  Debug<<"PRA "<<id++<<endl<<flush;
413#endif
414       
415        if (pruneInvalidRays)
416        {
417          if (!hitA.mObject && !hitB.mObject) {
418#if DEBUG_PROCESS_RAY
419                cout<<"I1 ";
420#endif
421                return 0;
422          }
423        }
424       
425        // regardless of the pruneInvalidRays setting reject
426        // rays whic degenerate to a point
427        if (EpsilonEqualV3(hitA.mPoint, hitB.mPoint, Limits::Small)) {
428#if DEBUG_PROCESS_RAY
429                cout<<"I2 ";
430#endif
431          return 0;
432        }
433       
434        const bool validA = ValidateRay(simpleRay.mOrigin, simpleRay.mDirection, box, hitA);
435        const bool validB = //castDoubleRay &&
436          ValidateRay(simpleRay.mOrigin, -simpleRay.mDirection, box, hitB);
437       
438       
439#if DEBUG_RAYCAST
440        Debug<<"PR1"<<flush;
441#endif
442       
443        // reset both contributions
444        if (!validA || !validB) {
445          if (0 || pruneInvalidRays)
446                return 0;
447
448#if DEBUG_PROCESS_RAY
449                cout<<"I2 ";
450#endif
451
452          // reset both contributions of this ray
453          hitA.mObject = NULL;
454          hitB.mObject = NULL;
455        }
456       
457        // 8.11. 2007 JB
458        // degenerate rays checked by geometrical constraint...
459        //      !pruneInvalidRays || (hitA.mObject != hitB.mObject);
460
461       
462#if DEBUG_RAYCAST
463        Debug<<"PR2"<<flush;
464#endif
465
466        // VH - I do not know what this is for, commented out temporarily
467        // !!!!!!!!!!!!!! !!!!!!!!!!!! !!!!!!!!!!
468#if 1
469        const bool validSample = true;
470        if (validSample) {
471          Vector3 clipA, clipB;
472          if (!ClipToViewSpaceBox(hitA.mPoint,
473                                                          hitB.mPoint,
474                                                          clipA,
475                                                          clipB)) {
476#if DEBUG_PROCESS_RAY
477                cout<<"I3 ";
478#endif
479
480                return 0;
481          }
482
483          if (!pruneInvalidRays || hitA.mObject) {
484
485                  VssRay *vssRay =
486                          RequestRay(!castDoubleRay ? simpleRay.mOrigin : clipB,
487                                                 hitA.mPoint,
488                                                 hitB.mObject,
489                                                 hitA.mObject,
490                                                 mPreprocessor.mPass,
491                                                 1.0f //simpleRay.mPdf
492                                                 );
493
494                if (validA)
495                  vssRay->mFlags |= VssRay::Valid;
496               
497                vssRay->mDistribution = simpleRay.mDistribution;
498                vssRay->mGeneratorId = simpleRay.mGeneratorId;
499
500                vssRays.push_back(vssRay);
501                ++ hits;
502                //cout << "vssray 1: " << *vssRay << " " << vssRay->mTermination - vssRay->mOrigin << endl;
503        }
504
505#if DEBUG_RAYCAST
506        Debug<<"PR3"<<flush;
507#endif
508
509        if (castDoubleRay && (!pruneInvalidRays || hitB.mObject))
510        {
511                VssRay *vssRay = RequestRay(
512                                            clipA,
513                                            hitB.mPoint,
514                                            hitA.mObject,
515                                            hitB.mObject,
516                                            mPreprocessor.mPass,
517                                            1.0f //simpleRay.mPdf
518                                            );
519
520                if (validB)
521                        vssRay->mFlags |= VssRay::Valid;
522
523                vssRay->mDistribution = simpleRay.mDistribution;
524                vssRay->mGeneratorId = simpleRay.mGeneratorId;
525                vssRays.push_back(vssRay);
526                ++ hits;
527                //cout << "vssray 2: " << *vssRay << endl;
528        }
529#if DEBUG_RAYCAST
530          Debug<<"PR4"<<flush;
531#endif
532        } // validSample
533#else
534        // Just pass if intersected or not
535        hits = (hitA.mObject != 0) ? 1 : 0;
536        intersect = hitA;
537#endif 
538       
539        return hits;
540}
541
542
543void
544RayCaster::CastRays(
545                    SimpleRayContainer &rays,
546                    VssRayContainer &vssRays,
547                    const AxisAlignedBox3 &sbox,
548                    const bool castDoubleRay,
549                    const bool pruneInvalidRays )
550{
551  SimpleRayContainer::const_iterator rit, rit_end = rays.end();
552 
553  for (rit = rays.begin(); rit != rit_end; ++ rit) {
554   
555    CastRay(
556            *rit,                               
557            vssRays,
558            sbox,
559            castDoubleRay,
560            pruneInvalidRays);
561        }
562}
563 
564
565}
Note: See TracBrowser for help on using the repository browser.