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

Revision 2603, 13.7 KB checked in by bittner, 16 years ago (diff)

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