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

Revision 2582, 10.5 KB checked in by bittner, 16 years ago (diff)

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