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

Revision 2606, 11.4 KB checked in by bittner, 16 years ago (diff)

merge on nemo
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  mlrtaTraverseGroupASEye4(&minBox.x, &maxBox.x, result4, dist4);
83
84  if (saveMutationRays) {
85    fprintf(fileOut, "I %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f\n",
86                        minBox.x, minBox.y, minBox.z, maxBox.x, maxBox.y, maxBox.z);
87    for (int i = 0; i < 4; i++) {
88      fprintf(fileOut, "%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
89                          cntSavedRays,
90                          origin4[i].x,
91                          origin4[i].y,
92                          origin4[i].z,
93                          dirs4[i].x,
94                          dirs4[i].y,
95                          dirs4[i].z,
96                          (result4[i] != -1) ? 1 : 0,
97                          (result4[i] != -1) ? dist4[i] : 0);
98      cntSavedRays ++;
99      cntSavedRaysFLUSH++;
100    }
101    if (cntSavedRaysFLUSH > intSavedLIMIT) {
102      fflush(fileOut);
103      cntSavedRaysFLUSH = 0;
104    }
105  } // if we are saving rays
106
107  return;
108}
109
110
111 
112int IntelRayCaster::CastRay(const SimpleRay &simpleRay,
113                                                        VssRayContainer &vssRays,
114                                                        const AxisAlignedBox3 &box,
115                                                        const bool castDoubleRay,
116                                                        const bool pruneInvalidRays
117                                                        )
118{ 
119  //cout << "intel ray" << endl;
120  VssRay *vssRay  = NULL;
121  int hits = 0;
122  int hittriangle;
123  Intersection hitA(simpleRay.mOrigin), hitB(simpleRay.mOrigin);
124 
125  float dist;
126  float dist1, dist2;
127  double normal[3];
128 
129  hittriangle = mlrtaIntersectAS(
130                                 &simpleRay.mOrigin.x,
131                                 &simpleRay.mDirection.x,
132                                 normal,
133                                 dist);
134  dist1 = dist;
135 
136  if (hittriangle != -1 ) {
137    Intersectable *intersect = mPreprocessor.GetParentObject(hittriangle);
138   
139    if (intersect)
140    {
141      hitA.mObject = intersect;
142      hitA.mNormal = Vector3(normal[0], normal[1], normal[2]);
143      // Get the normal of that face
144      //                Mesh *mesh = ((MeshInstance *)objectA)->GetMesh();
145      //                normalA = mesh->GetFacePlane(mFaceParents[forward_hit_triangles[i]].mFaceIndex).mNormal;
146      //-rays[index+i].mDirection; // $$ temporary
147      hitA.mPoint = simpleRay.Extrap(dist);
148    }
149  }
150 
151  if (castDoubleRay)
152  {
153    Vector3 dir = -simpleRay.mDirection;
154    hittriangle = mlrtaIntersectAS(
155                                   &simpleRay.mOrigin.x,
156                                   &dir.x,
157                                   normal,
158                                   dist);
159    dist2 = dist;
160
161    Intersectable *intersect = mPreprocessor.GetParentObject(hittriangle);
162               
163    if (intersect)
164    {
165      hitB.mObject = intersect;
166      hitB.mNormal = Vector3(normal[0], normal[1], normal[2]);
167      // Get the normal of that face
168      //                Mesh *mesh = ((MeshInstance *)objectB)->GetMesh();
169      //                normalA = mesh->GetFacePlane(mFaceParents[forward_hit_triangles[i]].mFaceIndex).mNormal;
170      //-rays[index+i].mDirection; // $$ temporary
171      hitB.mPoint = simpleRay.Extrap(dist);
172    }
173  }
174
175  if (saveRays && preprocessor->mTotalRaysCast > saveRaysStart) {   
176    if (castDoubleRay)
177        {
178      fprintf(fileOut, "D %d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f %d %4.7f\n",
179              cntSavedRays,
180              simpleRay.mOrigin.x,
181              simpleRay.mOrigin.y,
182              simpleRay.mOrigin.z,
183              simpleRay.mDirection.x,
184              simpleRay.mDirection.y,
185              simpleRay.mDirection.z,
186              hitA.mObject ? 1 : 0,
187              hitA.mObject ? dist1 : 0,
188              hitB.mObject ? 1 : 0,
189              hitB.mObject ? dist2 : 0
190              );
191        }
192    else
193      fprintf(fileOut, "S %d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
194              cntSavedRays,
195              simpleRay.mOrigin.x,
196              simpleRay.mOrigin.y,
197              simpleRay.mOrigin.z,
198              simpleRay.mDirection.x,
199              simpleRay.mDirection.y,
200              simpleRay.mDirection.z,
201              hitA.mObject ? 1 : 0,
202              hitA.mObject ? dist1 : 0
203              );
204    cntSavedRays++;
205    cntSavedRaysFLUSH++;
206    if (cntSavedRaysFLUSH > intSavedLIMIT) {
207      fflush(fileOut);
208      cntSavedRaysFLUSH = 0;
209    }
210  }     
211 
212  return ProcessRay(
213                    simpleRay,
214                    hitA,
215                    hitB,
216                    vssRays,
217                    box,
218                    castDoubleRay,
219                    pruneInvalidRays
220                    );
221}
222 
223
224void IntelRayCaster::CastRays16(
225                                SimpleRayContainer &rays,
226                                VssRayContainer &vssRays,
227                                const AxisAlignedBox3 &sbox,
228                                const bool castDoubleRay,
229                                const bool pruneInvalidRays)
230{
231  CastRays16(rays, 0, vssRays, sbox, castDoubleRay, pruneInvalidRays);
232}
233 
234void IntelRayCaster::CastRays16(
235                                SimpleRayContainer &rays,
236                                const int offset,
237                                VssRayContainer &vssRays,
238                                const AxisAlignedBox3 &sbox,
239                                const bool castDoubleRay,
240                                const bool pruneInvalidRays)
241{ 
242  int i, k;
243  const int num = 16;
244 
245#if DEBUG_RAYCAST
246  Debug<<"C16 "<<flush;
247  static int counter=0;
248  Debug<<counter++<<endl;
249#endif
250 
251  static int forward_hit_triangles[16];
252  static float forward_dist[16];
253 
254  static int backward_hit_triangles[16];
255  static float backward_dist[16];
256 
257 
258  Vector3 min = mPreprocessor.mSceneGraph->GetBox().Min();
259  Vector3 max = mPreprocessor.mSceneGraph->GetBox().Max();
260
261  for (k=offset, i=0; i < num; i++, k++) {
262#if DEBUG_RAYCAST
263    if (counter == 43964) {
264      Debug<<rays[k].mOrigin<<" "<<rays[k].mDirection<<endl;
265    }
266#endif
267    mlrtaStoreRayAS16(&rays[k].mOrigin.x,
268                      &rays[k].mDirection.x,
269                      i);
270  }
271 
272#if DEBUG_RAYCAST
273  Debug<<"TA\n"<<flush;
274#endif
275 
276  mlrtaTraverseGroupAS16(&min.x,
277                         &max.x,
278                         forward_hit_triangles,
279                         forward_dist);
280 
281#if DEBUG_RAYCAST
282  Debug<<"TAB\n"<<flush;
283#endif
284
285  if (castDoubleRay) {
286    for (k=offset, i=0; i < num; i++, k++)  {
287      Vector3 dir = -rays[k].mDirection;
288      mlrtaStoreRayAS16(&rays[k].mOrigin.x,
289                        &dir.x,
290                        i);
291    }
292   
293#if DEBUG_RAYCAST
294    Debug<<"TB\n"<<flush;
295#endif
296   
297    mlrtaTraverseGroupAS16(&min.x,
298                           &max.x,
299                           backward_hit_triangles,
300                           backward_dist);
301  }
302       
303#if DEBUG_RAYCAST
304  Debug<<"BBB\n"<<flush;
305#endif
306
307  if (saveRays && preprocessor->mTotalRaysCast > saveRaysStart) {
308    if (castDoubleRay)
309      fprintf(fileOut, "G\n");
310    else
311      fprintf(fileOut, "H\n");
312  }
313 
314  for (i=0, k=offset; i < num; i++, k++)
315  {
316    Intersection hitA(rays[k].mOrigin), hitB(rays[k].mOrigin);
317
318#if DEBUG_RAYCAST
319    Debug<<"FH\n"<<flush;
320#endif
321
322    Intersectable *intersect =
323      mPreprocessor.GetParentObject(forward_hit_triangles[i]);
324
325    if (intersect)
326    {
327      hitA.mObject = intersect;
328      // Get the normal of that face
329      hitA.mNormal = mPreprocessor.GetParentNormal(forward_hit_triangles[i]);
330     
331      //-rays[index+i].mDirection; // $$ temporary
332      hitA.mPoint = rays[k].Extrap(forward_dist[i]);
333    }
334   
335#if DEBUG_RAYCAST
336    Debug<<"BH\n"<<flush;
337#endif
338
339    if (castDoubleRay)
340    {
341      Intersectable *intersect =
342        mPreprocessor.GetParentObject(backward_hit_triangles[i]);
343     
344      if (intersect)
345      {
346        hitB.mObject = intersect;
347        hitB.mNormal = mPreprocessor.GetParentNormal(backward_hit_triangles[i]);
348       
349        // normalB = rays[i].mDirection; // $$ temporary
350        hitB.mPoint = rays[k].Extrap(-backward_dist[i]);
351      }
352    }
353
354
355    if (saveRays && preprocessor->mTotalRaysCast > saveRaysStart) {
356      if (castDoubleRay)     
357        fprintf(fileOut, "%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f %d %4.7f\n",
358                cntSavedRays,
359                rays[k].mOrigin.x,
360                rays[k].mOrigin.y,
361                rays[k].mOrigin.z,
362                rays[k].mDirection.x,
363                rays[k].mDirection.y,
364                rays[k].mDirection.z,
365                hitA.mObject ? 1 : 0,
366                hitA.mObject ? forward_dist[i] : 0,
367                hitB.mObject ? 1 : 0,
368                hitB.mObject ? -backward_dist[i] : 0
369                );
370      else
371        fprintf(fileOut, "%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
372                cntSavedRays,
373                rays[k].mOrigin.x,
374                rays[k].mOrigin.y,
375                rays[k].mOrigin.z,
376                rays[k].mDirection.x,
377                rays[k].mDirection.y,
378                rays[k].mDirection.z,
379                hitA.mObject ? 1 : 0,
380                hitA.mObject ? forward_dist[i] : 0
381                );
382      cntSavedRays++;
383      cntSavedRaysFLUSH++;
384    }
385
386   
387#if DEBUG_RAYCAST
388    Debug<<"PR\n"<<flush;
389#endif
390
391#if 1
392    ProcessRay(rays[k],
393               hitA,
394               hitB,
395               vssRays,
396               sbox,
397               castDoubleRay,
398               pruneInvalidRays
399               );
400#endif
401  } // for
402
403  if (saveRays) {
404    if (cntSavedRaysFLUSH > intSavedLIMIT) {
405      fflush(fileOut);
406      cntSavedRaysFLUSH = 0;
407    }
408  }
409
410#if DEBUG_RAYCAST
411  Debug<<"C16F\n"<<flush;
412#endif
413}
414
415
416
417
418
419void
420IntelRayCaster::CastSimpleForwardRays(
421                                      SimpleRayContainer &rays,
422                                      const AxisAlignedBox3 &sbox
423                                      )
424{
425  int hit_triangles[16];
426  float dist[16];
427  Vector3 normals[16];
428  Vector3 min = sbox.Min();
429  Vector3 max = sbox.Max();
430 
431  int packets = rays.size() / 16;
432 
433  int i, j, k = 0;
434  Vector3 dir;
435 
436  for (i=0; i < packets; i++) {
437    for (j=0; j < 16; j++, k++)
438      mlrtaStoreRayAS16(&rays[k].mOrigin.x,
439                        &rays[k].mDirection.x,
440                        j);
441   
442    mlrtaTraverseGroupAS16(&min.x,
443                           &max.x,
444                           hit_triangles,
445                           dist);       
446  } // for
447
448
449  for (; k < rays.size(); k++) {
450        double normal[3];
451        hit_triangles[0] = mlrtaIntersectAS(
452                                            &rays[k].mOrigin.x,
453                                            &rays[k].mDirection.x,
454                                            normal,
455                                            dist[0]);
456  }
457 
458}
459
460
461void
462IntelRayCaster::CastRays(
463                         SimpleRayContainer &rays,
464                         VssRayContainer &vssRays,
465                         const AxisAlignedBox3 &sbox,
466                         const bool castDoubleRay,
467                         const bool pruneInvalidRays )
468{
469
470  int buckets = rays.size()/16;
471  int offset = 0;
472
473#if 0
474  int time = GetTime();
475  CastSimpleForwardRays(rays, sbox);
476  cout<<1e-3*2*rays.size()/TimeDiff(time, GetTime())<<" Mrays/sec"<<endl;
477#endif
478 
479  for (int i=0; i < buckets; i++, offset+=16) {
480        CastRays16(rays, offset, vssRays, sbox, castDoubleRay, pruneInvalidRays);
481
482        if ((int)rays.size() > 100000 && i % (100000/16) == 0)
483          cout<<"\r"<<offset<<"/"<<(int)rays.size()<<"\r";
484  }
485
486  for (; offset < (int)rays.size(); offset++)
487        CastRay(rays[offset], vssRays, sbox, castDoubleRay, pruneInvalidRays);
488
489}
490
491}
492
493#endif
Note: See TracBrowser for help on using the repository browser.