DGtal  1.2.0
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 #include "DGtal/geometry/volumes/distance/LpMetric.h"
37 
38 //////////////////////////////////////////////////////////////////////////////
39 
40 ///////////////////////////////////////////////////////////////////////////////
41 // IMPLEMENTATION of inline methods.
42 ///////////////////////////////////////////////////////////////////////////////
43 
44 ///////////////////////////////////////////////////////////////////////////////
45 // ----------------------- Standard services ------------------------------
46 
47 //-----------------------------------------------------------------------------
48 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
49 inline
50 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
51 ~VoronoiCovarianceMeasureOnDigitalSurface()
52 {
53 }
54 //-----------------------------------------------------------------------------
55 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
56 inline
57 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
58 VoronoiCovarianceMeasureOnDigitalSurface( ConstAlias< Surface > _surface,
59  Surfel2PointEmbedding _surfelEmbedding,
60  Scalar _R, Scalar _r,
61  KernelFunction chi_r,
62  Scalar t, Metric aMetric, bool verbose )
63  : mySurface( _surface ), mySurfelEmbedding( _surfelEmbedding ), myChi( chi_r ),
64  myVCM( _R, _r, aMetric, verbose ), myRadiusTrivial( t )
65 {
66  if ( verbose ) trace.beginBlock( "Computing VCM on digital surface." );
67  const KSpace & ks = mySurface->container().space();
68  std::vector<Point> vectPoints;
69 
70  // Get points.
71  if ( verbose ) trace.beginBlock( "Getting points." );
72  std::set<Point> pointSet;
73  for ( ConstIterator it = mySurface->begin(), itE = mySurface->end(); it != itE; ++it )
74  getPoints( std::inserter( pointSet, pointSet.begin() ), *it );
75  vectPoints.resize( pointSet.size() );
76  std::copy( pointSet.begin(), pointSet.end(), vectPoints.begin() );
77  pointSet.clear();
78  if ( verbose ) trace.endBlock();
79 
80  // Compute Voronoi Covariance Matrix for all points.
81  myVCM.init( vectPoints.begin(), vectPoints.end() );
82 
83  // Compute VCM( chi_r ) for each point.
84  if ( verbose ) trace.beginBlock ( "Integrating VCM( chi_r(p) ) for each point." );
85  int i = 0;
86  // HatPointFunction< Point, Scalar > chi_r( 1.0, r );
87  for ( typename std::vector<Point>::const_iterator it = vectPoints.begin(), itE = vectPoints.end();
88  it != itE; ++it )
89  {
90  if ( verbose ) trace.progressBar( ++i, vectPoints.size() );
91  Point p = *it;
92  MatrixNN measure = myVCM.measure( myChi, p );
93  // On diagonalise le résultat.
94  EigenStructure & evcm = myPt2EigenStructure[ p ];
95  LinearAlgebraTool::getEigenDecomposition( measure, evcm.vectors, evcm.values );
96  }
97  myVCM.clean(); // free some memory.
98  if ( verbose ) trace.endBlock();
99 
100  if ( verbose ) trace.beginBlock ( "Computing average orientation for each surfel." );
101  typedef functors::HatFunction<Scalar> Functor;
102  Functor fct( 1.0, myRadiusTrivial );
103  LpMetric<Space> l2(2.0); //L2 metric in R^3 for surface propagation.
104  typedef functors::ElementaryConvolutionNormalVectorEstimator< Surfel, CanonicSCellEmbedder<KSpace> >
105  SurfelFunctor;
106  typedef LocalEstimatorFromSurfelFunctorAdapter< DigitalSurfaceContainer, LpMetric<Space>, SurfelFunctor, Functor>
107  NormalEstimator;
108 
109  CanonicSCellEmbedder<KSpace> canonic_embedder( ks );
110  SurfelFunctor surfelFct( canonic_embedder, 1.0 );
111  NormalEstimator estimator;
112  estimator.attach( *mySurface);
113  estimator.setParams( l2, surfelFct, fct , myRadiusTrivial);
114  estimator.init( 1.0, mySurface->begin(), mySurface->end());
115  i = 0;
116  std::vector<Point> pts;
117  int surf_size = mySurface->size();
118  for ( ConstIterator it = mySurface->begin(), itE = mySurface->end(); it != itE; ++it )
119  {
120  if ( verbose ) trace.progressBar(++i, surf_size );
121  Surfel s = *it;
122  Normals & normals = mySurfel2Normals[ s ];
123  // get rough estimation of normal
124  normals.trivialNormal = estimator.eval( it );
125  // get points associated with surfel s
126  getPoints( std::back_inserter( pts ), s );
127  for ( typename std::vector<Point>::const_iterator itPts = pts.begin(), itPtsE = pts.end();
128  itPts != itPtsE; ++itPts )
129  {
130  Point p = *itPts;
131  const EigenStructure& evcm = myPt2EigenStructure[ p ];
132  VectorN n = evcm.vectors.column( Space::dimension-1 );
133  if ( n.dot( normals.trivialNormal ) < 0 ) normals.vcmNormal -= n;
134  else normals.vcmNormal += n;
135  }
136  if ( pts.size() > 1 ) normals.vcmNormal /= pts.size();
137  pts.clear();
138  }
139  if ( verbose ) trace.endBlock();
140 
141  if ( verbose ) trace.endBlock();
142 }
143 
144 //-----------------------------------------------------------------------------
145 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
146 inline
147 DGtal::CountedConstPtrOrConstPtr< typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Surface >
148 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
149 surface() const
150 {
151  return mySurface;
152 }
153 //-----------------------------------------------------------------------------
154 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
155 inline
156 DGtal::Surfel2PointEmbedding
157 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
158 surfelEmbedding() const
159 {
160  return mySurfelEmbedding;
161 }
162 //-----------------------------------------------------------------------------
163 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
164 inline
165 typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Scalar
166 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
167 R() const
168 {
169  return myVCM.R();
170 }
171 //-----------------------------------------------------------------------------
172 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
173 inline
174 typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Scalar
175 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
176 r() const
177 {
178  return myVCM.r();
179 }
180 //-----------------------------------------------------------------------------
181 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
182 inline
183 typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Scalar
184 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
185 radiusTrivial() const
186 {
187  return myRadiusTrivial;
188 }
189 //-----------------------------------------------------------------------------
190 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
191 inline
192 const typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Surfel2Normals&
193 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
194 mapSurfel2Normals() const
195 {
196  return mySurfel2Normals;
197 }
198 //-----------------------------------------------------------------------------
199 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
200 inline
201 const typename DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::Point2EigenStructure&
202 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
203 mapPoint2ChiVCM() const
204 {
205  return myPt2EigenStructure;
206 }
207 
208 //-----------------------------------------------------------------------------
209 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
210 inline
211 bool
212 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
213 getChiVCMEigenvalues( VectorN& values, Surfel s ) const
214 {
215  std::vector<Point> pts;
216  getPoints( std::back_inserter( pts ), s );
217  bool ok = true;
218  int i = 0;
219  values = VectorN(); // Setting values to 0 before averaging.
220  for ( typename std::vector<Point>::const_iterator itPts = pts.begin(), itPtsE = pts.end();
221  itPts != itPtsE; ++itPts, ++i )
222  {
223  Point p = *itPts;
224  typename Point2EigenStructure::const_iterator itEigen = myPt2EigenStructure.find( p );
225  if ( itEigen == myPt2EigenStructure.end() )
226  {
227  ok = false;
228  break;
229  }
230  const EigenStructure& evcm = itEigen->second;
231  values += evcm.values;
232  }
233  if ( i > 1 ) values /= i;
234  return ok;
235 }
236 
237 //-----------------------------------------------------------------------------
238 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
239 inline
240 bool
241 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
242 getChiVCMEigenStructure( VectorN& values, MatrixNN& vectors, Surfel s ) const
243 {
244  std::vector<Point> pts;
245  getPoints( std::back_inserter( pts ), s );
246  bool ok = true;
247  int i = 0;
248  values = VectorN(); // Setting values to 0 before averaging.
249  vectors = MatrixNN(); // Setting values to 0 before averaging.
250  for ( typename std::vector<Point>::const_iterator itPts = pts.begin(), itPtsE = pts.end();
251  itPts != itPtsE; ++itPts, ++i )
252  {
253  Point p = *itPts;
254  typename Point2EigenStructure::const_iterator itEigen = myPt2EigenStructure.find( p );
255  if ( itEigen == myPt2EigenStructure.end() )
256  {
257  ok = false;
258  break;
259  }
260  const EigenStructure& evcm = itEigen->second;
261  values += evcm.values;
262  vectors += evcm.vectors;
263  }
264  if ( i > 1 ) {
265  values /= i;
266  vectors /= i;
267  }
268  return ok;
269 }
270 
271 //-----------------------------------------------------------------------------
272 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
273 template <typename PointOutputIterator>
274 inline
275 PointOutputIterator
276 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
277 getPoints( PointOutputIterator outIt, Surfel s ) const
278 {
279  BOOST_CONCEPT_ASSERT(( boost::OutputIterator< PointOutputIterator, Point > ));
280  const KSpace & ks = mySurface->container().space();
281  Dimension k = ks.sOrthDir( s );
282  switch ( mySurfelEmbedding ) {
283  case Pointels:
284  {
285  typename KSpace::Cells faces = ks.uFaces( ks.unsigns( s ) );
286  for ( typename KSpace::Cells::const_iterator it = faces.begin(), itE = faces.end();
287  it != itE; ++it )
288  {
289  if ( ks.uDim( *it ) == 0 ) // get only pointels (cell of dim 0)
290  *outIt++ = ks.uCoords( *it );
291  }
292  // Dimension i = (k+1)%3;
293  // Dimension j = (i+1)%3;
294  // SCell l1 = ks.sIncident( s, i, true );
295  // SCell l2 = ks.sIncident( s, i, false );
296  // *outIt++ = ks.sCoords( ks.sIncident( l1, j, true ) );
297  // *outIt++ = ks.sCoords( ks.sIncident( l1, j, false ) );
298  // *outIt++ = ks.sCoords( ks.sIncident( l2, j, true ) );
299  // *outIt++ = ks.sCoords( ks.sIncident( l2, j, false ) );
300  }
301  break;
302  case InnerSpel:
303  *outIt++ = ks.sCoords( ks.sDirectIncident( s, k ) );
304  break;
305  case OuterSpel:
306  *outIt++ = ks.sCoords( ks.sIndirectIncident( s, k ) );
307  break;
308  }
309  return outIt;
310 }
311 
312 ///////////////////////////////////////////////////////////////////////////////
313 // Interface - public :
314 
315 /**
316  * Writes/Displays the object on an output stream.
317  * @param out the output stream where the object is written.
318  */
319 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
320 inline
321 void
322 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
323 selfDisplay ( std::ostream & out ) const
324 {
325  out << "[VoronoiCovarianceMeasureOnDigitalSurface"
326  << " #pts=" << myPt2EigenStructure.size()
327  << " #surf=" << mySurfel2Normals.size()
328  << "]";
329 }
330 
331 /**
332  * Checks the validity/consistency of the object.
333  * @return 'true' if the object is valid, 'false' otherwise.
334  */
335 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
336 inline
337 bool
338 DGtal::VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction>::
339 isValid() const
340 {
341  return true;
342 }
343 
344 
345 
346 ///////////////////////////////////////////////////////////////////////////////
347 // Implementation of inline functions //
348 
349 template <typename TDigitalSurfaceContainer, typename TSeparableMetric, typename TKernelFunction>
350 inline
351 std::ostream&
352 DGtal::operator<< ( std::ostream & out,
353  const VoronoiCovarianceMeasureOnDigitalSurface<TDigitalSurfaceContainer, TSeparableMetric, TKernelFunction> & object )
354 {
355  object.selfDisplay( out );
356  return out;
357 }
358 
359 // //
360 ///////////////////////////////////////////////////////////////////////////////
361 
362 
Generated on Tue Jun 1 2021 17:06:10 for DGtal by doxygen 1.9.1