RegionOfInterest.cpp
Go to the documentation of this file.
1 // (C) Copyright Renaud Detry 2007-2015.
2 // Distributed under the GNU General Public License and under the
3 // BSD 3-Clause License (See accompanying file LICENSE.txt).
4 
5 /** @file */
6 
7 /**
8  * @file
9  */
10 
11 #include <nuklei/RegionOfInterest.h>
12 #include <nuklei/KernelCollection.h>
13 #include <sstream>
14 #include <boost/assign/std/vector.hpp>
15 
16 #ifdef NUKLEI_USE_CGAL
17 
18 #ifndef CGAL_EIGEN3_ENABLED
19 # define CGAL_EIGEN3_ENABLED 1
20 #endif
21 
22 #include <CGAL/Simple_cartesian.h>
23 #include <CGAL/AABB_tree.h>
24 #include <CGAL/AABB_traits.h>
25 #include <CGAL/AABB_triangle_primitive.h>
26 #include <CGAL/property_map.h>
27 #include <CGAL/IO/read_off_points.h>
28 #include <CGAL/IO/read_xyz_points.h>
29 #include <CGAL/IO/write_xyz_points.h>
30 #include <CGAL/IO/Polyhedron_iostream.h>
31 
32 #include <CGAL/squared_distance_3.h>
33 #include <CGAL/Point_with_normal_3.h>
34 
35 #endif
36 
37 namespace nuklei
38 {
39 #ifdef NUKLEI_USE_CGAL
40  typedef CGAL::Simple_cartesian<double> K;
41  typedef K::FT FT;
42  typedef K::Ray_3 Ray;
43  typedef K::Line_3 Line;
44  typedef K::Point_3 Point;
45  typedef K::Plane_3 Plane;
46  typedef K::Vector_3 Vector;
47  typedef K::Segment_3 Segment;
48  typedef K::Triangle_3 Triangle;
49  typedef std::list<Triangle>::iterator Iterator;
50  typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
51  typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
52  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
53 #if CGAL_VERSION_NR >= 1040300000
54  typedef boost::optional< Tree::Intersection_and_primitive_id<Plane>::Type > Plane_intersection;
55 #else
56  typedef boost::optional< Tree::Object_and_primitive_id > Plane_intersection;
57 #endif
58 #endif
59 
60  SphereROI::SphereROI(const std::string &centerAndRadius)
61  {
62  std::istringstream is(centerAndRadius);
63  for (unsigned int i = 0; i < 3; ++i)
64  NUKLEI_ASSERT(is >> center_[i]);
65  NUKLEI_ASSERT(is >> radius_);
66  NUKLEI_ASSERT(radius_ >= 0);
67  }
68 
69  BoxROI::BoxROI(const std::string &s)
70  {
71  setCenterOriSize(s);
72  }
73 
74  void BoxROI::setCenterAxesSize(const std::string &centerSize)
75  {
76  std::istringstream is(centerSize);
77  Vector3 center, axis1, axis2;
78  for (unsigned int i = 0; i < 3; ++i)
79  NUKLEI_ASSERT(is >> center[i]);
80  for (unsigned int i = 0; i < 3; ++i)
81  NUKLEI_ASSERT(is >> axis1[i]);
82  for (unsigned int i = 0; i < 3; ++i)
83  NUKLEI_ASSERT(is >> axis2[i]);
84  NUKLEI_ASSERT(is >> edgeLengths_[2]);
85 
86  axis1 -= center; axis2 -= center;
87  edgeLengths_[0] = axis1.Length()*2*1.2;
88  edgeLengths_[1] = axis2.Length()*2*1.2;
89  edgeLengths_[2] *= 2*1.0;
90  axis1.Normalize(); axis2.Normalize();
91 
92  Matrix3 ori;
93  ori.SetColumn(0, axis1);
94  ori.SetColumn(1, axis2);
95  ori.SetColumn(2, la::normalized(axis1.Cross(axis2)));
96  ori.Orthonormalize();
97 
98  centerLoc_ = center;
99  la::copyRotation(centerOri_, ori);
100  }
101 
102  void BoxROI::setCenterOriSize(const std::string &centerSize)
103  {
104  std::istringstream is(centerSize);
105  Quaternion ori;
106  for (unsigned int i = 0; i < 3; ++i)
107  NUKLEI_ASSERT(is >> centerLoc_[i]);
108  for (unsigned int i = 0; i < 4; ++i)
109  NUKLEI_ASSERT(is >> ori[i]);
110  for (unsigned int i = 0; i < 3; ++i)
111  NUKLEI_ASSERT(is >> edgeLengths_[i]);
112  la::copyRotation(centerOri_, ori);
113  }
114 
115  bool SphereROI::contains_(const Vector3 &v) const
116  {
117  return dist<groupS::r3>::d(v, center_) <= radius_;
118  }
119 
120  bool BoxROI::contains_(const Vector3 &v) const
121  {
122  Vector3 projection = la::project(centerLoc_, centerOri_, v);
123 
124  return std::fabs(projection[0]) <= edgeLengths_[0]/2. &&
125  std::fabs(projection[1]) <= edgeLengths_[1]/2. &&
126  std::fabs(projection[2]) <= edgeLengths_[2]/2.;
127  }
128 
129  struct BoxFace
130  {
131  public:
132  BoxFace(const int a,
133  const int b,
134  const int c,
135  const int d) :
136  a(a), b(b), c(c), d(d)
137  {
138 
139  }
140 
141  int a, b, c, d;
142  };
143 
144  BoxFace make_boxface(const int a,
145  const int b,
146  const int c,
147  const int d)
148  {
149  return BoxFace(a, b, c, d);
150  }
151 
152  void BoxROI::pushTriangles_(std::list<boost::any>& triangles) const
153  {
154 #ifdef NUKLEI_USE_CGAL
155  KernelCollection corners;
156  {
157  kernel::r3 k;
158  k.loc_ = Vector3(edgeLengths_.X()/2 , + edgeLengths_.Y()/2 , + edgeLengths_.Z()/2);
159  corners.add(k);
160  }
161  {
162  kernel::r3 k;
163  k.loc_ = Vector3(edgeLengths_.X()/2 , + edgeLengths_.Y()/2 , - edgeLengths_.Z()/2);
164  corners.add(k);
165  }
166  {
167  kernel::r3 k;
168  k.loc_ = Vector3(edgeLengths_.X()/2 , - edgeLengths_.Y()/2 , + edgeLengths_.Z()/2);
169  corners.add(k);
170  }
171  {
172  kernel::r3 k;
173  k.loc_ = Vector3(edgeLengths_.X()/2 , - edgeLengths_.Y()/2 , - edgeLengths_.Z()/2);
174  corners.add(k);
175  }
176  {
177  kernel::r3 k;
178  k.loc_ = Vector3(- edgeLengths_.X()/2 , + edgeLengths_.Y()/2 , + edgeLengths_.Z()/2);
179  corners.add(k);
180  }
181  {
182  kernel::r3 k;
183  k.loc_ = Vector3(- edgeLengths_.X()/2 , + edgeLengths_.Y()/2 , - edgeLengths_.Z()/2);
184  corners.add(k);
185  }
186  {
187  kernel::r3 k;
188  k.loc_ = Vector3(- edgeLengths_.X()/2 , - edgeLengths_.Y()/2 , + edgeLengths_.Z()/2);
189  corners.add(k);
190  }
191  {
192  kernel::r3 k;
193  k.loc_ = Vector3(- edgeLengths_.X()/2 , - edgeLengths_.Y()/2 , - edgeLengths_.Z()/2);
194  corners.add(k);
195  }
196 
197  corners.transformWith(centerLoc_, centerOri_);
198 
199  std::vector<BoxFace> faces;
200  {
201  using namespace boost::assign;
202  faces += make_boxface(0, 1, 3, 2),
203  make_boxface(2, 3, 7, 6),
204  make_boxface(6, 7, 5, 4),
205  make_boxface(4, 5, 1, 0),
206  make_boxface(0, 1, 3, 2),
207  make_boxface(0, 2, 6, 4),
208  make_boxface(1, 3, 7, 5);
209  }
210 
211  for (std::vector<BoxFace>::const_iterator i = faces.begin();
212  i != faces.end(); ++i)
213  {
214  const Vector3& aa = corners.at(i->a).getLoc();
215  const Vector3& bb = corners.at(i->b).getLoc();
216  const Vector3& cc = corners.at(i->c).getLoc();
217  const Vector3& dd = corners.at(i->d).getLoc();
218 
219  Point a(aa.X(), aa.Y(), aa.Z());
220  Point b(bb.X(), bb.Y(), bb.Z());
221  Point c(cc.X(), cc.Y(), cc.Z());
222  Point d(dd.X(), dd.Y(), dd.Z());
223 
224  triangles.push_back(Triangle(a,b,c));
225  triangles.push_back(Triangle(a,d,c));
226  }
227 #else
228  NUKLEI_THROW("This function requires CGAL. See http://renaud-detry.net/nuklei/group__install.html");
229 #endif
230  }
231 
232 
233  void BoxROI::buildAABBTree_()
234  {
235 #ifdef NUKLEI_USE_CGAL
236  std::list<Triangle> triangles;
237  std::list<boost::any> t;
238  pushTriangles(t);
239  for (std::list<boost::any>::const_iterator i = t.begin(); i != t.end(); ++i)
240  {
241  triangles.push_back(boost::any_cast<Triangle>(*i));
242  }
243 
244  // constructs AABB tree
245  boost::shared_ptr<Tree> tree(new Tree(triangles.begin(),triangles.end()));
246  tree_ = tree;
247 #else
248  NUKLEI_THROW("This function requires CGAL. See http://renaud-detry.net/nuklei/group__install.html");
249 #endif
250  }
251 
252 #ifdef NUKLEI_USE_CGAL
253  bool intersects(Plane p, const boost::any& t)
254  {
255  NUKLEI_ASSERT(!p.is_degenerate());
256  NUKLEI_ASSERT(!t.empty());
257  boost::shared_ptr<Tree> tree = boost::any_cast< boost::shared_ptr<Tree> >(t);
258  Plane_intersection plane_intersection = tree->any_intersection(p);
259  return static_cast<bool>(plane_intersection);
260  }
261 #endif
262 
263  bool BoxROI::intersectsPlane(const Vector3& p, const Vector3& q, const Vector3& r) const
264  {
265 #ifdef NUKLEI_USE_CGAL
266  Point pp(p.X(), p.Y(), p.Z());
267  Point qq(q.X(), q.Y(), q.Z());
268  Point rr(r.X(), r.Y(), r.Z());
269  Plane plane_query(pp, qq, rr);
270  return intersects(plane_query, tree_);
271 #else
272  NUKLEI_THROW("This function requires CGAL. See http://renaud-detry.net/nuklei/group__install.html");
273 #endif
274  }
275 
276  bool BoxROI::intersectsPlane(const Vector3& p, const Vector3& v) const
277  {
278 #ifdef NUKLEI_USE_CGAL
279  Vector vec(v.X(), v.Y(), v.Z());
280  Point a(p.X(), p.Y(), p.Z());
281  Plane plane_query(a,vec);
282  return intersects(plane_query, tree_);
283 #else
284  NUKLEI_THROW("This function requires CGAL. See http://renaud-detry.net/nuklei/group__install.html");
285 #endif
286  }
287 
288 }
nuklei
Public namespace.
Definition: Color.cpp:9
RegionOfInterest.h
NUKLEI_ASSERT
#define NUKLEI_ASSERT(expression)
Throws an Error if expression is not true.
Definition: Common.h:113
KernelCollection.h
nuklei::la::project
Vector3 project(const Vector3 &x, const Matrix3 &X, const Vector3 &z)
Returns .
Definition: LinearAlgebra.h:434
NUKLEI_THROW
#define NUKLEI_THROW(x)
Throws an Error.
Definition: Common.h:94
© Copyright 2007-2013 Renaud Detry.
Distributed under the terms of the GNU General Public License (GPL).
(See accompanying file LICENSE.txt or copy at http://www.gnu.org/copyleft/gpl.html.)
Revised Sun Sep 13 2020 19:10:06.