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

Revision 2592, 13.3 KB checked in by bittner, 16 years ago (diff)

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