KernelCollectionPartialView.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 #ifdef NUKLEI_HAS_PARTIAL_VIEW
8 
9 #ifndef CGAL_EIGEN3_ENABLED
10 # define CGAL_EIGEN3_ENABLED 1
11 #endif
12 
13 #include <CGAL/Simple_cartesian.h>
14 #include <CGAL/AABB_tree.h>
15 #include <CGAL/AABB_traits.h>
16 #include <CGAL/Polyhedron_3.h>
17 #if CGAL_VERSION_NR >= 1040300000
18 #include <CGAL/AABB_face_graph_triangle_primitive.h>
19 #else
20 #include <CGAL/AABB_polyhedron_triangle_primitive.h>
21 #endif
22 #include <CGAL/property_map.h>
23 #include <CGAL/IO/read_off_points.h>
24 #include <CGAL/IO/read_xyz_points.h>
25 #include <CGAL/IO/write_xyz_points.h>
26 #include <CGAL/IO/Polyhedron_iostream.h>
27 
28 #include <CGAL/squared_distance_3.h>
29 #include <CGAL/Point_with_normal_3.h>
30 
31 #endif
32 
33 #include <nuklei/KernelCollection.h>
34 #include <nuklei/ObservationIO.h>
35 
36 #include <boost/shared_ptr.hpp>
37 
38 #include <nuklei/KernelCollection.h>
39 
40 #ifdef NUKLEI_HAS_PARTIAL_VIEW
41 
42 typedef CGAL::Simple_cartesian<double> K;
43 typedef K::Point_3 Point;
44 typedef K::Plane_3 Plane;
45 typedef K::Vector_3 Vector;
46 typedef K::Segment_3 Segment;
47 typedef K::Line_3 Line;
48 typedef CGAL::Polyhedron_3<K> Polyhedron;
49 #if CGAL_VERSION_NR >= 1040300000
50 typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
51 #else
52 typedef CGAL::AABB_polyhedron_triangle_primitive<K,Polyhedron> Primitive;
53 #endif
54 typedef CGAL::AABB_traits<K, Primitive> Traits;
55 typedef CGAL::AABB_tree<Traits> Tree;
56 typedef Tree::Object_and_primitive_id Object_and_primitive_id;
57 typedef Tree::Primitive_id Primitive_id;
58 
59 typedef CGAL::Point_with_normal_3<K> Point_with_normal;
60 typedef std::vector<Point_with_normal> PointList;
61 
62 #endif
63 
64 namespace nuklei {
65 
66  typedef std::vector< std::pair< Vector3, std::vector<int> > > viewcache_t;
67 
68 #ifdef NUKLEI_HAS_PARTIAL_VIEW
69 
70 
71  inline bool isVisible(const Vector3& wtarget,
72  const Vector3& wnormal,
73  const Vector3& wcamera,
74  const Polyhedron& poly,
75  const Tree& tree,
76  const coord_t& tolerance)
77  {
78  if (std::fabs(wnormal.SquaredLength()-1) < FLOATTOL)
79  {
80  Vector3 ctot = la::normalized(wtarget-wcamera);
81  double dot = wnormal.Dot(ctot);
82  if (std::acos(std::fabs(dot)) > (80./180*M_PI))
83  {
84  return false;
85  }
86  }
87 
88  Point camera(wcamera.X(), wcamera.Y(), wcamera.Z());
89  Point target(wtarget.X(), wtarget.Y(), wtarget.Z());
90 
91  //double d = dot();
92  Segment segment_query(camera,target);
93  bool visible = true;
94 
95  std::list<Object_and_primitive_id> intersections;
96  tree.all_intersections(segment_query, std::back_inserter(intersections));
97 
98  int self_intersect = 0;
99  for (std::list<Object_and_primitive_id>::const_iterator intersection = intersections.begin();
100  intersection != intersections.end(); ++intersection)
101  {
102  Object_and_primitive_id op = *intersection;
103  CGAL::Object object = op.first;
104  Point intersectionPoint;
105  if(CGAL::assign(intersectionPoint,object))
106  {
107  double sq2 = squared_distance(intersectionPoint, target);
108  if (sq2 < tolerance*tolerance) self_intersect++;
109  }
110  else
111  {
112  // This is very unlikely (p<<0.001)
113  //NUKLEI_THROW("intersection is not a point!");
114  }
115  }
116 
117  int num_inter = intersections.size();
118  if (self_intersect > 1)
119  NUKLEI_LOG("Number of self-intersections (" << self_intersect << ") greater than 1");
120 
121  if (num_inter == 0)
122  visible = true;
123  else
124  {
125  if (num_inter - self_intersect > 0)
126  visible = false;
127  else
128  visible = true;
129  }
130  if (num_inter > 1 && visible)
131  NUKLEI_INFO("Strange visibility at" << target);
132 
133  return visible;
134 
135  }
136 
137  inline bool isVisible(const Vector3& wtarget,
138  const Vector3& wcamera,
139  const Polyhedron& poly,
140  const Tree& tree,
141  const coord_t& tolerance)
142  {
143  return isVisible(wtarget,
144  Vector3::ZERO,
145  wcamera,
146  poly,
147  tree,
148  tolerance);
149  }
150 
151 #endif
152 
153  bool KernelCollection::isVisibleFrom(const Vector3& p, const Vector3& viewpoint,
154  const coord_t& tolerance) const
155  {
156  NUKLEI_TRACE_BEGIN();
157 #ifdef NUKLEI_HAS_PARTIAL_VIEW
158  if (!deco_.has_key(MESH_KEY))
159  NUKLEI_THROW("Undefined mesh. Call buildMesh() first.");
160  if (!deco_.has_key(AABBTREE_KEY))
161  NUKLEI_THROW("Undefined AABB tree. Call buildMesh() first.");
162 
163  return isVisible(p, viewpoint,
164  *deco_.get< boost::shared_ptr<Polyhedron> >(MESH_KEY),
165  *deco_.get< boost::shared_ptr<Tree> >(AABBTREE_KEY),
166  tolerance);
167 #else
168  NUKLEI_THROW("This function requires the partial view build of Nuklei. See http://renaud-detry.net/nuklei/group__install.html");
169 #endif
170  NUKLEI_TRACE_END();
171  }
172 
173  bool KernelCollection::isVisibleFrom(const kernel::r3xs2p& p, const Vector3& viewpoint,
174  const coord_t& tolerance) const
175  {
176  NUKLEI_TRACE_BEGIN();
177 #ifdef NUKLEI_HAS_PARTIAL_VIEW
178  if (!deco_.has_key(MESH_KEY))
179  NUKLEI_THROW("Undefined mesh. Call buildMesh() first.");
180  if (!deco_.has_key(AABBTREE_KEY))
181  NUKLEI_THROW("Undefined AABB tree. Call buildMesh() first.");
182 
183  return isVisible(p.loc_, p.dir_, viewpoint,
184  *deco_.get< boost::shared_ptr<Polyhedron> >(MESH_KEY),
185  *deco_.get< boost::shared_ptr<Tree> >(AABBTREE_KEY),
186  tolerance);
187 #else
188  NUKLEI_THROW("This function requires the partial view build of Nuklei. See http://renaud-detry.net/nuklei/group__install.html");
189 #endif
190  NUKLEI_TRACE_END();
191  }
192 
193  template<typename C>
194  C KernelCollection::partialView(const Vector3& viewpoint,
195  const coord_t& tolerance,
196  const bool useViewcache,
197  const bool useRayToSurfacenormalAngle) const
198  {
199  NUKLEI_TRACE_BEGIN();
200 
201  C index_collection;
202 
203  if (!useViewcache)
204  {
205 #ifdef NUKLEI_HAS_PARTIAL_VIEW
206  if (!deco_.has_key(MESH_KEY))
207  NUKLEI_THROW("Undefined mesh. Call buildMesh() first.");
208  if (!deco_.has_key(AABBTREE_KEY))
209  NUKLEI_THROW("Undefined AABB tree. Call buildMesh() first.");
210 
211  const Polyhedron& poly = *deco_.get< boost::shared_ptr<Polyhedron> >(MESH_KEY);
212  const Tree& tree = *deco_.get< boost::shared_ptr<Tree> >(AABBTREE_KEY);
213 
214  for (const_iterator v = begin(); v != end(); ++v)
215  {
216  Vector3 p = v->getLoc();
217  Vector3 normal = Vector3::ZERO;
218  if (useRayToSurfacenormalAngle)
219  normal = kernel::r3xs2p(*v).dir_;
220  if (isVisible(p, normal, viewpoint, poly, tree, tolerance))
221  index_collection.push_back(std::distance(begin(), v));
222  }
223 #else
224  NUKLEI_THROW("This function requires the partial view build of Nuklei. See http://renaud-detry.net/nuklei/group__install.html");
225 #endif
226  }
227  else
228  {
229  if (!deco_.has_key(VIEWCACHE_KEY))
230  NUKLEI_THROW("Undefined view cache. Call buildPartialViewCache() first.");
231 
232  const viewcache_t &viewIndex = *deco_.get< boost::shared_ptr<viewcache_t> >(VIEWCACHE_KEY);
233 
234  double minDist = 1e6;
235  viewcache_t::const_iterator closest = as_const(viewIndex).begin();
236  for (viewcache_t::const_iterator oo = as_const(viewIndex).begin(); oo != viewIndex.end(); ++oo)
237  {
238  double tmp = (viewpoint-oo->first).SquaredLength();
239  if (tmp < minDist*minDist)
240  {
241  minDist = tmp;
242  closest = oo;
243  }
244  }
245  for (std::vector<int>::const_iterator i = closest->second.begin(); i != closest->second.end(); ++i)
246  index_collection.push_back(*i);
247  }
248  return index_collection;
249  NUKLEI_TRACE_END();
250  }
251 
252 
253  std::vector<int> KernelCollection::partialView(const Vector3& viewpoint,
254  const coord_t& tolerance,
255  const bool useViewcache,
256  const bool useRayToSurfacenormalAngle) const
257  {
258  NUKLEI_TRACE_BEGIN();
259 #ifdef NUKLEI_HAS_PARTIAL_VIEW
260  return partialView< std::vector<int> >(viewpoint, tolerance, useViewcache, useRayToSurfacenormalAngle);
261 #else
262  NUKLEI_THROW("This function requires the partial view build of Nuklei. See http://renaud-detry.net/nuklei/group__install.html");
263 #endif
264  NUKLEI_TRACE_END();
265  }
266 
267 
268  KernelCollection::const_partialview_iterator KernelCollection::partialViewBegin(const Vector3& viewpoint,
269  const coord_t& tolerance,
270  const bool useViewcache,
271  const bool useRayToSurfacenormalAngle) const
272  {
273  NUKLEI_TRACE_BEGIN();
274 #ifdef NUKLEI_HAS_PARTIAL_VIEW
275 
276  typedef const_partialview_iterator::index_container
277  index_container;
278  typedef const_partialview_iterator::index_container_ptr
279  index_container_ptr;
280  typedef const_partialview_iterator::index_t
281  index_t;
282 
283  index_container_ptr index_collection(new index_container);
284  *index_collection = partialView< index_container >(viewpoint, tolerance, useViewcache, useRayToSurfacenormalAngle);
285 
286  return const_partialview_iterator(begin(), index_collection);
287 #else
288  NUKLEI_THROW("This function requires the partial view build of Nuklei. See http://renaud-detry.net/nuklei/group__install.html");
289 #endif
290  NUKLEI_TRACE_END();
291  }
292 
293  void KernelCollection::buildPartialViewCache(const double meshTol, const bool useRayToSurfacenormalAngle)
294  {
295  NUKLEI_TRACE_BEGIN();
296 #ifdef NUKLEI_HAS_PARTIAL_VIEW
297  Vector3 mean = as_const(*this).moments()->getLoc();
298  double stdev = as_const(*this).moments()->getLocH();
299 
300  std::vector< Vector3 > keys;
301  for (int o = 0; o < 10000; ++o)
302  {
303  Vector3 key = Random::uniformDirection3d();
304  keys.push_back(key);
305  }
306 
307  {
308  double minDist2 = .15;
309  std::vector< Vector3 > tmp = keys;
310  keys.clear();
311  for (std::vector<Vector3>::const_iterator
312  i = tmp.begin();
313  i != tmp.end(); ++i)
314  {
315  bool cnt = false;
316  Vector3 iLoc = (*i);
317  for (std::vector<Vector3>::const_iterator
318  j = keys.begin();
319  j != keys.end(); ++j)
320  {
321  if ( (iLoc - *j).SquaredLength() < minDist2*minDist2 )
322  {
323  cnt = true;
324  break;
325  }
326  }
327  if (!cnt)
328  {
329  keys.push_back(iLoc);
330  }
331  }
332  }
333 
334  boost::shared_ptr<viewcache_t> viewIndex(new viewcache_t());
335 
336  for (unsigned int o = 0; o < keys.size(); ++o)
337  {
338  Vector3 key = keys.at(o);
339  Vector3 vp = mean + key*stdev*20;
340  std::vector<int> vi = as_const(*this).partialView(vp, meshTol, false, useRayToSurfacenormalAngle);
341 #if 0
342  // debug - delete when code is considered stable
343  KernelCollection v;
344  for (std::vector<int>::iterator i = vi.begin(); i != vi.end(); ++i)
345  v.add(as_const(*this).at(*i));
346  v.computeKernelStatistics();
347  kernel::base::ptr k = v.randomKernel().create();
348  k->setLoc(vp);
349  v.add(*k);
350  v.computeKernelStatistics();
351 #endif
352  viewIndex->push_back(std::make_pair(key, vi));
353  }
354 
355  if (deco_.has_key(VIEWCACHE_KEY)) deco_.erase(VIEWCACHE_KEY);
356  deco_.insert(VIEWCACHE_KEY, viewIndex);
357 
358 #if 0
359  // debug - delete when code is considered stable
360 
361  for (viewcache_t::iterator o = viewIndex->begin(); o != viewIndex->end(); ++o)
362  {
363  double d = 1e6;
364  KernelCollection near;
365 
366  for (viewcache_t::iterator oo = viewIndex->begin(); oo != viewIndex->end(); ++oo)
367  {
368  if (oo == o) continue;
369  double dd = (o->first-oo->first).Length();
370  if (dd < d) {
371  d = dd;
372  near.clear();
373  for (std::vector<int>::iterator j = oo->second.begin(); j != oo->second.end(); ++j)
374  near.add(as_const(*this).at(*j));
375  near.computeKernelStatistics();
376  if (near.size() == 0) continue;
377  kernel::base::ptr k = near.randomKernel().create();
378  k->setLoc(mean + oo->first*stdev*20);
379  near.add(*k);
380  near.computeKernelStatistics();
381  }
382  }
383  for (KernelCollection::iterator i = near.begin(); i != near.end(); ++i)
384  {
385  ColorDescriptor cd;
386  cd.setColor(RGBColor(1, 0, 0));
387  i->setDescriptor(cd);
388  }
389  for (std::vector<int>::iterator j = o->second.begin(); j != o->second.end(); ++j)
390  {
391  near.add(as_const(*this).at(*j));
392  near.back().setLoc(near.back().getLoc()+Vector3(0.001, 0, 0));
393  }
394  near.computeKernelStatistics();
395  if (near.size() == 0) continue;
396  kernel::base::ptr k = near.randomKernel().create();
397  k->setLoc(mean + o->first*stdev*20);
398  near.add(*k);
399  near.computeKernelStatistics();
400  writeObservations("/tmp/v/" + stringify(std::distance(viewIndex->begin(), o)), near, Observation::SERIAL);
401  }
402 #endif
403 
404 #else
405  NUKLEI_THROW("This function requires the partial view build of Nuklei. See http://renaud-detry.net/nuklei/group__install.html");
406 #endif
407  NUKLEI_TRACE_END();
408 
409  }
410 
411 }
412 
nuklei::RGBColor
Definition: Color.h:56
nuklei::KernelCollection::add
void add(const kernel::base &f)
Adds a copy of f.
Definition: KernelCollection.cpp:64
nuklei
Public namespace.
Definition: Color.cpp:9
nuklei::KernelCollection::const_iterator
Container::const_iterator const_iterator
KernelCollection iterator.
Definition: KernelCollection.h:165
nuklei::KernelCollection::partialViewBegin
const_partialview_iterator partialViewBegin(const Vector3 &viewpoint, const coord_t &tolerance=FLOATTOL, const bool useViewcache=false, const bool useRayToSurfacenormalAngle=false) const
Assuming that the points in this collection form the surface of an object, this function returns an i...
Definition: KernelCollectionPartialView.cpp:268
nuklei::writeObservations
void writeObservations(ObservationWriter &w, const KernelCollection &kc)
Writes the content of kc using the provided writer w.
Definition: ObservationIO.cpp:454
nuklei::stringify
std::string stringify(const T &x, int precision=-1, int width=0)
Converts an object to a std::string using operator<<.
Definition: Common.h:333
nuklei::KernelCollection::buildPartialViewCache
void buildPartialViewCache(const double meshTol, const bool useRayToSurfacenormalAngle=false)
Builds set of partial views of the object. See Intermediary Results.
Definition: KernelCollectionPartialView.cpp:293
nuklei::kernel::base::create
virtual base::ptr create() const =0
Create a new kernel of the same type.
nuklei::KernelCollection
This class acts as a vector-like container for kernels. It also provides methods related to kernel de...
Definition: KernelCollection.h:153
nuklei::kernel::r3xs2p
r3xs2_base< groupS::s2p > r3xs2p
Definition: Kernel.h:623
nuklei::ColorDescriptor
Definition: Descriptor.h:68
nuklei::KernelCollection::back
Container::reference back()
Returns the kernel at index size()-1.
Definition: KernelCollection.h:186
nuklei::kernel::base::ptr
NUKLEI_UNIQUE_PTR< kernel::base > ptr
NUKLEI_UNIQUE_PTR for kernel::base.
Definition: Kernel.h:50
nuklei::kernel::r3xs2_base
Definition: Kernel.h:505
ObservationIO.h
nuklei::coord_t
double coord_t
Type for point coordinates.
Definition: Definitions.h:25
nuklei::KernelCollection::end
Container::iterator end()
Returns an iterator pointing to the last kernel.
Definition: KernelCollection.h:207
nuklei::KernelCollection::randomKernel
const kernel::base & randomKernel() const
Returns a kernel from the collection.
Definition: KernelCollection.cpp:301
nuklei::KernelCollection::mean
kernel::base::ptr mean() const
Returns a kernel holding the mean position and orientation of the data.
Definition: KernelCollection.cpp:201
nuklei::KernelCollection::begin
Container::iterator begin()
Returns an iterator pointing to the first kernel.
Definition: KernelCollection.h:201
nuklei::KernelCollection::const_partialview_iterator
nuklei_trsl::reorder_iterator< const_iterator > const_partialview_iterator
Partial View Iterator type.
Definition: KernelCollection.h:491
nuklei::KernelCollection::isVisibleFrom
bool isVisibleFrom(const Vector3 &p, const Vector3 &viewpoint, const coord_t &tolerance=FLOATTOL) const
Assuming that the points in this collection form the surface of an object, this function computes whe...
Definition: KernelCollectionPartialView.cpp:153
nuklei::KernelCollection::size
Container::size_type size() const
Returns the number of kernels.
Definition: KernelCollection.h:193
nuklei::KernelCollection::computeKernelStatistics
void computeKernelStatistics()
Computes the sum of all kernel weights (total weight), and the maximum kernel cut point.
Definition: KernelCollection.cpp:139
nuklei::KernelCollection::partialView
std::vector< int > partialView(const Vector3 &viewpoint, const coord_t &tolerance=FLOATTOL, const bool useViewcache=false, const bool useRayToSurfacenormalAngle=false) const
Assuming that the points in this collection form the surface of an object, this function returns the ...
Definition: KernelCollectionPartialView.cpp:253
nuklei::KernelCollection::iterator
Container::iterator iterator
KernelCollection iterator.
Definition: KernelCollection.h:163
nuklei::KernelCollection::clear
void clear()
Resets the class to its initial state.
Definition: KernelCollection.cpp:55
nuklei::Random::uniformDirection3d
static Vector3 uniformDirection3d()
This function returns a random direction vector in three dimensions.
Definition: Random.cpp:320
KernelCollection.h
nuklei::as_const
T const & as_const(T const &t)
Provides a const reference to an object.
Definition: Common.h:197
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.