DGtal  1.5.beta
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)
38 #define SphericalHoughNormalVectorEstimator_RECURSES
39 
40 #if !defined SphericalHoughNormalVectorEstimator_h
42 #define SphericalHoughNormalVectorEstimator_h
43 
45 // Inclusions
46 #include <cmath>
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;
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 = 0u; 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 
153  {
154  std::default_random_engine generator;
155  std::uniform_int_distribution<int> distribution(0,
156  static_cast<int>(myPoints.size()) - 1 );
157  double aspect;
158 
159  for(auto t = 0u; t < myNbTrials ; ++t)
160  {
161  unsigned int i,j,k;
162 
163  //We pick 3 distinct point indices.
164  i = distribution(generator);
165  j = distribution(generator);
166  while ( (j = distribution(generator)) == i);
167  while (( (k = distribution(generator)) == i) || (k == j) );
168 
169  RealPoint vector = getNormal(i,j,k,aspect);
170  if ((vector.norm() > 0.00001) && (aspect > myAspectRatio))
171  {
172  //we have an admissible triangle, we push both normal vectors
173  for(auto acc = 0u; acc < myNbAccumulators; ++acc)
174  {
175  RealPoint shifted = myRotations[acc]*vector;
176  myAccumulators[acc].addDirection( shifted );
177  myAccumulators[acc].addDirection( -shifted );
178  }
179  }
180  }
181  //We return the max bin orientation summing up all accumulators vote
182  typename SphericalAccumulator<RealPoint>::Size posPhi,posTheta;
183  RealPoint vote;
184  for(auto acc = 0u; acc < myNbAccumulators; ++acc)
185  {
186  myAccumulators[acc].maxCountBin(posPhi, posTheta);
187  RealPoint dir = myInverseRotations[acc]*myAccumulators[acc].representativeDirection(posPhi, posTheta).getNormalized() ;
188 
189  //We only consider z-oriented normals (since we pushed vector and -vector)
190  if ( dir.dot(RealPoint(0,0,1)) > 0.0 )
191  vote += dir;
192  else
193  vote += -dir;
194  }
195  return vote.getNormalized();
196  }
197 
202  void reset()
203  {
204  myPoints.clear();
205  //accumulators cleanup
206  for(auto i = 0u; i < myNbAccumulators; ++i)
207  myAccumulators[i].clear();
208  }
209 
210  private:
211 
216  {
217  const double theta = (rand()+0.f)/RAND_MAX * 2* M_PI;
218  const double phi = (rand()+0.f)/RAND_MAX * 2* M_PI;
219  const double psi = (rand()+0.f)/RAND_MAX * 2* M_PI;
220  Matrix Rt;
221  Rt.setComponent(0,0,1);
222  Rt.setComponent(1,0,0);
223  Rt.setComponent(2,0,0);
224  Rt.setComponent(0,1,0);
225  Rt.setComponent(1,1,cos(theta));
226  Rt.setComponent(2,1,-sin(theta));
227  Rt.setComponent(0,2,0);
228  Rt.setComponent(1,2,sin(theta));
229  Rt.setComponent(2,2,cos(theta));
230 
231  Matrix Rph;
232  Rph.setComponent(0,0,cos(phi));
233  Rph.setComponent(1,0,0);
234  Rph.setComponent(2,0,sin(phi));
235  Rph.setComponent(0,1,0);
236  Rph.setComponent(1,1,1);
237  Rph.setComponent(2,1,0);
238  Rph.setComponent(0,2,-sin(phi));
239  Rph.setComponent(1,2,0);
240  Rph.setComponent(2,2,cos(phi));
241 
242  Matrix Rps;
243  Rps.setComponent(0,0,cos(psi));
244  Rps.setComponent(1,0,-sin(psi));
245  Rps.setComponent(2,0,0);
246  Rps.setComponent(0,1,sin(psi));
247  Rps.setComponent(1,1,cos(psi));
248  Rps.setComponent(2,1,0);
249  Rps.setComponent(0,2,0);
250  Rps.setComponent(1,2,0);
251  Rps.setComponent(2,2,1);
252 
253  return Rt*Rph*Rps;
254  }
255 
268  Quantity getNormal(const unsigned int i,
269  const unsigned int j,
270  const unsigned int k,
271  double &aspect) const
272  {
273  ASSERT( i < myPoints.size());
274  ASSERT( j < myPoints.size());
275  ASSERT( k < myPoints.size());
276 
277  const RealPoint v = myPoints[i] - myPoints[j];
278  const RealPoint u = myPoints[i] - myPoints[k];
279  const RealPoint w = myPoints[j] - myPoints[k];
280 
281  //aspect ratio
282  const double a = u.norm() , b = v.norm();
283  const double c = w.norm();
284  const double s = (a+b+c)/2.0;
285  double denom = (8.0*(s-a)*(s-b)*(s-c));
286  if ( std::abs( denom ) <= std::abs( denom ) * 1e-6 )
287  aspect = 0.;
288  else
289  aspect = a*b*c / denom;
290 
291  return v.crossProduct(u);
292  }
293 
296 
298  const double myH;
299 
301  const double myAspectRatio;
302 
304  const unsigned int myNbTrials;
305 
307  const unsigned int mySize;
308 
310  const unsigned int myNbAccumulators;
311 
313  std::vector<RealPoint> myPoints;
314 
316  std::vector< SphericalAccumulator<RealPoint> > myAccumulators;
317 
319  std::vector< Matrix > myRotations;
320 
322  std::vector< Matrix > myInverseRotations;
323 
324 
325  }; // end of class SphericalHoughNormalVectorEstimator
326  } // namespace functors
327 } // namespace DGtal
328 
329 
330 // //
332 
333 #endif // !defined SphericalHoughNormalVectorEstimator_h
334 
335 #undef SphericalHoughNormalVectorEstimator_RECURSES
336 #endif // else defined(SphericalHoughNormalVectorEstimator_RECURSES)
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
void setComponent(const DGtal::Dimension i, const DGtal::Dimension j, const Component &aValue)
Aim: This functor estimates normal vector for a collection of surfels using spherical accumulator bas...
std::vector< RealPoint > myPoints
vector of embedded surfels
void pushSurfel(const Surfel &aSurf, const double aDistance)
std::vector< SphericalAccumulator< RealPoint > > myAccumulators
Spherical Accumulators.
std::vector< Matrix > myInverseRotations
Random inverse rotations.
const SCellEmbedder * myEmbedder
Alias of the geometrical embedder.
Quantity getNormal(const unsigned int i, const unsigned int j, const unsigned int k, double &aspect) const
const unsigned int myNbAccumulators
Number of randomly shifted spherical accumulators to consider.
const unsigned int myNbTrials
Number of trials in the neignborhood.
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)
const double myAspectRatio
Minimal aspect ratio (norm of the cross-product) to consider a given triangle.
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi(const Face f)
DGtal is the top-level namespace which contains all DGtal functions and types.
PointVector< 3, double > RealPoint