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

// ============================================================================
// $Id: raypack.h $
//
// raypack.h
//        CRayPacket class - core of the ray-packet traversal routines
//
// Class: CRayPacket2x2
//
// REPLACEMENT_STRING
//
// Initial coding by Vlastimil Havran, 2006. The data design is in fact
// Jakko Biker layout as proposed in the article on Intel Web Site in year 2005
// http://www.intel.com/cd/ids/developer/asmo-na/eng/245711.htm?page=1

#ifndef __RAYPACK_H__
#define __RAYPACK_H__

#include <cassert>

namespace GtpVisibilityPreprocessor {

#include "Vector3.h"
#include "Matrix4x4.h"
#include "ktbconf.h"
  
#ifdef _USE_HAVRAN_SSE
  
#ifdef __SSE__

// System headers for SSE
#ifdef __INTEL_COMPILER
#include <xmmintrin.h>
#else
// We assume GNU GCC compiler 3.4 or higher
#include <xmmintrin.h>
#endif


// forward declarations
#define SSE_INTRINSIC
#ifdef SSE_INTRINSIC
#define ALIGN16             __declspec(align(16)) 

#ifdef _MSC_VER
#define NEW_ALIGN16(type,n)  ((type*)_aligned_malloc((n)*sizeof(type),16))
#define FREE_ALIGN16(array)  if(array){_aligned_free(array);(array)=0;}
#else
#define NEW_ALIGN16(type,n)  ((type*)malloc((n)*sizeof(type)))
#define FREE_ALIGN16(array)  if (array) { free(array);(array)=0;}
#endif // _MSC_VC

#define PAD_FOUR(h)           mulFour(h)
union extract_m128
{
  __m128 m;
  float f[4];
};
#else
#define NEW_ALIGN16(type,n)  ((type*)malloc((n)*sizeof(type)))
#define ALIGN16              
#define FREE_ALIGN16(array)  if(array){free(array);(array)=0;}
#define PAD_FOUR(h)           (h)
#endif
  
// -------------------------------------------------------------------
// RayPacket2x2 class.  A set of 4 ray is defined by a location and a
// direction. The direction is always normalized (length == 1) during
// use.
// -------------------------------------------------------------------
class RayPacket2x2
{
public:
  enum {
    // The number of rays in one packet
    PACKSIZE = 4
  };
  
  // constructor
  RayPacket2x2(
    // origin and the direction of rays
    const Vector3 orf[],
    const Vector3 dirf[],
    // The same type for all four rays
    const int _type,
    // All the rays has to start at the same cell
    const void *_originCell = NULL,
    // All the rays has to start at the same object
    const Intersectable *_startObject = NULL,
    // and also for shadow rays finish at the same object
    const Intersectable *_stopObject = NULL
  ) {
    // location
    ox[0] = orf[0].x; ox[1] = orf[1].x; ox[2] = orf[2].x; ox[3] = orf[3].x;
    oy[0] = orf[0].y; oy[1] = orf[1].y; oy[2] = orf[2].y; oy[3] = orf[3].y;
    oz[0] = orf[0].z; oz[1] = orf[1].z; oz[2] = orf[2].z; oz[3] = orf[3].z;
    // direction, we assume to be normalized
    dx[0] = dirf[0].x; dx[1] = dirf[1].x; dx[2] = dirf[2].x; dx[3] = dirf[3].x;
    dy[0] = dirf[0].y; dy[1] = dirf[1].y; dy[2] = dirf[2].y; dy[3] = dirf[3].y;
    dz[0] = dirf[0].z; dz[1] = dirf[1].z; dz[2] = dirf[2].z; dz[3] = dirf[3].z;
    // Other components
    ttype = _type;
    origin = _originCell;
    termination = NULL;
    startObject = _startObject;
    stopObject = _stopObject;
    Init();
  }
  // dummy constructor
  RayPacket2x2() {}

  // Inititalize the ray again when already constructed
  void Init(
    // origin and the direction of rays
    const Vector3 orf[],
    const Vector3 dirf[],
    // The same type for all four rays
    const int _type,
    // All the rays has to start at the same cell
    const void *_originCell = NULL,
    // All the rays has to start at the same object
    const Intersectable *_startObject = NULL,
    // and also for shadow rays finish at the same object
    const Intersectable *_stopObject = NULL,
    // if the direction of rays is normalized or not
    bool dirNormalized = false)
  {
    // location
    ox[0] = orf[0].x; ox[1] = orf[1].x; ox[2] = orf[2].x; ox[3] = orf[3].x;
    oy[0] = orf[0].y; oy[1] = orf[1].y; oy[2] = orf[2].y; oy[3] = orf[3].y;
    oz[0] = orf[0].z; oz[1] = orf[1].z; oz[2] = orf[2].z; oz[3] = orf[3].z;
    // direction
    if (dirNormalized) {
      // already normalized
      dx[0] = dirf[0].x; dx[1] = dirf[1].x; dx[2] = dirf[2].x; dx[3] = dirf[3].x;
      dy[0] = dirf[0].y; dy[1] = dirf[1].y; dy[2] = dirf[2].y; dy[3] = dirf[3].y;
      dz[0] = dirf[0].z; dz[1] = dirf[1].z; dz[2] = dirf[2].z; dz[3] = dirf[3].z;
    }
    else {
      dx[0] = dirf[0].x; dx[1] = dirf[1].x; dx[2] = dirf[2].x; dx[3] = dirf[3].x;
      dy[0] = dirf[0].y; dy[1] = dirf[1].y; dy[2] = dirf[2].y; dy[3] = dirf[3].y;
      dz[0] = dirf[0].z; dz[1] = dirf[1].z; dz[2] = dirf[2].z; dz[3] = dirf[3].z;
      std::cerr << "Normalization not yet implemented" << std::endl;
      abort();
    }
    // other components
    ttype = _type;
    origin = _originCell;
    termination = NULL;
    startObject = _startObject;
    stopObject = _stopObject;
    Init();
  }

  // Computes the inverted direction of the rays, used optionally by
  // a ray traversal algorithm. This has to be reconsidered, if it
  // is really valuable. !!!
  void ComputeInvertedDir() const;

  // Computes the sign of the rays and returns false if all the directions
  // for all three axes are the same, but it could be different among axes,
  // but for one axis all 4 rays must have the same direction
  bool ComputeDirSign() const;
  
  // the cell in the ASDS, where ray starts from
  void SetOrigin(const void *c) {origin = c;}
  const void *GetOrigin() const  { return origin; }

  // the cell in the ASDS, where ray finishes the walk
  void SetTermination(const void *c) {termination = c; }
  const void* GetTermination() const { return termination;}

  // the object on which the ray starts at
  const Intersectable* GetStartObject() const { return startObject;}
  const Intersectable* GetStopObject() const { return stopObject;}

  void SetStartObject(const Intersectable *newStartObject) {
    startObject = newStartObject;
  }
  void SetStopObject(const Intersectable *newStopObject) {
    stopObject = newStopObject;
  }
  int GetType() const { return ttype; }

  void ApplyTransform(const Matrix4x4 &tform) {
    for (int i = 0; i < 4; i++) {
      Vector3 o_orig(ox[i], oy[i], oz[i]);
      Vector3 t_orig = tform * o_orig;
      ox[i] = t_orig.x; oy[i] = t_orig.y; oz[i] = t_orig.z;
      Vector3 o_dir(dx[i], dy[i], dz[i]);
      Vector3 t_dir = RotateOnly(tform, o_dir);
      dx[i] = t_dir.x; dy[i] = t_dir.y; dz[i] = t_dir.z;

      // ?? note that normalization to the unit size of the direction
      // ?? is NOT computed -- this is what we want.
    }
    Precompute();
  }
  
  // Reading and Setting origin of the ray and direction
  // Ray origin
  inline void SetLoc(int i, const Vector3 &l);
  Vector3  GetLoc(int i) const;
  // Direction
  void SetDir(int i, const Vector3 &ndir);
  Vector3  GetDir(int i) const;

  // Retuns an object that is intersected by i-th ray
  Intersectable* GetObject(int i) const;
  void SetObject(int i, Intersectable* obj);
  // Retuns a signed distance that is intersected by i-th ray
  float GetT(int i) const;
  void SetT(int i, float t);
  // Retuns maximum signed distance that is intersected by i-th ray
  float GetMaxT(int i) const;
  void SetMaxT(int i, float t);
  
  // make such operation to slightly change the ray direction
  // in case any component of ray direction is zero. This is
  // carried out for all the rays in a packet
  void CorrectZeroComponents();

  // Returns the sign of the direction if this was precomputed
  const int &GetSign(int axis) const { return sign_dir[axis];}

  // Reset the result of intersection
  void ResetObjects() {
    obj[0] = obj[1] = obj[2] = obj[3] = 0;
  }
  
private:
  // Here it is crucial the layout of the rays in memory
  // The layout by Jakko Biker is used
  typedef float real;
  union
  {
    struct
    {
      // ox[N],oy[N],oz[N] - origin of the ray N
      union { real ox[4]; __m128 ox4; };
      union { real oy[4]; __m128 oy4; };
      union { real oz[4]; __m128 oz4; };
    };
    __m128 orig[3];
  };
  union
  {
    struct
    {
      // dx[N],dy[N],dz[N] - direction of the ray N
      union { real dx[4]; __m128 dx4; };
      union { real dy[4]; __m128 dy4; };
      union { real dz[4]; __m128 dz4; };
    };
    __m128 dir[3];
  };

#define _USE_INVDIR_RP
#ifdef _USE_INVDIR_RP
  union
  {
    struct
    {
      // idx[N],idy[N],idz[N] - direction of the ray N
      // inverted dir - maybe an overkill for SSE
      // to be checked !
      union { real idx[4]; __m128 idx4; };
      union { real idy[4]; __m128 idy4; };
      union { real idz[4]; __m128 idz4; };
    };
    __m128 idir[3];
  };
#endif

  // The auxiliary and result values for traversal
      
  // Here is the result - currently computed signed distance
  union { real t[4]; __m128 t4; };
  // and so far minimum signed distance computed. This is required
  // for computing ray object intersections
  union { real tmax[4]; __m128 tmax4; };
  // Here are the pointers to the objects that correspond to tmax[4]
  // above. They can be different for all the rays !
  union { Intersectable* obj[4]; __m128 obj4; };

  friend class CKTBTraversal;
  
  // Type of the ray: primary, shadow, dummy etc., see ERayType above
  int ttype;

  // The sign of direction to be used in some algorithms. The sign
  // has to be the same for all the rays in all the components of the
  // direction vector !!!!
  mutable int      sign_dir[3];
 
  // I should have some abstract cell data type !!! here
  // corresponds to the spatial elementary cell
  const void *origin;
  const void *termination;
  
  // the object on which surface a ray starts from
  const Intersectable *startObject;
  // the object on which surface a ray ends, for computation
  // of the visibility queries between two points
  const Intersectable *stopObject;

  /// Precompute some CRay parameters. Most of them used for ropes traversal.
  inline void Init();

  // Precompute some values that are necessary.
  inline void Precompute();
};

// --------------------------------------------------------------------------
//  RayPacket2x2::SetLoc()
// --------------------------------------------------------------------------
inline void
RayPacket2x2::SetLoc(int i, const Vector3 &l)
{
  assert( (i>=0) && (i<4));
  ox[i] = l.x;
  oy[i] = l.y;
  oz[i] = l.z;
}

inline void
RayPacket2x2::SetDir(int i, const Vector3 &ndr)
{
  // We assume that the direction is normalized !!!
  assert( (i>=0) && (i<4));
  dx[i] = ndr.x;
  dy[i] = ndr.y;
  dz[i] = ndr.z;
}

inline Vector3
RayPacket2x2::GetLoc(int i) const
{
  assert( (i>=0) && (i<4));
  return Vector3(ox[i],oy[i],oz[i]);
}

inline Vector3
RayPacket2x2::GetDir(int i) const
{
  assert( (i>=0) && (i<4));
  return Vector3(dx[i],dy[i],dz[i]);
}

// --------------------------------------------------------------------------
//  RayPacket2x2::Precompute()
// --------------------------------------------------------------------------
inline void
RayPacket2x2::Precompute()
{
  // initialize inverted dir ?
#ifdef _USE_INVDIR_RP
  //
  const float epsAdd = -1e-25;
  idx[0] = 1.0f / (dx[0] + epsAdd);
  idx[1] = 1.0f / (dx[1] + epsAdd);
  idx[2] = 1.0f / (dx[2] + epsAdd);
  idx[3] = 1.0f / (dx[3] + epsAdd);
  
  idy[0] = 1.0f / (dy[0] + epsAdd);
  idy[1] = 1.0f / (dy[1] + epsAdd);
  idy[2] = 1.0f / (dy[2] + epsAdd);
  idy[3] = 1.0f / (dy[3] + epsAdd);
  
  idz[0] = 1.0f / (dz[0] + epsAdd);
  idz[1] = 1.0f / (dz[1] + epsAdd);
  idz[2] = 1.0f / (dz[2] + epsAdd);
  idz[3] = 1.0f / (dz[3] + epsAdd);
#endif
}

// --------------------------------------------------------------------------
//  RayPacket2x2::Init()
// --------------------------------------------------------------------------
inline void
RayPacket2x2::Init()
{
  // apply the standard precomputation
  Precompute();
}

// Computes the sign of the rays and returns false if all the directions
// for all three axes are the same, but it could be different among axes,
// but for one axis all 4 rays must have the same direction
inline bool
RayPacket2x2::ComputeDirSign() const
{
  // Set the sign of the direction 1 when negative
  sign_dir[0] = (dx[0] < 0.0f);
  sign_dir[1] = (dy[0] < 0.0f);
  sign_dir[2] = (dz[0] < 0.0f);
  for (int i = 1; i < 4; i++) {
    if (sign_dir[0] != (dx[i] < 0.0f))
      return true; // different direction in x 
    if (sign_dir[1] != (dy[i] < 0.0f))
      return true; // different direction in y
    if (sign_dir[2] != (dz[i] < 0.0f))
      return true; // different direction in z
  }// for

  // Returns false if all 4 rays have the consistent direction
  return false;
}

inline Intersectable*
RayPacket2x2::GetObject(int i) const
{
  assert( (i>=0) && (i<4));
  return obj[i];
}

inline void
RayPacket2x2::SetObject(int i, Intersectable* object)
{
  assert( (i>=0) && (i<4));
  obj[i] = object;
}

inline float
RayPacket2x2::GetT(int i) const
{
  assert( (i>=0) && (i<4));
  return t[i];
}

inline void
RayPacket2x2::SetT(int i, float tnew)
{
  assert( (i>=0) && (i<4));
  t[i] = tnew;
}

inline float
RayPacket2x2::GetMaxT(int i) const
{
  assert( (i>=0) && (i<4));
  return tmax[i];
}

inline void
RayPacket2x2::SetMaxT(int i, float tmaxnew)
{
  assert( (i>=0) && (i<4));
  tmax[i] = tmaxnew;
}
#endif // __SSE__

#endif // _USE_HAVRAN_SSE
  
} // namespace

#endif // __RAYPACK_H__