DGtal 1.3.0
Loading...
Searching...
No Matches
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
55namespace 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;
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 {
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
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:76
void setComponent(const DGtal::Dimension i, const DGtal::Dimension j, const Component &aValue)
Aim: implements an accumulator (as histograms for 1D scalars) adapted to spherical point samples.
size_t Size
Type to represent bin indexes.
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.
DGtal is the top-level namespace which contains all DGtal functions and types.