DGtal  0.9.2
SphericalHoughNormalVectorEstimator.h
1 
17 #pragma once
18 
34 #if defined(SphericalHoughNormalVectorEstimator_RECURSES)
35 #error Recursive header files inclusion detected in SphericalHoughNormalVectorEstimator.h
36 #else // defined(SphericalHoughNormalVectorEstimator_RECURSES)
37 
38 #define SphericalHoughNormalVectorEstimator_RECURSES
39 
40 #if !defined SphericalHoughNormalVectorEstimator_h
41 
42 #define SphericalHoughNormalVectorEstimator_h
43 
45 // Inclusions
46 #include <iostream>
47 #include <DGtal/base/Common.h>
48 #include <DGtal/topology/SCellsFunctors.h>
49 #include <vector>
50 #include "DGtal/geometry/tools/SphericalAccumulator.h"
51 #include <random>
52 #include "DGtal/math/linalg/SimpleMatrix.h"
54 
55 namespace DGtal
56 {
57  namespace functors
58  {
60  // template class SphericalHoughNormalVectorEstimator
83  template <typename TSurfel, typename TEmbedder>
85  {
86  public:
87 
88  typedef TSurfel Surfel;
89  typedef TEmbedder SCellEmbedder;
90  typedef typename SCellEmbedder::RealPoint RealPoint;
91  typedef RealPoint Quantity;
93 
105  const double h,
106  const double minimalAspectRatio = 0.001,
107  const unsigned int nbTrials = 100,
108  const unsigned int accumulatorSize = 10,
109  const unsigned int nbAccumulators = 5) :
110  myEmbedder(&anEmbedder),myH(h), myAspectRatio(minimalAspectRatio),
111  myNbTrials( nbTrials), mySize(accumulatorSize) , myNbAccumulators(nbAccumulators)
112  {
114 
115  //We precompute the random rotations and accumulators
116  for(auto i=0; i < myNbAccumulators; ++i)
117  {
118  Matrix m = randomRotation();
119  myAccumulators.push_back( accum );
120  myRotations.push_back( m );
121  myInverseRotations.push_back( m.inverse() );
122  }
123  }
124 
129 
130 
138  void pushSurfel(const Surfel & aSurf,
139  const double aDistance)
140  {
141  BOOST_VERIFY(aDistance == aDistance);
142  RealPoint p = myH * ( myEmbedder->operator()(aSurf) );
143  myPoints.push_back(p);
144  }
145 
152  Quantity eval( )
153  {
154  std::default_random_engine generator;
155  std::uniform_int_distribution<int> distribution(0, myPoints.size() - 1 );
156  double aspect;
157 
158  for(auto t = 0; t < myNbTrials ; ++t)
159  {
160  unsigned int i,j,k;
161 
162  //We pick 3 distinct point indices.
163  i = distribution(generator);
164  j = distribution(generator);
165  while ( (j = distribution(generator)) == i);
166  while (( (k = distribution(generator)) == i) || (k == j) );
167 
168  RealPoint vector = getNormal(i,j,k,aspect);
169  if ((vector.norm() > 0.00001) && (aspect > myAspectRatio))
170  {
171  //we have an admissible triangle, we push both normal vectors
172  for(auto acc=0; acc < myNbAccumulators; ++acc)
173  {
174  RealPoint shifted = myRotations[acc]*vector;
175  myAccumulators[acc].addDirection( shifted );
176  myAccumulators[acc].addDirection( -shifted );
177  }
178  }
179  }
180  //We return the max bin orientation summing up all accumulators vote
181  typename SphericalAccumulator<RealPoint>::Size posPhi,posTheta;
182  RealPoint vote;
183  for(auto acc=0; acc < myNbAccumulators; ++acc)
184  {
185  myAccumulators[acc].maxCountBin(posPhi, posTheta);
186  RealPoint dir = myInverseRotations[acc]*myAccumulators[acc].representativeDirection(posPhi, posTheta).getNormalized() ;
187 
188  //We only consider z-oriented normals (since we pushed vector and -vector)
189  if ( dir.dot(RealPoint(0,0,1)) > 0.0 )
190  vote += dir;
191  else
192  vote += -dir;
193  }
194  return vote.getNormalized();
195  }
196 
201  void reset()
202  {
203  myPoints.clear();
204  //accumulators cleanup
205  for(auto i=0; i < myNbAccumulators; ++i)
206  myAccumulators[i].clear();
207  }
208 
209  private:
210 
214  Matrix randomRotation() const
215  {
216  const double theta = (rand()+0.f)/RAND_MAX * 2* M_PI;
217  const double phi = (rand()+0.f)/RAND_MAX * 2* M_PI;
218  const double psi = (rand()+0.f)/RAND_MAX * 2* M_PI;
219  Matrix Rt;
220  Rt.setComponent(0,0,1);
221  Rt.setComponent(1,0,0);
222  Rt.setComponent(2,0,0);
223  Rt.setComponent(0,1,0);
224  Rt.setComponent(1,1,cos(theta));
225  Rt.setComponent(2,1,-sin(theta));
226  Rt.setComponent(0,2,0);
227  Rt.setComponent(1,2,sin(theta));
228  Rt.setComponent(2,2,cos(theta));
229 
230  Matrix Rph;
231  Rph.setComponent(0,0,cos(phi));
232  Rph.setComponent(1,0,0);
233  Rph.setComponent(2,0,sin(phi));
234  Rph.setComponent(0,1,0);
235  Rph.setComponent(1,1,1);
236  Rph.setComponent(2,1,0);
237  Rph.setComponent(0,2,-sin(phi));
238  Rph.setComponent(1,2,0);
239  Rph.setComponent(2,2,cos(phi));
240 
241  Matrix Rps;
242  Rps.setComponent(0,0,cos(psi));
243  Rps.setComponent(1,0,-sin(psi));
244  Rps.setComponent(2,0,0);
245  Rps.setComponent(0,1,sin(psi));
246  Rps.setComponent(1,1,cos(psi));
247  Rps.setComponent(2,1,0);
248  Rps.setComponent(0,2,0);
249  Rps.setComponent(1,2,0);
250  Rps.setComponent(2,2,1);
251 
252  return Rt*Rph*Rps;
253  }
254 
267  Quantity getNormal(const unsigned int i,
268  const unsigned int j,
269  const unsigned int k,
270  double &aspect) const
271  {
272  ASSERT( i < myPoints.size());
273  ASSERT( j < myPoints.size());
274  ASSERT( k < myPoints.size());
275 
276  const RealPoint v = myPoints[i] - myPoints[j];
277  const RealPoint u = myPoints[i] - myPoints[k];
278  const RealPoint w = myPoints[j] - myPoints[k];
279 
280  //aspect ratio
281  const double a = u.norm() , b = v.norm();
282  const double c = w.norm();
283  const double s = (a+b+c)/2.0;
284  aspect = a*b*c/(8.0*(s-a)*(s-b)*(s-c));
285 
286  return v.crossProduct(u);
287  }
288 
290  const SCellEmbedder * myEmbedder;
291 
293  const double myH;
294 
296  const double myAspectRatio;
297 
299  const unsigned int myNbTrials;
300 
302  const unsigned int mySize;
303 
305  const unsigned int myNbAccumulators;
306 
308  std::vector<RealPoint> myPoints;
309 
311  std::vector< SphericalAccumulator<RealPoint> > myAccumulators;
312 
314  std::vector< Matrix > myRotations;
315 
317  std::vector< Matrix > myInverseRotations;
318 
319 
320  }; // end of class SphericalHoughNormalVectorEstimator
321  } // namespace functors
322 } // namespace DGtal
323 
324 
325 // //
327 
328 #endif // !defined SphericalHoughNormalVectorEstimator_h
329 
330 #undef SphericalHoughNormalVectorEstimator_RECURSES
331 #endif // else defined(SphericalHoughNormalVectorEstimator_RECURSES)
Aim: This functor estimates normal vector for a collection of surfels using spherical accumulator bas...
SimpleMatrix< Component, TM, TN > inverse() const
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:186
std::vector< Matrix > myInverseRotations
Random inverse rotations.
Quantity getNormal(const unsigned int i, const unsigned int j, const unsigned int k, double &aspect) const
void setComponent(const DGtal::Dimension i, const DGtal::Dimension j, const Component &aValue)
void pushSurfel(const Surfel &aSurf, const double aDistance)
DGtal is the top-level namespace which contains all DGtal functions and types.
std::vector< SphericalAccumulator< RealPoint > > myAccumulators
Spherical Accumulators.
const unsigned int myNbTrials
Number of trials in the neignborhood.
const SCellEmbedder * myEmbedder
Alias of the geometrical embedder.
const double myAspectRatio
Minimal aspect ratio (norm of the cross-product) to consider a given triangle.
const unsigned int myNbAccumulators
Number of randomly shifted spherical accumulators to consider.
SphericalHoughNormalVectorEstimator(ConstAlias< SCellEmbedder > anEmbedder, const double h, const double minimalAspectRatio=0.001, const unsigned int nbTrials=100, const unsigned int accumulatorSize=10, const unsigned int nbAccumulators=5)
std::vector< RealPoint > myPoints
vector of embedded surfels