source: GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp @ 2588

Revision 2588, 11.7 KB checked in by bittner, 16 years ago (diff)

sceneBox bugfix
Line 
1
2// Written by Vlastimil Havran, December 2007
3
4// This macro allows to use the ray shooting written by
5// Vlastimil Havran, 2007-2008
6
7
8#include "HavranRayCaster.h"
9#include "VssRay.h"
10#include "KdTree.h"
11#include "Preprocessor.h"
12#include "IntersectableWrapper.h"
13
14#ifdef USE_HAVRAN_RAYCASTER
15#include "timer.h"
16#include "ktball.h"
17#include "ktb.h"
18#include "raypack.h"
19#endif // USE_HAVRAN_RAYCASTER
20
21#define DEBUG_RAYCAST 0
22
23namespace GtpVisibilityPreprocessor {
24
25
26#ifdef USE_HAVRAN_RAYCASTER
27// static rays
28RayPacket2x2
29HavranRayCaster::raypack;
30#endif // USE_HAVRAN_RAYCASTER
31 
32HavranRayCaster::HavranRayCaster(const Preprocessor &preprocessor):
33  RayCaster(preprocessor), mKtbtree(0)
34{
35#ifdef USE_HAVRAN_RAYCASTER
36  mKtbtree = new CKTB();
37#endif 
38}
39 
40
41HavranRayCaster::~HavranRayCaster()
42{
43#ifdef USE_HAVRAN_RAYCASTER
44  delete mKtbtree;
45  mKtbtree = 0;
46#endif // USE_HAVRAN_RAYCASTER
47}
48
49void
50HavranRayCaster::Build(ObjectContainer &objlist)
51{
52#ifdef USE_HAVRAN_RAYCASTER
53
54  CTimer timer;
55
56  cout << "\nBuilding up kd-tree for Havran ray caster ..."<<endl<<flush;
57
58  timer.Start();
59  mKtbtree->BuildUp(objlist);
60  timer.Stop();
61  cout << "\nBuilding up kd-tree is finished, user time = "
62       << timer.UserTime() << " real time = " << timer.RealTime() <<
63        endl <<flush;
64#endif
65}
66
67
68// Using packet of 4 rays supposing that these are coherent
69// We give a box to which each ray is clipped to before the
70// ray shooting is computed !
71void HavranRayCaster::CastRaysPacket4(Vector3 &boxmax,
72                                                                          Vector3 &boxmin,
73                                                                          Vector3 origin4[],
74                                                                          Vector3 direction4[],
75                                                                          int     result4[],
76                                                                          float   dist4[])
77{
78#ifdef USE_HAVRAN_RAYCASTER
79  for (int i = 0; i < 4; i++) {
80    result4[i] = -1;
81    raypack.SetLoc(i, origin4[i]);
82    raypack.SetDir(i, direction4[i]);
83  }
84 
85  // The same operations for packets of rays, if implemented by
86  // a particular ASDS, otherwise it is emulated by decomposition
87  // of a packet to individual rays and traced individually.
88  mKtbtree->FindNearestI(raypack, boxmin, boxmax);
89
90  for (int i = 0; i < 4; i++) {
91    // if ray intersects an object, we set the pointer to
92    // this object
93    Intersectable* intersectable;
94    if ( (intersectable = raypack.GetObject(i)) )
95          // $$JB this is very hacky - > should be replaced with fetching the index of the triangle
96          result4[i] = (int)intersectable;
97  }
98 
99  return;
100#endif
101}
102
103int HavranRayCaster::CastRay(const SimpleRay &simpleRay,
104                             VssRayContainer &vssRays,
105                             const AxisAlignedBox3 &box,
106                             const bool castDoubleRay,
107                             const bool pruneInvalidRays)
108{
109#ifdef USE_HAVRAN_RAYCASTER
110
111  int hits = 0;
112  Intersection hitA(simpleRay.mOrigin), hitB(simpleRay.mOrigin);
113 
114  // inside test for bounding box
115  // enlarge box slightly so the view point fits for sure
116  //  AxisAlignedBox3 sbox = box;
117  //  sbox.Enlarge(Vector3(-Limits::Small));
118  // $$ JB moved here from Validate routine
119
120//   if (!box.IsInside(simpleRay.mOrigin)) {
121//      cout<<"out of box "<<simpleRay.mOrigin<<" "<<box<<endl;
122//      return 0;
123//   }
124 
125  // ray.mFlags &= ~Ray::CULL_BACKFACES;
126
127  if (mKtbtree->FindNearestI(simpleRay))  {
128    hitA.mObject = SimpleRay::IntersectionRes[0].intersectable;
129    float tdist = SimpleRay::IntersectionRes[0].tdist;
130    hitA.mPoint = simpleRay.Extrap(tdist);
131    hitA.mNormal = SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
132    //hitA.mNormal = (dynamic_cast< TriangleIntersectable *>
133    //              (hitA.mObject))->GetNormal(0);
134  }
135 
136 
137  if (castDoubleRay) {
138    Vector3 *v = (Vector3*)(&simpleRay.mDirection);
139    *v = -(*v);
140    // ray.mFlags &= ~Ray::CULL_BACKFACES;
141    if (mKtbtree->FindNearestI(simpleRay))  {
142      hitB.mObject = SimpleRay::IntersectionRes[0].intersectable;
143      float tdist = SimpleRay::IntersectionRes[0].tdist;
144      hitB.mPoint = simpleRay.Extrap(tdist);
145      hitB.mNormal = SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
146      //hitB.mNormal = (dynamic_cast< TriangleIntersectable *>
147      //              (hitB.mObject))->GetNormal(0);
148    }
149    // restore the direction to the original
150    *v = -simpleRay.mDirection;
151  }
152
153  // This code is also in IntelRayCaster.cpp
154  return ProcessRay(
155                    simpleRay,
156                    hitA,
157                    hitB,
158                    vssRays,
159                    box,
160                    castDoubleRay,
161                    pruneInvalidRays
162                    );
163#else // USE_HAVRAN_RAYCASTER
164  return 0;
165#endif
166}
167
168void HavranRayCaster::CastRays16(SimpleRayContainer &rays,
169                                 VssRayContainer &vssRays,
170                                 const AxisAlignedBox3 &sbox,
171                                 const bool castDoubleRay,
172                                 const bool pruneInvalidRays)
173{
174  CastRays16(rays, 0, vssRays, sbox, castDoubleRay, pruneInvalidRays);
175}
176
177
178void HavranRayCaster::CastRays16(SimpleRayContainer &rays,
179                                 int offset,
180                                 VssRayContainer &vssRays,
181                                 const AxisAlignedBox3 &sbox,
182                                 const bool castDoubleRays,
183                                 const bool pruneInvalidRays)
184{
185#ifdef USE_HAVRAN_RAYCASTER
186 
187#if DEBUG_RAYCAST
188  Debug << "C16 " << flush;
189#endif
190
191#if 0
192  SimpleRayContainer::const_iterator sit = rays.begin() + offset;
193  SimpleRayContainer::const_iterator sit_end = sit + 16;
194 
195  // no acceleration for ray bundles implemented right now
196  for (sit = rays.begin(); sit != sit_end; ++ sit)
197  {
198    CastRay(*sit, vssRays, sbox, castDoubleRays, pruneInvalidRays);
199  }
200#else
201  if (castDoubleRays) {
202#if 0
203    // Special routine to cast double sided rays
204    mKtbtree->SetOffset(0);
205    mKtbtree->FindNearestI_16twoDir(rays);
206#else
207    // Here we decompose shooting into two phases
208#if 1   
209    // Here we shoot first forward rays and then backward ones
210    mKtbtree->SetOffset(0);
211    mKtbtree->FindNearestI_16oneDir(rays, offset);
212    SimpleRayContainer::iterator sit = rays.begin() + offset;
213    SimpleRayContainer::const_iterator sit_end = rays.begin() + offset + 16;
214    for (; sit != sit_end; ++ sit)
215    {
216      (*sit).mDirection = - (*sit).mDirection;
217    }
218    // We store the results at different place
219    mKtbtree->SetOffset(16);
220    mKtbtree->FindNearestI_16oneDir(rays, offset);
221    sit = rays.begin() + offset;
222    for ( ; sit != sit_end; ++ sit)
223    {
224      (*sit).mDirection = - (*sit).mDirection;
225    }
226#else
227    // Here we shoot first backward rays and forward ones
228    SimpleRayContainer::iterator sit = rays.begin() + offset;
229    SimpleRayContainer::const_iterator sit_end = rays.begin() + offset + 16;
230    for ( ; sit != sit_end; ++ sit)
231    {
232      (*sit).mDirection = - (*sit).mDirection;
233    }
234    // backward rays to be shot
235    mKtbtree->SetOffset(16);
236    mKtbtree->FindNearestI_16oneDir(rays, offset);
237    sit = rays.begin() + offset;
238    for ( ; sit != sit_end; ++ sit)
239    {
240      (*sit).mDirection = - (*sit).mDirection;
241    }
242    // forward rays to be shot
243    mKtbtree->SetOffset(0);
244    mKtbtree->FindNearestI_16oneDir(rays, offset);
245#endif   
246#endif   
247  }
248  else {
249    // Shoot all 16 rays  at the same time using a special algorithm
250    mKtbtree->SetOffset(0);
251    mKtbtree->FindNearestI_16oneDir(rays, offset);   
252  }
253#endif
254
255 
256  for (int i=0, k=offset; i < 16; i++, k++)
257  {
258    Intersection hitA(rays[k].mOrigin), hitB(rays[k].mOrigin);
259
260#if DEBUG_RAYCAST
261    Debug<<"FH\n"<<flush;
262#endif
263
264    Intersectable *intersect = SimpleRay::IntersectionRes[i].intersectable;
265
266    if (intersect)
267    {
268      hitA.mObject = intersect;
269      // Get the normal of that face
270      hitA.mNormal = intersect->GetNormal(0);
271     
272      //-rays[index+i].mDirection; // $$ temporary
273      float tdist = SimpleRay::IntersectionRes[i].tdist;
274      hitA.mPoint = rays[k].Extrap(tdist);
275    }
276   
277#if DEBUG_RAYCAST
278    Debug<<"BH\n"<<flush;
279#endif
280
281    if (castDoubleRays)
282    {
283      Intersectable *intersect =
284                SimpleRay::IntersectionRes[i+16].intersectable;
285
286      if (intersect)
287      {
288                hitB.mObject = intersect;
289                hitB.mNormal = intersect->GetNormal(0);
290                // $$ JB : i -> i+16
291                float tdist = SimpleRay::IntersectionRes[i+16].tdist;
292                hitB.mPoint = rays[k].Extrap(-tdist);
293      }
294    }
295   
296#if DEBUG_RAYCAST
297    Debug<<"PR\n"<<flush;
298#endif
299
300#if 1
301    ProcessRay(rays[k],
302                           hitA,
303                           hitB,
304                           vssRays,
305                           sbox,
306                           castDoubleRays,
307                           pruneInvalidRays
308                           );
309#endif
310  } // for
311
312 
313#if DEBUG_RAYCAST
314  Debug<<"C16F\n"<<flush;
315#endif
316
317#endif // USE_HAVRAN_RAYCASTER
318}
319
320
321
322void
323HavranRayCaster::CastSimpleForwardRays(
324                                       SimpleRayContainer &rays,
325                                       const AxisAlignedBox3 &sbox
326                                      )
327{
328  int hit_triangles[16];
329  float dist[16];
330  Vector3 normals[16];
331  Vector3 min = sbox.Min();
332  Vector3 max = sbox.Max();
333 
334  int packets = rays.size() / 16;
335 
336  int i, j, k = 0;
337  Vector3 dir;
338 
339  // By groups of rays
340  for (i=0; i < packets; i++) {
341    int offset = i * 16;
342    mKtbtree->FindNearestI_16oneDir(rays, offset);
343    // ??? What to do with the results ? These are
344    // not used at the moment also in IntelRayCaster.cpp
345  } // for
346
347
348  for (; k < rays.size(); k++) {
349    double normal[3];
350    mKtbtree->FindNearestI(rays[k]);
351    // ??? What to do with the results ? These are
352    // not used at the moment also in IntelRayCaster.cpp
353  }
354
355
356  return;
357}
358
359void HavranRayCaster::CastRays(
360                               SimpleRayContainer &rays,
361                               VssRayContainer &vssRays,
362                               const AxisAlignedBox3 &sbox,
363                               const bool castDoubleRay,
364                               const bool pruneInvalidRays )
365{
366  int buckets = rays.size()/16;
367  int offset = 0;
368
369#if 0
370  int time = GetTime();
371  CastSimpleForwardRays(rays, sbox);
372  cout<<1e-3*2*rays.size()/TimeDiff(time, GetTime())<<" Mrays/sec"<<endl;
373#endif
374 
375  for (int i=0; i < buckets; i++, offset+=16) {
376    CastRays16(rays, offset, vssRays, sbox,
377               castDoubleRay, pruneInvalidRays);
378
379    if ((int)rays.size() > 100000 && i % (100000/16) == 0)
380      cout<<"\r"<<offset<<"/"<<(int)rays.size()<<"\r";
381  }
382
383  for (; offset < (int)rays.size(); offset++)
384    CastRay(rays[offset], vssRays, sbox, castDoubleRay, pruneInvalidRays);
385}
386
387#if 0
388// BUG1 41579  196948 1064111
389// BUG2 254    1672   10869
390 
391// Just for testing concept
392void
393HavranRayCaster::CastRaysPacket2x2(RayPacket2x2 &raysPack,
394                                                                   bool castDoubleRay,
395                                                                   const bool pruneInvalidRays)
396{
397#ifdef USE_HAVRAN_RAYCASTER
398
399  if (castDoubleRay) {
400    // cast forward rays
401    mKtbtree->FindNearestI(raysPack);
402    for (int i = 0; i < 4; i++)
403      raysPack.SetDir(i, -raysPack.GetDir(i));
404    // cast backward rays
405    mKtbtree->FindNearestI(raysPack);
406    // reverse the rays back
407    for (int i = 0; i < 4; i++)
408      raysPack.SetDir(i, -raysPack.GetDir(i));
409  }
410  else {
411    // only rays forwards
412    mKtbtree->FindNearestI(raysPack);
413#if 0
414    // Only verification of correctness by casting single rays
415    static int cntBugs = 0;
416    SimpleRay ray;
417    int cntErrors = 0, cntDistErrors = 0;
418    bool newBug = false;
419    for (int i = 0; i < 4; i++) {
420      ray.mOrigin = raysPack.GetLoc(i);
421      ray.mDirection = raysPack.GetDir(i);
422      mKtbtree->FindNearestI(ray);
423      if (raysPack.GetObject(i) != SimpleRay::IntersectionRes[0].intersectable) {
424        float dist = (raysPack.GetT(i) - SimpleRay::IntersectionRes[0].tdist);
425        if (fabs(dist) > 0.001f) {
426          cntErrors++; newBug = true;
427          cntBugs++;
428          cout << " BUG1 d= " << dist;
429        }
430      }
431      else {
432        float dist = 0.f;
433        if (raysPack.GetObject(i) && SimpleRay::IntersectionRes[0].intersectable)
434          if (fabs((dist=(fabs (raysPack.GetT(i) - SimpleRay::IntersectionRes[0].tdist)))) > 1.f) {
435            cntDistErrors++; newBug = true; cntBugs++;
436            cout << " BUG2 distdf= " << dist ;     
437          }
438      }
439    } // for
440    if (newBug) cout << " CB= " << cntBugs << "\n";
441#endif
442  }
443
444  return;
445#endif // USE_HAVRAN_RAYCASTER
446}
447#endif
448
449
450} // the namespace
Note: See TracBrowser for help on using the repository browser.