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

Revision 2602, 13.4 KB checked in by bittner, 17 years ago (diff)

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