RayIntersectionPredicates.h

 
#pragma once
 
#if defined(RayIntersectionPredicates_RECURSES)
#error Recursive header files inclusion detected in RayIntersectionPredicates.h
#else // defined(RayIntersectionPredicates_RECURSES)
#define RayIntersectionPredicates_RECURSES
 
#if !defined RayIntersectionPredicates_h
#define RayIntersectionPredicates_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
 
namespace DGtal
{
 
  template <typename TPoint>
  struct RayIntersectionPredicate
  {
 
    BOOST_STATIC_ASSERT( TPoint::dimension == 3);
 
    typedef TPoint Point;
    
    typedef TPoint Vector;
    
    typedef typename TPoint::Component Component;
 
    RayIntersectionPredicate( const Point &origin, 
                              const Vector &dest)
      : myOrigin(origin), myDest(dest) 
    {
      ASSERT_MSG( dest.norm1() != NumberTraits<typename Point::UnsignedComponent>::ZERO, 
                  "Direction must be non-null vector"); 
    }
 
    bool operator()(const Point &v1, 
                    const Point &v2,
                    const Point &v3) const
    {
      
      ASSERT((v1 != v2 ) && (v1 != v3) && (v2 != v3));
 
      Point e1, e2;  //Edge1, Edge2
      Point P, Q, T;
      Component det, u, v;
 
      //Find vectors for two edges sharing V1
      e1 = v2 - v1;
      e2 = v3 - v1;
 
      //Begin calculating determinant - also used to calculate u parameter
      P = myDest.crossProduct( e2 );
      
      //if determinant is near zero, ray lies in plane of triangle
      det = e1.dot( P );
      if(det == NumberTraits<Component>::ZERO) 
        {
          return false;
        }
      
      //calculate distance from V1 to ray origin
      T =  myOrigin -  v1;
      
      //Calculate u parameter and test bound
      u = T.dot( P );  //* inv_det;
      
      if (det >  NumberTraits<Component>::ZERO)
        {
          if ((u < NumberTraits<Component>::ZERO) ||
              ( u > det))
            {
              return false;
            }
        }
      else
        {
          if ((u > NumberTraits<Component>::ZERO) ||
              ( u < det))
            {
              return false;
            }
        }
 
      //Prepare to test v parameter
      Q = T.crossProduct( e1 );
      
      //Calculate V parameter and test bound
      v = myDest.dot( Q ); 
      
      //The intersection lies outside of the triangle
      if (det >  NumberTraits<Component>::ZERO)
        {
          if ((v < NumberTraits<Component>::ZERO)  ||
              ((u+v) > det)) 
            {
              return false;
            }       
        }
      else
        {
          if ((v > NumberTraits<Component>::ZERO)  ||
              ((u+v) < det)) 
            {
              return false;
            }       
        }
 
      //distance to triangle must be positive
      Component t = e2.dot( Q ) ;
      if (t*det < NumberTraits<Component>::ZERO)
        return false;
       
      return true;
    }
 
    bool operator()(const Point &v1, 
                    const Point &v2,
                    const Point &v3,
                    const Point &v4) const
    {
      return (this->operator()(v1,v2,v3) ||
              this->operator()(v1,v4,v3) );
    }
 
    template < typename Surfel >
    bool operator()(const Surfel &aSurfel) const
    {
      auto const & aPreSurfel = aSurfel.preCell();
 
      Component x1,x2,x3,x4;
      Component y1,y2,y3,y4;
      Component z1,z2,z3,z4;
      Component ONE = NumberTraits<Component>::ONE;
      
      Point baseQuadCenter =  aPreSurfel.coordinates;
      
      bool yodd = ( NumberTraits<Component>::castToInt64_t(aPreSurfel.coordinates[ 1 ]) & 1 );
      bool zodd = ( NumberTraits<Component>::castToInt64_t(aPreSurfel.coordinates[ 2 ]) & 1 ); 
      
      if(!zodd)
        {
          //zsurfel
          x1= baseQuadCenter[0]-ONE; y1= baseQuadCenter[1]-ONE; z1= baseQuadCenter[2];
          x2= baseQuadCenter[0]+ONE; y2= baseQuadCenter[1]-ONE; z2= baseQuadCenter[2];
          x3= baseQuadCenter[0]+ONE; y3= baseQuadCenter[1]+ONE; z3= baseQuadCenter[2];
          x4= baseQuadCenter[0]-ONE; y4= baseQuadCenter[1]+ONE; z4= baseQuadCenter[2];
        }
      else if(!yodd)
        {
          //ysurfel
          x1= baseQuadCenter[0]-ONE; y1= baseQuadCenter[1]; z1= baseQuadCenter[2]-ONE;
          x2= baseQuadCenter[0]-ONE; y2= baseQuadCenter[1]; z2= baseQuadCenter[2]+ONE;
          x3= baseQuadCenter[0]+ONE; y3= baseQuadCenter[1]; z3= baseQuadCenter[2]+ONE;
          x4= baseQuadCenter[0]+ONE; y4= baseQuadCenter[1]; z4= baseQuadCenter[2]-ONE;
        }
      else
        {
          //xsurfel
          x1= baseQuadCenter[0]; y1= baseQuadCenter[1]-ONE; z1= baseQuadCenter[2]-ONE;
          x2= baseQuadCenter[0]; y2= baseQuadCenter[1]+ONE; z2= baseQuadCenter[2]-ONE;
          x3= baseQuadCenter[0]; y3= baseQuadCenter[1]+ONE; z3= baseQuadCenter[2]+ONE;
          x4= baseQuadCenter[0]; y4= baseQuadCenter[1]-ONE; z4= baseQuadCenter[2]+ONE;
        }
      return this->operator()(Point(x1, y1, z1), Point(x2 ,y2, z2),
                              Point(x3, y3, z3), Point(x4, y4, z4)); 
    }
    
    Point myOrigin;
    Point myDest;    
  }; 
 
 
 
} // namespace DGtal
 
//                                                                           //
 
#endif // !defined RayIntersectionPredicates_h
 
#undef RayIntersectionPredicates_RECURSES
#endif // else defined(RayIntersectionPredicates_RECURSES)