source: GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h @ 2629

Revision 2629, 13.4 KB checked in by bittner, 17 years ago (diff)

commit after merge with vlastimil

Line 
1// ============================================================================
2// $Id: raypack.h $
3//
4// raypack.h
5//        CRayPacket class - core of the ray-packet traversal routines
6//
7// Class: CRayPacket2x2
8//
9// REPLACEMENT_STRING
10//
11// Initial coding by Vlastimil Havran, 2006. The data design is in fact
12// Jakko Biker layout as proposed in the article on Intel Web Site in year 2005
13// http://www.intel.com/cd/ids/developer/asmo-na/eng/245711.htm?page=1
14
15#ifndef __RAYPACK_H__
16#define __RAYPACK_H__
17
18#include <cassert>
19
20namespace GtpVisibilityPreprocessor {
21
22#include "Vector3.h"
23#include "Matrix4x4.h"
24#include "ktbconf.h"
25 
26#ifdef _USE_HAVRAN_SSE
27 
28#ifdef __SSE__
29
30// System headers for SSE
31#ifdef __INTEL_COMPILER
32#include <xmmintrin.h>
33#else
34// We assume GNU GCC compiler 3.4 or higher
35#include <xmmintrin.h>
36#endif
37
38
39// forward declarations
40#define SSE_INTRINSIC
41#ifdef SSE_INTRINSIC
42#define ALIGN16             __declspec(align(16))
43
44#ifdef _MSC_VER
45#define NEW_ALIGN16(type,n)  ((type*)_aligned_malloc((n)*sizeof(type),16))
46#define FREE_ALIGN16(array)  if(array){_aligned_free(array);(array)=0;}
47#else
48#define NEW_ALIGN16(type,n)  ((type*)malloc((n)*sizeof(type)))
49#define FREE_ALIGN16(array)  if (array) { free(array);(array)=0;}
50#endif // _MSC_VC
51
52#define PAD_FOUR(h)           mulFour(h)
53union extract_m128
54{
55  __m128 m;
56  float f[4];
57};
58#else
59#define NEW_ALIGN16(type,n)  ((type*)malloc((n)*sizeof(type)))
60#define ALIGN16             
61#define FREE_ALIGN16(array)  if(array){free(array);(array)=0;}
62#define PAD_FOUR(h)           (h)
63#endif
64 
65// -------------------------------------------------------------------
66// RayPacket2x2 class.  A set of 4 ray is defined by a location and a
67// direction. The direction is always normalized (length == 1) during
68// use.
69// -------------------------------------------------------------------
70class RayPacket2x2
71{
72public:
73  enum {
74    // The number of rays in one packet
75    PACKSIZE = 4
76  };
77 
78  // constructor
79  RayPacket2x2(
80    // origin and the direction of rays
81    const Vector3 orf[],
82    const Vector3 dirf[],
83    // The same type for all four rays
84    const int _type,
85    // All the rays has to start at the same cell
86    const void *_originCell = NULL,
87    // All the rays has to start at the same object
88    const Intersectable *_startObject = NULL,
89    // and also for shadow rays finish at the same object
90    const Intersectable *_stopObject = NULL
91  ) {
92    // location
93    ox[0] = orf[0].x; ox[1] = orf[1].x; ox[2] = orf[2].x; ox[3] = orf[3].x;
94    oy[0] = orf[0].y; oy[1] = orf[1].y; oy[2] = orf[2].y; oy[3] = orf[3].y;
95    oz[0] = orf[0].z; oz[1] = orf[1].z; oz[2] = orf[2].z; oz[3] = orf[3].z;
96    // direction, we assume to be normalized
97    dx[0] = dirf[0].x; dx[1] = dirf[1].x; dx[2] = dirf[2].x; dx[3] = dirf[3].x;
98    dy[0] = dirf[0].y; dy[1] = dirf[1].y; dy[2] = dirf[2].y; dy[3] = dirf[3].y;
99    dz[0] = dirf[0].z; dz[1] = dirf[1].z; dz[2] = dirf[2].z; dz[3] = dirf[3].z;
100    // Other components
101    ttype = _type;
102    origin = _originCell;
103    termination = NULL;
104    startObject = _startObject;
105    stopObject = _stopObject;
106    Init();
107  }
108  // dummy constructor
109  RayPacket2x2() {}
110
111  // Inititalize the ray again when already constructed
112  void Init(
113    // origin and the direction of rays
114    const Vector3 orf[],
115    const Vector3 dirf[],
116    // The same type for all four rays
117    const int _type,
118    // All the rays has to start at the same cell
119    const void *_originCell = NULL,
120    // All the rays has to start at the same object
121    const Intersectable *_startObject = NULL,
122    // and also for shadow rays finish at the same object
123    const Intersectable *_stopObject = NULL,
124    // if the direction of rays is normalized or not
125    bool dirNormalized = false)
126  {
127    // location
128    ox[0] = orf[0].x; ox[1] = orf[1].x; ox[2] = orf[2].x; ox[3] = orf[3].x;
129    oy[0] = orf[0].y; oy[1] = orf[1].y; oy[2] = orf[2].y; oy[3] = orf[3].y;
130    oz[0] = orf[0].z; oz[1] = orf[1].z; oz[2] = orf[2].z; oz[3] = orf[3].z;
131    // direction
132    if (dirNormalized) {
133      // already normalized
134      dx[0] = dirf[0].x; dx[1] = dirf[1].x; dx[2] = dirf[2].x; dx[3] = dirf[3].x;
135      dy[0] = dirf[0].y; dy[1] = dirf[1].y; dy[2] = dirf[2].y; dy[3] = dirf[3].y;
136      dz[0] = dirf[0].z; dz[1] = dirf[1].z; dz[2] = dirf[2].z; dz[3] = dirf[3].z;
137    }
138    else {
139      dx[0] = dirf[0].x; dx[1] = dirf[1].x; dx[2] = dirf[2].x; dx[3] = dirf[3].x;
140      dy[0] = dirf[0].y; dy[1] = dirf[1].y; dy[2] = dirf[2].y; dy[3] = dirf[3].y;
141      dz[0] = dirf[0].z; dz[1] = dirf[1].z; dz[2] = dirf[2].z; dz[3] = dirf[3].z;
142      std::cerr << "Normalization not yet implemented" << std::endl;
143      abort();
144    }
145    // other components
146    ttype = _type;
147    origin = _originCell;
148    termination = NULL;
149    startObject = _startObject;
150    stopObject = _stopObject;
151    Init();
152  }
153
154  // Computes the inverted direction of the rays, used optionally by
155  // a ray traversal algorithm. This has to be reconsidered, if it
156  // is really valuable. !!!
157  void ComputeInvertedDir() const;
158
159  // Computes the sign of the rays and returns false if all the directions
160  // for all three axes are the same, but it could be different among axes,
161  // but for one axis all 4 rays must have the same direction
162  bool ComputeDirSign() const;
163 
164  // the cell in the ASDS, where ray starts from
165  void SetOrigin(const void *c) {origin = c;}
166  const void *GetOrigin() const  { return origin; }
167
168  // the cell in the ASDS, where ray finishes the walk
169  void SetTermination(const void *c) {termination = c; }
170  const void* GetTermination() const { return termination;}
171
172  // the object on which the ray starts at
173  const Intersectable* GetStartObject() const { return startObject;}
174  const Intersectable* GetStopObject() const { return stopObject;}
175
176  void SetStartObject(const Intersectable *newStartObject) {
177    startObject = newStartObject;
178  }
179  void SetStopObject(const Intersectable *newStopObject) {
180    stopObject = newStopObject;
181  }
182  int GetType() const { return ttype; }
183
184  void ApplyTransform(const Matrix4x4 &tform) {
185    for (int i = 0; i < 4; i++) {
186      Vector3 o_orig(ox[i], oy[i], oz[i]);
187      Vector3 t_orig = tform * o_orig;
188      ox[i] = t_orig.x; oy[i] = t_orig.y; oz[i] = t_orig.z;
189      Vector3 o_dir(dx[i], dy[i], dz[i]);
190      Vector3 t_dir = RotateOnly(tform, o_dir);
191      dx[i] = t_dir.x; dy[i] = t_dir.y; dz[i] = t_dir.z;
192
193      // ?? note that normalization to the unit size of the direction
194      // ?? is NOT computed -- this is what we want.
195    }
196    Precompute();
197  }
198 
199  // Reading and Setting origin of the ray and direction
200  // Ray origin
201  inline void SetLoc(int i, const Vector3 &l);
202  Vector3  GetLoc(int i) const;
203  // Direction
204  void SetDir(int i, const Vector3 &ndir);
205  Vector3  GetDir(int i) const;
206
207  // Retuns an object that is intersected by i-th ray
208  Intersectable* GetObject(int i) const;
209  void SetObject(int i, Intersectable* obj);
210  // Retuns a signed distance that is intersected by i-th ray
211  float GetT(int i) const;
212  void SetT(int i, float t);
213  // Retuns maximum signed distance that is intersected by i-th ray
214  float GetMaxT(int i) const;
215  void SetMaxT(int i, float t);
216 
217  // make such operation to slightly change the ray direction
218  // in case any component of ray direction is zero. This is
219  // carried out for all the rays in a packet
220  void CorrectZeroComponents();
221
222  // Returns the sign of the direction if this was precomputed
223  const int &GetSign(int axis) const { return sign_dir[axis];}
224
225  // Reset the result of intersection
226  void ResetObjects() {
227    obj[0] = obj[1] = obj[2] = obj[3] = 0;
228  }
229 
230private:
231  // Here it is crucial the layout of the rays in memory
232  // The layout by Jakko Biker is used
233  typedef float real;
234  union
235  {
236    struct
237    {
238      // ox[N],oy[N],oz[N] - origin of the ray N
239      union { real ox[4]; __m128 ox4; };
240      union { real oy[4]; __m128 oy4; };
241      union { real oz[4]; __m128 oz4; };
242    };
243    __m128 orig[3];
244  };
245  union
246  {
247    struct
248    {
249      // dx[N],dy[N],dz[N] - direction of the ray N
250      union { real dx[4]; __m128 dx4; };
251      union { real dy[4]; __m128 dy4; };
252      union { real dz[4]; __m128 dz4; };
253    };
254    __m128 dir[3];
255  };
256
257#define _USE_INVDIR_RP
258#ifdef _USE_INVDIR_RP
259  union
260  {
261    struct
262    {
263      // idx[N],idy[N],idz[N] - direction of the ray N
264      // inverted dir - maybe an overkill for SSE
265      // to be checked !
266      union { real idx[4]; __m128 idx4; };
267      union { real idy[4]; __m128 idy4; };
268      union { real idz[4]; __m128 idz4; };
269    };
270    __m128 idir[3];
271  };
272#endif
273
274  // The auxiliary and result values for traversal
275     
276  // Here is the result - currently computed signed distance
277  union { real t[4]; __m128 t4; };
278  // and so far minimum signed distance computed. This is required
279  // for computing ray object intersections
280  union { real tmax[4]; __m128 tmax4; };
281  // Here are the pointers to the objects that correspond to tmax[4]
282  // above. They can be different for all the rays !
283  union { Intersectable* obj[4]; __m128 obj4; };
284
285  friend class CKTBTraversal;
286 
287  // Type of the ray: primary, shadow, dummy etc., see ERayType above
288  int ttype;
289
290  // The sign of direction to be used in some algorithms. The sign
291  // has to be the same for all the rays in all the components of the
292  // direction vector !!!!
293  mutable int      sign_dir[3];
294 
295  // I should have some abstract cell data type !!! here
296  // corresponds to the spatial elementary cell
297  const void *origin;
298  const void *termination;
299 
300  // the object on which surface a ray starts from
301  const Intersectable *startObject;
302  // the object on which surface a ray ends, for computation
303  // of the visibility queries between two points
304  const Intersectable *stopObject;
305
306  /// Precompute some CRay parameters. Most of them used for ropes traversal.
307  inline void Init();
308
309  // Precompute some values that are necessary.
310  inline void Precompute();
311};
312
313// --------------------------------------------------------------------------
314//  RayPacket2x2::SetLoc()
315// --------------------------------------------------------------------------
316inline void
317RayPacket2x2::SetLoc(int i, const Vector3 &l)
318{
319  assert( (i>=0) && (i<4));
320  ox[i] = l.x;
321  oy[i] = l.y;
322  oz[i] = l.z;
323}
324
325inline void
326RayPacket2x2::SetDir(int i, const Vector3 &ndr)
327{
328  // We assume that the direction is normalized !!!
329  assert( (i>=0) && (i<4));
330  dx[i] = ndr.x;
331  dy[i] = ndr.y;
332  dz[i] = ndr.z;
333}
334
335inline Vector3
336RayPacket2x2::GetLoc(int i) const
337{
338  assert( (i>=0) && (i<4));
339  return Vector3(ox[i],oy[i],oz[i]);
340}
341
342inline Vector3
343RayPacket2x2::GetDir(int i) const
344{
345  assert( (i>=0) && (i<4));
346  return Vector3(dx[i],dy[i],dz[i]);
347}
348
349// --------------------------------------------------------------------------
350//  RayPacket2x2::Precompute()
351// --------------------------------------------------------------------------
352inline void
353RayPacket2x2::Precompute()
354{
355  // initialize inverted dir ?
356#ifdef _USE_INVDIR_RP
357  //
358  const float epsAdd = -1e-25;
359  idx[0] = 1.0f / (dx[0] + epsAdd);
360  idx[1] = 1.0f / (dx[1] + epsAdd);
361  idx[2] = 1.0f / (dx[2] + epsAdd);
362  idx[3] = 1.0f / (dx[3] + epsAdd);
363 
364  idy[0] = 1.0f / (dy[0] + epsAdd);
365  idy[1] = 1.0f / (dy[1] + epsAdd);
366  idy[2] = 1.0f / (dy[2] + epsAdd);
367  idy[3] = 1.0f / (dy[3] + epsAdd);
368 
369  idz[0] = 1.0f / (dz[0] + epsAdd);
370  idz[1] = 1.0f / (dz[1] + epsAdd);
371  idz[2] = 1.0f / (dz[2] + epsAdd);
372  idz[3] = 1.0f / (dz[3] + epsAdd);
373#endif
374}
375
376// --------------------------------------------------------------------------
377//  RayPacket2x2::Init()
378// --------------------------------------------------------------------------
379inline void
380RayPacket2x2::Init()
381{
382  // apply the standard precomputation
383  Precompute();
384}
385
386// Computes the sign of the rays and returns false if all the directions
387// for all three axes are the same, but it could be different among axes,
388// but for one axis all 4 rays must have the same direction
389inline bool
390RayPacket2x2::ComputeDirSign() const
391{
392  // Set the sign of the direction 1 when negative
393  sign_dir[0] = (dx[0] < 0.0f);
394  sign_dir[1] = (dy[0] < 0.0f);
395  sign_dir[2] = (dz[0] < 0.0f);
396  for (int i = 1; i < 4; i++) {
397    if (sign_dir[0] != (dx[i] < 0.0f))
398      return true; // different direction in x
399    if (sign_dir[1] != (dy[i] < 0.0f))
400      return true; // different direction in y
401    if (sign_dir[2] != (dz[i] < 0.0f))
402      return true; // different direction in z
403  }// for
404
405  // Returns false if all 4 rays have the consistent direction
406  return false;
407}
408
409inline Intersectable*
410RayPacket2x2::GetObject(int i) const
411{
412  assert( (i>=0) && (i<4));
413  return obj[i];
414}
415
416inline void
417RayPacket2x2::SetObject(int i, Intersectable* object)
418{
419  assert( (i>=0) && (i<4));
420  obj[i] = object;
421}
422
423inline float
424RayPacket2x2::GetT(int i) const
425{
426  assert( (i>=0) && (i<4));
427  return t[i];
428}
429
430inline void
431RayPacket2x2::SetT(int i, float tnew)
432{
433  assert( (i>=0) && (i<4));
434  t[i] = tnew;
435}
436
437inline float
438RayPacket2x2::GetMaxT(int i) const
439{
440  assert( (i>=0) && (i<4));
441  return tmax[i];
442}
443
444inline void
445RayPacket2x2::SetMaxT(int i, float tmaxnew)
446{
447  assert( (i>=0) && (i<4));
448  tmax[i] = tmaxnew;
449}
450#endif // __SSE__
451
452#endif // _USE_HAVRAN_SSE
453 
454} // namespace
455
456#endif // __RAYPACK_H__
457
Note: See TracBrowser for help on using the repository browser.