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

Revision 2635, 13.6 KB checked in by bittner, 17 years ago (diff)

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