DGtal  0.9.2
RayIntersectionPredicates.h
1 
17 #pragma once
18 
31 #if defined(RayIntersectionPredicates_RECURSES)
32 #error Recursive header files inclusion detected in RayIntersectionPredicates.h
33 #else // defined(RayIntersectionPredicates_RECURSES)
34 
35 #define RayIntersectionPredicates_RECURSES
36 
37 #if !defined RayIntersectionPredicates_h
38 
39 #define RayIntersectionPredicates_h
40 
42 // Inclusions
43 #include <iostream>
44 #include "DGtal/base/Common.h"
46 
47 namespace DGtal
48 {
49 
80  template <typename TPoint>
82  {
83 
85  BOOST_STATIC_ASSERT( TPoint::dimension == 3);
86 
88  typedef TPoint Point;
89 
91  typedef TPoint Vector;
92 
94  typedef typename TPoint::Component Component;
95 
105  RayIntersectionPredicate( const Point &origin,
106  const Vector &dest)
107  : myOrigin(origin), myDest(dest)
108  {
109  ASSERT_MSG( dest.norm1() != NumberTraits<typename Point::UnsignedComponent>::ZERO,
110  "Direction must be non-null vector");
111  }
112 
125  bool operator()(const Point &v1,
126  const Point &v2,
127  const Point &v3) const
128  {
129 
130  ASSERT((v1 != v2 ) && (v1 != v3) && (v2 != v3));
131 
132  Point e1, e2; //Edge1, Edge2
133  Point P, Q, T;
134  Component det, u, v;
135 
136  //Find vectors for two edges sharing V1
137  e1 = v2 - v1;
138  e2 = v3 - v1;
139 
140  //Begin calculating determinant - also used to calculate u parameter
141  P = myDest.crossProduct( e2 );
142 
143  //if determinant is near zero, ray lies in plane of triangle
144  det = e1.dot( P );
146  {
147  return false;
148  }
149 
150  //calculate distance from V1 to ray origin
151  T = myOrigin - v1;
152 
153  //Calculate u parameter and test bound
154  u = T.dot( P ); //* inv_det;
155 
157  {
158  if ((u < NumberTraits<Component>::ZERO) ||
159  ( u > det))
160  {
161  return false;
162  }
163  }
164  else
165  {
166  if ((u > NumberTraits<Component>::ZERO) ||
167  ( u < det))
168  {
169  return false;
170  }
171  }
172 
173  //Prepare to test v parameter
174  Q = T.crossProduct( e1 );
175 
176  //Calculate V parameter and test bound
177  v = myDest.dot( Q );
178 
179  //The intersection lies outside of the triangle
181  {
182  if ((v < NumberTraits<Component>::ZERO) ||
183  ((u+v) > det))
184  {
185  return false;
186  }
187  }
188  else
189  {
190  if ((v > NumberTraits<Component>::ZERO) ||
191  ((u+v) < det))
192  {
193  return false;
194  }
195  }
196 
197  //distance to triangle must be positive
198  Component t = e2.dot( Q ) ;
199  if (t*det < NumberTraits<Component>::ZERO)
200  return false;
201 
202  return true;
203  }
204 
217  bool operator()(const Point &v1,
218  const Point &v2,
219  const Point &v3,
220  const Point &v4) const
221  {
222  return (this->operator()(v1,v2,v3) ||
223  this->operator()(v1,v4,v3) );
224  }
225 
239  template < typename Surfel >
240  bool operator()(const Surfel &aSurfel) const
241  {
242  auto const & aPreSurfel = aSurfel.preCell();
243 
244  Component x1,x2,x3,x4;
245  Component y1,y2,y3,y4;
246  Component z1,z2,z3,z4;
247  Component ONE = NumberTraits<Component>::ONE;
248 
249  Point baseQuadCenter = aPreSurfel.coordinates;
250 
251  bool yodd = ( NumberTraits<Component>::castToInt64_t(aPreSurfel.coordinates[ 1 ]) & 1 );
252  bool zodd = ( NumberTraits<Component>::castToInt64_t(aPreSurfel.coordinates[ 2 ]) & 1 );
253 
254  if(!zodd)
255  {
256  //zsurfel
257  x1= baseQuadCenter[0]-ONE; y1= baseQuadCenter[1]-ONE; z1= baseQuadCenter[2];
258  x2= baseQuadCenter[0]+ONE; y2= baseQuadCenter[1]-ONE; z2= baseQuadCenter[2];
259  x3= baseQuadCenter[0]+ONE; y3= baseQuadCenter[1]+ONE; z3= baseQuadCenter[2];
260  x4= baseQuadCenter[0]-ONE; y4= baseQuadCenter[1]+ONE; z4= baseQuadCenter[2];
261  }
262  else if(!yodd)
263  {
264  //ysurfel
265  x1= baseQuadCenter[0]-ONE; y1= baseQuadCenter[1]; z1= baseQuadCenter[2]-ONE;
266  x2= baseQuadCenter[0]-ONE; y2= baseQuadCenter[1]; z2= baseQuadCenter[2]+ONE;
267  x3= baseQuadCenter[0]+ONE; y3= baseQuadCenter[1]; z3= baseQuadCenter[2]+ONE;
268  x4= baseQuadCenter[0]+ONE; y4= baseQuadCenter[1]; z4= baseQuadCenter[2]-ONE;
269  }
270  else
271  {
272  //xsurfel
273  x1= baseQuadCenter[0]; y1= baseQuadCenter[1]-ONE; z1= baseQuadCenter[2]-ONE;
274  x2= baseQuadCenter[0]; y2= baseQuadCenter[1]+ONE; z2= baseQuadCenter[2]-ONE;
275  x3= baseQuadCenter[0]; y3= baseQuadCenter[1]+ONE; z3= baseQuadCenter[2]+ONE;
276  x4= baseQuadCenter[0]; y4= baseQuadCenter[1]-ONE; z4= baseQuadCenter[2]+ONE;
277  }
278  return this->operator()(Point(x1, y1, z1), Point(x2 ,y2, z2),
279  Point(x3, y3, z3), Point(x4, y4, z4));
280  }
281 
282  Point myOrigin;
283  Point myDest;
284  };
285 
286 
287 
288 } // namespace DGtal
289 
290 // //
292 
293 #endif // !defined RayIntersectionPredicates_h
294 
295 #undef RayIntersectionPredicates_RECURSES
296 #endif // else defined(RayIntersectionPredicates_RECURSES)
TPoint::Component Component
Type of point coordinates.
bool operator()(const Surfel &aSurfel) const
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:69
static DGtal::int64_t castToInt64_t(const T &aT)
Definition: NumberTraits.h:146
DGtal is the top-level namespace which contains all DGtal functions and types.
bool operator()(const Point &v1, const Point &v2, const Point &v3) const
RayIntersectionPredicate(const Point &origin, const Vector &dest)
This class implements various intersection predicates between a ray and a triangle, a quad or a surfel in dimension 3.
bool operator()(const Point &v1, const Point &v2, const Point &v3, const Point &v4) const
BOOST_STATIC_ASSERT(TPoint::dimension==3)
Only in dimension 3.