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

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