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

Revision 2705, 13.8 KB checked in by mattausch, 17 years ago (diff)

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