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

Revision 2662, 13.8 KB checked in by mattausch, 16 years ago (diff)

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