source: GTP/trunk/Lib/Vis/Preprocessing/src/IntelRayCaster.cpp @ 2743

Revision 2743, 11.6 KB checked in by mattausch, 16 years ago (diff)
Line 
1#include "IntelRayCaster.h"
2#include "VssRay.h"
3#include "Preprocessor.h"
4#include "SceneGraph.h"
5
6#ifdef GTP_INTERNAL
7#include "ArchModeler2MLRT.hxx"
8
9#define DEBUG_RAYCAST 0
10
11
12namespace GtpVisibilityPreprocessor {
13
14
15FILE *fileOut = 0;
16bool saveRays = false;
17bool saveMutationRays = false;
18const int saveRaysStart = 3000000;
19int cntSavedRaysFLUSH = 0;
20unsigned long int cntSavedRays = 0;
21const int intSavedLIMIT = 1024;
22
23void
24InitSaving()
25{
26  fileOut = fopen("fileRays.txt", "wb");
27  if (!fileOut) {
28    cout << "ERROR: file fileRays.txt cannot be opened .. exiting" << endl;
29    exit(3);
30  }
31}
32void
33FinishSaving()
34{
35  fclose(fileOut);
36}
37
38IntelRayCaster::IntelRayCaster(const Preprocessor &preprocessor,
39                               const string externKdTree):
40RayCaster(preprocessor)
41{
42  if (!InitRayCast(externKdTree))
43    cout << "warning: intel ray tracer could not be initialized!" << endl;
44  if (saveRays || saveMutationRays)
45    InitSaving();
46}
47
48
49IntelRayCaster::~IntelRayCaster()
50{
51  if (saveRays)
52    FinishSaving();
53}
54
55
56bool IntelRayCaster::InitRayCast(const string externKdTree)
57{
58  cout<<"Intel ray cast file: " << externKdTree << endl;
59       
60  return mlrtaLoadAS(externKdTree.c_str());
61}
62
63
64// Using packet of 4 rays supposing that these are coherent
65// We give a box to which each ray is clipped to before the
66// ray shooting is computed !
67
68// Using packet of 4 rays supposing that these are coherent
69// We give a box to which each ray is clipped to before the
70// ray shooting is computed !
71void IntelRayCaster::CastRaysPacket4(const Vector3 &minBox,
72                                                                         const Vector3 &maxBox,
73                                                                         const Vector3 origin4[],
74                                                                         const Vector3 dirs4[],
75                                                                         int     result4[],
76                                                                         float   dist4[])
77{
78  for (int i = 0; i < 4; i++) {
79    mlrtaStoreRayASEye4(&origin4[i].x, &dirs4[i].x, i);
80  }
81 
82  rawCastTimer.Entry();
83  mlrtaTraverseGroupASEye4(&minBox.x, &maxBox.x, result4, dist4);
84  rawCastTimer.Exit();
85
86  if (saveMutationRays) {
87    fprintf(fileOut, "I %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f\n",
88                        minBox.x, minBox.y, minBox.z, maxBox.x, maxBox.y, maxBox.z);
89    for (int i = 0; i < 4; i++) {
90      fprintf(fileOut, "%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
91                          cntSavedRays,
92                          origin4[i].x,
93                          origin4[i].y,
94                          origin4[i].z,
95                          dirs4[i].x,
96                          dirs4[i].y,
97                          dirs4[i].z,
98                          (result4[i] != -1) ? 1 : 0,
99                          (result4[i] != -1) ? dist4[i] : 0);
100      cntSavedRays ++;
101      cntSavedRaysFLUSH++;
102    }
103    if (cntSavedRaysFLUSH > intSavedLIMIT) {
104      fflush(fileOut);
105      cntSavedRaysFLUSH = 0;
106    }
107  } // if we are saving rays
108
109  return;
110}
111
112
113 
114int IntelRayCaster::CastRay(const SimpleRay &simpleRay,
115                                                        VssRayContainer &vssRays,
116                                                        const AxisAlignedBox3 &box,
117                                                        const bool castDoubleRay,
118                                                        const bool pruneInvalidRays
119                                                        )
120{ 
121  //cout << "intel ray" << endl;
122  VssRay *vssRay = NULL;
123  int hits = 0;
124  int hittriangle;
125  Intersection hitA(simpleRay.mOrigin), hitB(simpleRay.mOrigin);
126 
127  float dist;
128  double normal[3];
129
130  rawCastTimer.Entry();
131 
132  hittriangle = mlrtaIntersectAS(
133                                 &simpleRay.mOrigin.x,
134                                 &simpleRay.mDirection.x,
135                                 normal,
136                                 dist);
137
138  rawCastTimer.Exit();
139 
140 
141  if (hittriangle != -1 ) {
142    Intersectable *intersect = mPreprocessor.GetParentObject(hittriangle);
143   
144    if (intersect)
145    {
146                hitA.mObject = intersect;
147                hitA.mNormal = Vector3(normal[0], normal[1], normal[2]);
148
149                // Get the normal of that face
150                // Mesh *mesh = ((MeshInstance *)objectA)->GetMesh();
151                // normalA = mesh->GetFacePlane(mFaceParents[forward_hit_triangles[i]].mFaceIndex).mNormal;
152                // -rays[index+i].mDirection; // $$ temporary
153                hitA.mPoint = simpleRay.Extrap(dist);
154    }
155  }
156 
157  if (castDoubleRay)
158  {
159    Vector3 dir = -simpleRay.mDirection;
160
161        rawCastTimer.Entry();
162
163        hittriangle = mlrtaIntersectAS(
164                                   &simpleRay.mOrigin.x,
165                                   &dir.x,
166                                   normal,
167                                   dist);
168       
169        rawCastTimer.Exit();
170
171        Intersectable *intersect = mPreprocessor.GetParentObject(hittriangle);
172               
173    if (intersect)
174    {
175                hitB.mObject = intersect;
176                hitB.mNormal = Vector3(normal[0], normal[1], normal[2]);
177                // Get the normal of that face
178                // Mesh *mesh = ((MeshInstance *)objectB)->GetMesh();
179        //  normalA = mesh->GetFacePlane(mFaceParents[forward_hit_triangles[i]].mFaceIndex).mNormal;
180                //-rays[index+i].mDirection; // $$ temporary
181                hitB.mPoint = simpleRay.Extrap(dist);
182    }
183  }
184/*
185  if (saveRays && preprocessor->mTotalRaysCast > saveRaysStart) {   
186    if (castDoubleRay)
187        {
188      fprintf(fileOut, "D %d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f %d %4.7f\n",
189              cntSavedRays,
190              simpleRay.mOrigin.x,
191              simpleRay.mOrigin.y,
192              simpleRay.mOrigin.z,
193              simpleRay.mDirection.x,
194              simpleRay.mDirection.y,
195              simpleRay.mDirection.z,
196              hitA.mObject ? 1 : 0,
197              hitA.mObject ? dist1 : 0,
198              hitB.mObject ? 1 : 0,
199              hitB.mObject ? dist2 : 0
200              );
201        }
202    else
203      fprintf(fileOut, "S %d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
204              cntSavedRays,
205              simpleRay.mOrigin.x,
206              simpleRay.mOrigin.y,
207              simpleRay.mOrigin.z,
208              simpleRay.mDirection.x,
209              simpleRay.mDirection.y,
210              simpleRay.mDirection.z,
211              hitA.mObject ? 1 : 0,
212              hitA.mObject ? dist1 : 0
213              );
214    cntSavedRays++;
215    cntSavedRaysFLUSH++;
216    if (cntSavedRaysFLUSH > intSavedLIMIT) {
217      fflush(fileOut);
218      cntSavedRaysFLUSH = 0;
219    }
220  }     
221  */
222  return ProcessRay(
223                    simpleRay,
224                    hitA,
225                    hitB,
226                    vssRays,
227                    box,
228                    castDoubleRay,
229                    pruneInvalidRays
230                    );
231}
232 
233
234void IntelRayCaster::CastRays16(
235                                SimpleRayContainer &rays,
236                                VssRayContainer &vssRays,
237                                const AxisAlignedBox3 &sbox,
238                                const bool castDoubleRay,
239                                const bool pruneInvalidRays)
240{
241  CastRays16(rays, 0, vssRays, sbox, castDoubleRay, pruneInvalidRays);
242}
243 
244void IntelRayCaster::CastRays16(
245                                SimpleRayContainer &rays,
246                                const int offset,
247                                VssRayContainer &vssRays,
248                                const AxisAlignedBox3 &sbox,
249                                const bool castDoubleRay,
250                                const bool pruneInvalidRays)
251{ 
252  int i, k;
253  const int num = 16;
254 
255#if DEBUG_RAYCAST
256  Debug<<"C16 "<<flush;
257  static int counter=0;
258  Debug<<counter++<<endl;
259#endif
260 
261  static int forward_hit_triangles[16];
262  static float forward_dist[16];
263 
264  static int backward_hit_triangles[16];
265  static float backward_dist[16];
266 
267 
268  Vector3 min = mPreprocessor.mSceneGraph->GetBox().Min();
269  Vector3 max = mPreprocessor.mSceneGraph->GetBox().Max();
270
271  for (k=offset, i=0; i < num; i++, k++) {
272#if DEBUG_RAYCAST
273    if (counter == 43964) {
274      Debug<<rays[k].mOrigin<<" "<<rays[k].mDirection<<endl;
275    }
276#endif
277
278    mlrtaStoreRayAS16(&rays[k].mOrigin.x,
279                      &rays[k].mDirection.x,
280                      i);
281  }
282 
283#if DEBUG_RAYCAST
284  Debug<<"TA\n"<<flush;
285#endif
286  rawCastTimer.Entry();
287
288  mlrtaTraverseGroupAS16(&min.x,
289                         &max.x,
290                         forward_hit_triangles,
291                         forward_dist);
292 
293  rawCastTimer.Exit();
294 
295#if DEBUG_RAYCAST
296  Debug<<"TAB\n"<<flush;
297#endif
298
299  if (castDoubleRay) {
300    for (k=offset, i=0; i < num; i++, k++)  {
301      Vector3 dir = -rays[k].mDirection;
302      mlrtaStoreRayAS16(&rays[k].mOrigin.x,
303                        &dir.x,
304                        i);
305    }
306   
307#if DEBUG_RAYCAST
308    Debug<<"TB\n"<<flush;
309#endif
310    rawCastTimer.Entry();
311    mlrtaTraverseGroupAS16(&min.x,
312                           &max.x,
313                           backward_hit_triangles,
314                           backward_dist);
315
316        rawCastTimer.Exit();
317  }
318       
319#if DEBUG_RAYCAST
320  Debug<<"BBB\n"<<flush;
321#endif
322
323  if (saveRays && preprocessor->mTotalRaysCast > saveRaysStart) {
324    if (castDoubleRay)
325      fprintf(fileOut, "G\n");
326    else
327      fprintf(fileOut, "H\n");
328  }
329 
330  for (i=0, k=offset; i < num; i++, k++)
331  {
332    Intersection hitA(rays[k].mOrigin), hitB(rays[k].mOrigin);
333
334#if DEBUG_RAYCAST
335    Debug<<"FH\n"<<flush;
336#endif
337
338    Intersectable *intersect =
339      mPreprocessor.GetParentObject(forward_hit_triangles[i]);
340
341    if (intersect)
342    {
343      hitA.mObject = intersect;
344      // Get the normal of that face
345      hitA.mNormal = mPreprocessor.GetParentNormal(forward_hit_triangles[i]);
346     
347      //-rays[index+i].mDirection; // $$ temporary
348      hitA.mPoint = rays[k].Extrap(forward_dist[i]);
349    }
350   
351#if DEBUG_RAYCAST
352    Debug<<"BH\n"<<flush;
353#endif
354
355    if (castDoubleRay)
356    {
357      Intersectable *intersect =
358                  mPreprocessor.GetParentObject(backward_hit_triangles[i]);
359     
360      if (intersect)
361      {
362                  hitB.mObject = intersect;
363                  hitB.mNormal = mPreprocessor.GetParentNormal(backward_hit_triangles[i]);
364                 
365                  // normalB = rays[i].mDirection; // $$ temporary
366                  hitB.mPoint = rays[k].Extrap(-backward_dist[i]);
367      }
368    }
369
370
371    if (saveRays && preprocessor->mTotalRaysCast > saveRaysStart) {
372      if (castDoubleRay)     
373        fprintf(fileOut, "%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f %d %4.7f\n",
374                cntSavedRays,
375                rays[k].mOrigin.x,
376                rays[k].mOrigin.y,
377                rays[k].mOrigin.z,
378                rays[k].mDirection.x,
379                rays[k].mDirection.y,
380                rays[k].mDirection.z,
381                hitA.mObject ? 1 : 0,
382                hitA.mObject ? forward_dist[i] : 0,
383                hitB.mObject ? 1 : 0,
384                hitB.mObject ? -backward_dist[i] : 0
385                );
386      else
387        fprintf(fileOut, "%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
388                cntSavedRays,
389                rays[k].mOrigin.x,
390                rays[k].mOrigin.y,
391                rays[k].mOrigin.z,
392                rays[k].mDirection.x,
393                rays[k].mDirection.y,
394                rays[k].mDirection.z,
395                hitA.mObject ? 1 : 0,
396                hitA.mObject ? forward_dist[i] : 0
397                );
398      cntSavedRays++;
399      cntSavedRaysFLUSH++;
400    }
401
402   
403#if DEBUG_RAYCAST
404    Debug<<"PR\n"<<flush;
405#endif
406
407#if 1
408    ProcessRay(rays[k],
409               hitA,
410               hitB,
411               vssRays,
412               sbox,
413               castDoubleRay,
414               pruneInvalidRays
415               );
416#endif
417  } // for
418
419  if (saveRays) {
420    if (cntSavedRaysFLUSH > intSavedLIMIT) {
421      fflush(fileOut);
422      cntSavedRaysFLUSH = 0;
423    }
424  }
425
426#if DEBUG_RAYCAST
427  Debug<<"C16F\n"<<flush;
428#endif
429}
430
431
432
433
434
435void
436IntelRayCaster::CastSimpleForwardRays(
437                                      SimpleRayContainer &rays,
438                                      const AxisAlignedBox3 &sbox
439                                      )
440{
441  int hit_triangles[16];
442  float dist[16];
443  Vector3 normals[16];
444  Vector3 min = sbox.Min();
445  Vector3 max = sbox.Max();
446 
447  int packets = (int)rays.size() / 16;
448 
449  int i, j, k = 0;
450  Vector3 dir;
451 
452  for (i=0; i < packets; i++) {
453    for (j=0; j < 16; j++, k++)
454      mlrtaStoreRayAS16(&rays[k].mOrigin.x,
455                        &rays[k].mDirection.x,
456                        j);
457    rawCastTimer.Entry();
458
459    mlrtaTraverseGroupAS16(&min.x,
460                           &max.x,
461                           hit_triangles,
462                           dist);       
463        rawCastTimer.Exit();
464  } // for
465
466
467  for (; k < rays.size(); k++) {
468        double normal[3];
469        rawCastTimer.Entry();
470        hit_triangles[0] = mlrtaIntersectAS(
471                                            &rays[k].mOrigin.x,
472                                            &rays[k].mDirection.x,
473                                            normal,
474                                            dist[0]);
475        rawCastTimer.Exit();
476  }
477 
478}
479
480
481void
482IntelRayCaster::CastRays(SimpleRayContainer &rays,
483                                                 VssRayContainer &vssRays,
484                                                 const AxisAlignedBox3 &sbox,
485                                                 bool castDoubleRay,
486                                                 bool pruneInvalidRays)
487{
488
489  int buckets = (int)rays.size() / 16;
490  int offset = 0;
491
492#if 0
493  int time = GetTime();
494  CastSimpleForwardRays(rays, sbox);
495  cout<<1e-3*2*rays.size()/TimeDiff(time, GetTime())<<" Mrays/sec"<<endl;
496#endif
497 
498  for (int i=0; i < buckets; i++, offset+=16) {
499        CastRays16(rays, offset, vssRays, sbox, castDoubleRay, pruneInvalidRays);
500
501        if ((int)rays.size() > 100000 && i % (100000/16) == 0)
502          cout<<"\r"<<offset<<"/"<<(int)rays.size()<<"\r";
503  }
504
505  for (; offset < (int)rays.size(); offset++)
506        CastRay(rays[offset], vssRays, sbox, castDoubleRay, pruneInvalidRays);
507
508}
509
510}
511
512#endif
Note: See TracBrowser for help on using the repository browser.