DGtal  0.9.3
VoronoiCovarianceMeasureOnDigitalSurface.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file VoronoiCovarianceMeasureOnDigitalSurface.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2014/02/13
23  *
24  * Implementation of inline methods defined in VoronoiCovarianceMeasureOnDigitalSurface.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include "DGtal/topology/CanonicSCellEmbedder.h"
33 #include "DGtal/math/ScalarFunctors.h"
34 #include "DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h"
35 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/ElementaryConvolutionNormalVectorEstimator.h"
36 
37 //////////////////////////////////////////////////////////////////////////////
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 // IMPLEMENTATION of inline methods.
41 ///////////////////////////////////////////////////////////////////////////////
42 
43 ///////////////////////////////////////////////////////////////////////////////
44 // ----------------------- Standard services ------------------------------
45 
46 //-----------------------------------------------------------------------------
47 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
48 inline
49 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
50 ~VoronoiCovarianceMeasureOnDigitalSurface()
51 {
52 }
53 //-----------------------------------------------------------------------------
54 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
55 inline
56 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
57 VoronoiCovarianceMeasureOnDigitalSurface( ConstAlias< Surface > _surface,
58  Surfel2PointEmbedding _surfelEmbedding,
59  Scalar _R, Scalar _r,
60  KernelFunction chi_r,
61  Scalar t, Metric aMetric, bool verbose )
62  : mySurface( _surface ), mySurfelEmbedding( _surfelEmbedding ), myChi( chi_r ),
63  myVCM( _R, _r, aMetric, verbose ), myRadiusTrivial( t )
64 {
65  if ( verbose ) trace.beginBlock( "Computing VCM on digital surface." );
66  const KSpace & ks = mySurface->container().space();
67  std::vector<Point> vectPoints;
68 
69  // Get points.
70  if ( verbose ) trace.beginBlock( "Getting points." );
71  std::set<Point> pointSet;
72  for ( ConstIterator it = mySurface->begin(), itE = mySurface->end(); it != itE; ++it )
73  getPoints( std::inserter( pointSet, pointSet.begin() ), *it );
74  vectPoints.resize( pointSet.size() );
75  std::copy( pointSet.begin(), pointSet.end(), vectPoints.begin() );
76  pointSet.clear();
77  if ( verbose ) trace.endBlock();
78 
79  // Compute Voronoi Covariance Matrix for all points.
80  myVCM.init( vectPoints.begin(), vectPoints.end() );
81 
82  // Compute VCM( chi_r ) for each point.
83  if ( verbose ) trace.beginBlock ( "Integrating VCM( chi_r(p) ) for each point." );
84  int i = 0;
85  // HatPointFunction< Point, Scalar > chi_r( 1.0, r );
86  for ( typename std::vector<Point>::const_iterator it = vectPoints.begin(), itE = vectPoints.end();
87  it != itE; ++it )
88  {
89  if ( verbose ) trace.progressBar( ++i, vectPoints.size() );
90  Point p = *it;
91  MatrixNN measure = myVCM.measure( myChi, p );
92  // On diagonalise le résultat.
93  EigenStructure & evcm = myPt2EigenStructure[ p ];
94  LinearAlgebraTool::getEigenDecomposition( measure, evcm.vectors, evcm.values );
95  }
96  myVCM.clean(); // free some memory.
97  if ( verbose ) trace.endBlock();
98 
99  if ( verbose ) trace.beginBlock ( "Computing average orientation for each surfel." );
100  typedef functors::HatFunction<Scalar> Functor;
101  Functor fct( 1.0, myRadiusTrivial );
102  typedef functors::ElementaryConvolutionNormalVectorEstimator< Surfel, CanonicSCellEmbedder<KSpace> >
103  SurfelFunctor;
104  typedef LocalEstimatorFromSurfelFunctorAdapter< DigitalSurfaceContainer, Metric, SurfelFunctor, Functor>
105  NormalEstimator;
106 
107  CanonicSCellEmbedder<KSpace> canonic_embedder( ks );
108  SurfelFunctor surfelFct( canonic_embedder, 1.0 );
109  NormalEstimator estimator;
110  estimator.attach( *mySurface);
111  estimator.setParams( aMetric, surfelFct, fct , myRadiusTrivial);
112  estimator.init( 1.0, mySurface->begin(), mySurface->end());
113  i = 0;
114  std::vector<Point> pts;
115  int surf_size = mySurface->size();
116  for ( ConstIterator it = mySurface->begin(), itE = mySurface->end(); it != itE; ++it )
117  {
118  if ( verbose ) trace.progressBar(++i, surf_size );
119  Surfel s = *it;
120  Normals & normals = mySurfel2Normals[ s ];
121  // get rough estimation of normal
122  normals.trivialNormal = estimator.eval( it );
123  // get points associated with surfel s
124  getPoints( std::back_inserter( pts ), s );
125  for ( typename std::vector<Point>::const_iterator itPts = pts.begin(), itPtsE = pts.end();
126  itPts != itPtsE; ++itPts )
127  {
128  Point p = *itPts;
129  const EigenStructure& evcm = myPt2EigenStructure[ p ];
130  VectorN n = evcm.vectors.column( Space::dimension-1 );
131  if ( n.dot( normals.trivialNormal ) < 0 ) normals.vcmNormal -= n;
132  else normals.vcmNormal += n;
133  }
134  if ( pts.size() > 1 ) normals.vcmNormal /= pts.size();
135  pts.clear();
136  }
137  if ( verbose ) trace.endBlock();
138 
139  if ( verbose ) trace.endBlock();
140 }
141 
142 //-----------------------------------------------------------------------------
143 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
144 inline
145 DGtal::CountedConstPtrOrConstPtr< typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Surface >
146 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
147 surface() const
148 {
149  return mySurface;
150 }
151 //-----------------------------------------------------------------------------
152 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
153 inline
154 DGtal::Surfel2PointEmbedding
155 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
156 surfelEmbedding() const
157 {
158  return mySurfelEmbedding;
159 }
160 //-----------------------------------------------------------------------------
161 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
162 inline
163 typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Scalar
164 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
165 R() const
166 {
167  return myVCM.R();
168 }
169 //-----------------------------------------------------------------------------
170 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
171 inline
172 typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Scalar
173 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
174 r() const
175 {
176  return myVCM.r();
177 }
178 //-----------------------------------------------------------------------------
179 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
180 inline
181 typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Scalar
182 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
183 radiusTrivial() const
184 {
185  return myRadiusTrivial;
186 }
187 //-----------------------------------------------------------------------------
188 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
189 inline
190 const typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Surfel2Normals&
191 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
192 mapSurfel2Normals() const
193 {
194  return mySurfel2Normals;
195 }
196 //-----------------------------------------------------------------------------
197 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
198 inline
199 const typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Point2EigenStructure&
200 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
201 mapPoint2ChiVCM() const
202 {
203  return myPt2EigenStructure;
204 }
205 
206 //-----------------------------------------------------------------------------
207 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
208 inline
209 bool
210 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
211 getChiVCMEigenvalues( VectorN& values, Surfel s ) const
212 {
213  std::vector<Point> pts;
214  getPoints( std::back_inserter( pts ), s );
215  bool ok = true;
216  int i = 0;
217  values = VectorN(); // Setting values to 0 before averaging.
218  for ( typename std::vector<Point>::const_iterator itPts = pts.begin(), itPtsE = pts.end();
219  itPts != itPtsE; ++itPts, ++i )
220  {
221  Point p = *itPts;
222  typename Point2EigenStructure::const_iterator itEigen = myPt2EigenStructure.find( p );
223  if ( itEigen == myPt2EigenStructure.end() )
224  {
225  ok = false;
226  break;
227  }
228  const EigenStructure& evcm = itEigen->second;
229  values += evcm.values;
230  }
231  if ( i > 1 ) values /= i;
232  return ok;
233 }
234 
235 //-----------------------------------------------------------------------------
236 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
237 inline
238 bool
239 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
240 getChiVCMEigenStructure( VectorN& values, MatrixNN& vectors, Surfel s ) const
241 {
242  std::vector<Point> pts;
243  getPoints( std::back_inserter( pts ), s );
244  bool ok = true;
245  int i = 0;
246  values = VectorN(); // Setting values to 0 before averaging.
247  vectors = MatrixNN(); // Setting values to 0 before averaging.
248  for ( typename std::vector<Point>::const_iterator itPts = pts.begin(), itPtsE = pts.end();
249  itPts != itPtsE; ++itPts, ++i )
250  {
251  Point p = *itPts;
252  typename Point2EigenStructure::const_iterator itEigen = myPt2EigenStructure.find( p );
253  if ( itEigen == myPt2EigenStructure.end() )
254  {
255  ok = false;
256  break;
257  }
258  const EigenStructure& evcm = itEigen->second;
259  values += evcm.values;
260  vectors += evcm.vectors;
261  }
262  if ( i > 1 ) {
263  values /= i;
264  vectors /= i;
265  }
266  return ok;
267 }
268 
269 //-----------------------------------------------------------------------------
270 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
271 template <typename PointOutputIterator>
272 inline
273 PointOutputIterator
274 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
275 getPoints( PointOutputIterator outIt, Surfel s ) const
276 {
277  BOOST_CONCEPT_ASSERT(( boost::OutputIterator< PointOutputIterator, Point > ));
278  const KSpace & ks = mySurface->container().space();
279  Dimension k = ks.sOrthDir( s );
280  switch ( mySurfelEmbedding ) {
281  case Pointels:
282  {
283  typename KSpace::Cells faces = ks.uFaces( ks.unsigns( s ) );
284  for ( typename KSpace::Cells::const_iterator it = faces.begin(), itE = faces.end();
285  it != itE; ++it )
286  {
287  if ( ks.uDim( *it ) == 0 ) // get only pointels (cell of dim 0)
288  *outIt++ = ks.uCoords( *it );
289  }
290  // Dimension i = (k+1)%3;
291  // Dimension j = (i+1)%3;
292  // SCell l1 = ks.sIncident( s, i, true );
293  // SCell l2 = ks.sIncident( s, i, false );
294  // *outIt++ = ks.sCoords( ks.sIncident( l1, j, true ) );
295  // *outIt++ = ks.sCoords( ks.sIncident( l1, j, false ) );
296  // *outIt++ = ks.sCoords( ks.sIncident( l2, j, true ) );
297  // *outIt++ = ks.sCoords( ks.sIncident( l2, j, false ) );
298  }
299  break;
300  case InnerSpel:
301  *outIt++ = ks.sCoords( ks.sDirectIncident( s, k ) );
302  break;
303  case OuterSpel:
304  *outIt++ = ks.sCoords( ks.sIndirectIncident( s, k ) );
305  break;
306  }
307  return outIt;
308 }
309 
310 ///////////////////////////////////////////////////////////////////////////////
311 // Interface - public :
312 
313 /**
314  * Writes/Displays the object on an output stream.
315  * @param out the output stream where the object is written.
316  */
317 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
318 inline
319 void
320 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
321 selfDisplay ( std::ostream & out ) const
322 {
323  out << "[VoronoiCovarianceMeasureOnDigitalSurface"
324  << " #pts=" << myPt2EigenStructure.size()
325  << " #surf=" << mySurfel2Normals.size()
326  << "]";
327 }
328 
329 /**
330  * Checks the validity/consistency of the object.
331  * @return 'true' if the object is valid, 'false' otherwise.
332  */
333 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
334 inline
335 bool
336 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
337 isValid() const
338 {
339  return true;
340 }
341 
342 
343 
344 ///////////////////////////////////////////////////////////////////////////////
345 // Implementation of inline functions //
346 
347 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
348 inline
349 std::ostream&
350 DGtal::operator<< ( std::ostream & out,
351  const VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction> & object )
352 {
353  object.selfDisplay( out );
354  return out;
355 }
356 
357 // //
358 ///////////////////////////////////////////////////////////////////////////////
359 
360 
Generated on Sun Feb 12 2017 10:25:13 for DGtal by   doxygen 1.8.13