DGtal  1.5.beta
IntegralInvariantNormalVectorEstimator.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 IntegralInvariantNormalVectorEstimator.ih
19  * @author Jeremy Levallois (\c jeremy.levallois@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), INSA-Lyon, France
21  * LAboratoire de MAthématiques - LAMA (CNRS, UMR 5127), Université de Savoie, France
22  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
23  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
24  *
25  * @date 2014/05/12
26  *
27  * Implementation of inline methods defined in IntegralInvariantNormalVectorEstimator.h
28  *
29  * This file is part of the DGtal library.
30  */
31 
32 
33 //////////////////////////////////////////////////////////////////////////////
34 #include <cstdlib>
35 #include "DGtal/math/BasicMathFunctions.h"
36 //////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // IMPLEMENTATION of inline methods.
40 ///////////////////////////////////////////////////////////////////////////////
41 
42 ///////////////////////////////////////////////////////////////////////////////
43 // ----------------------- Standard services ------------------------------
44 
45 //-----------------------------------------------------------------------------
46 template <typename TKSpace, typename TPointPredicate>
47 inline
48 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
49 ~IntegralInvariantNormalVectorEstimator()
50 {
51  clear();
52 }
53 
54 //-----------------------------------------------------------------------------
55 template <typename TKSpace, typename TPointPredicate>
56 inline
57 void
58 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
59 clear()
60 {
61  for( unsigned int i = 0; i < myKernelsSet.size(); ++i )
62  if ( myKernelsSet[ i ] != 0 ) delete myKernelsSet[ i ];
63  myH = 1.0;
64  myRadius = 0.0;
65 }
66 //-----------------------------------------------------------------------------
67 template <typename TKSpace, typename TPointPredicate>
68 inline
69 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
70 IntegralInvariantNormalVectorEstimator()
71  : myKernelFunctor(NumberTraits<Value>::ONE),
72  myKernels(), myKernelsSet(),
73  myKernel( 0 ), myDigKernel( 0 ),
74  myPointPredicate( 0 ), myShapeDomain( 0 ),
75  myShapePointFunctor( 0 ), myShapeSpelFunctor( 0 ),
76  myConvolver( 0 ),
77  myH( 1.0 ), myRadius( 0.0 )
78 {
79 }
80 
81 //-----------------------------------------------------------------------------
82 template <typename TKSpace, typename TPointPredicate>
83 inline
84 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
85 IntegralInvariantNormalVectorEstimator
86 ( ConstAlias< KSpace > K, ConstAlias< PointPredicate > aPointPredicate )
87  : myKernelFunctor(NumberTraits<Value>::ONE),
88  myKernels(), myKernelsSet(),
89  myKernel( 0 ), myDigKernel( 0 ),
90  myPointPredicate( aPointPredicate ), myShapeDomain( 0 ),
91  myShapePointFunctor( 0 ), myShapeSpelFunctor( 0 ),
92  myConvolver( 0 ),
93  myH( 1.0 ), myRadius( 0.0 )
94 {
95  CountedConstPtrOrConstPtr<KSpace> ptrK( K );
96  myShapeDomain = CountedPtr<Domain>( new Domain( ptrK->lowerBound(), ptrK->upperBound() ) );
97  myShapePointFunctor = CountedPtr<ShapePointFunctor>( new ShapePointFunctor( *myPointPredicate, *myShapeDomain, 1, 0 ) );
98  myShapeSpelFunctor = CountedPtr<ShapeSpelFunctor>( new ShapeSpelFunctor( *myShapePointFunctor, K ) );
99  myConvolver = CountedPtr<Convolver>( new Convolver( *myShapeSpelFunctor, myKernelFunctor, K ) );
100 }
101 
102 //-----------------------------------------------------------------------------
103 template <typename TKSpace, typename TPointPredicate>
104 inline
105 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
106 IntegralInvariantNormalVectorEstimator
107 ( const Self& other )
108  : myKernelFunctor( other.myKernelFunctor ),
109  myKernels( other.myKernels ), myKernelsSet( other.myKernelsSet ),
110  myKernel( other.myKernel ), myDigKernel( other.myDigKernel ),
111  myPointPredicate( other.myPointPredicate ), myShapeDomain( other.myShapeDomain ),
112  myShapePointFunctor( other.myShapePointFunctor ), myShapeSpelFunctor( other.myShapeSpelFunctor ),
113  myConvolver( other.myConvolver ),
114  myH( other.myH ), myRadius( other.myRadius )
115 {}
116 //-----------------------------------------------------------------------------
117 template <typename TKSpace, typename TPointPredicate>
118 inline
119 typename DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::Self&
120 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
121 operator= ( const Self& other )
122 {
123  if ( this != &other )
124  {
125  // myKernelFunctor = other.myKernelFunctor;
126  myKernels = other.myKernels;
127  myKernelsSet = other.myKernelsSet;
128  myKernel = other.myKernel;
129  myDigKernel = other.myDigKernel;
130  myPointPredicate = other.myPointPredicate;
131  myShapeDomain = other.myShapeDomain;
132  myShapePointFunctor = other.myShapePointFunctor;
133  myShapeSpelFunctor = other.myShapeSpelFunctor;
134  myConvolver = other.myConvolver;
135  myH = other.myH;
136  myRadius = other.myRadius;
137  }
138  return *this;
139 }
140 //-----------------------------------------------------------------------------
141 template <typename TKSpace, typename TPointPredicate>
142 inline
143 typename DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::Scalar
144 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
145 h() const
146 {
147  return myH;
148 }
149 //-----------------------------------------------------------------------------
150 template <typename TKSpace, typename TPointPredicate>
151 inline
152 void
153 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
154 attach
155 ( ConstAlias< KSpace > K,
156  ConstAlias<PointPredicate> aPointPredicate )
157 {
158  myPointPredicate = aPointPredicate;
159  CountedConstPtrOrConstPtr<KSpace> ptrK( K );
160  myShapeDomain = CountedPtr<Domain>( new Domain( ptrK->lowerBound(), ptrK->upperBound() ) );
161  myShapePointFunctor = CountedPtr<ShapePointFunctor>( new ShapePointFunctor( *myPointPredicate, *myShapeDomain, 1, 0 ) );
162  myShapeSpelFunctor = CountedPtr<ShapeSpelFunctor>( new ShapeSpelFunctor( *myShapePointFunctor, K ) );
163  myConvolver = CountedPtr<Convolver>( new Convolver( *myShapeSpelFunctor, myKernelFunctor, K ) );
164 }
165 //-----------------------------------------------------------------------------
166 template <typename TKSpace, typename TPointPredicate>
167 inline
168 void
169 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
170 setParams
171 ( const double dRadius )
172 {
173  ASSERT( ( dRadius > 0.0 )
174  && "[DGtal::deprecated::IntegralInvariantNormalVectorEstimator:setParams] Radius parameter dRadius must be positive." );
175  myRadius = dRadius;
176 }
177 
178 //-----------------------------------------------------------------------------
179 template <typename TKSpace, typename TPointPredicate>
180 template <typename SurfelConstIterator>
181 inline
182 void
183 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
184 init
185 ( const double _h, SurfelConstIterator /* itb */, SurfelConstIterator /* ite */ )
186 {
187  ASSERT( ( _h > 0.0 )
188  && "[DGtal::deprecated::IntegralInvariantNormalVectorEstimator:init] Gridstep parameter h must be positive." );
189  ASSERT( ( myRadius > 0.0 )
190  && "[DGtal::deprecated::IntegralInvariantNormalVectorEstimator:init] Radius parameter dRadius must have been initialized with a call to 'setParams'." );
191  ASSERT( ( myConvolver != 0 )
192  && "[DGtal::deprecated::IntegralInvariantNormalVectorEstimator:init] Shape of interest must have been initialized with a call to 'attach'." );
193 
194  typedef typename RealPoint::Component Scalar;
195  // Clear stuff
196  for( unsigned int i = 0; i < myKernelsSet.size(); ++i )
197  if ( myKernelsSet[ i ] != 0 ) delete myKernelsSet[ i ];
198 
199  myH = _h;
200  double eRadius = myRadius * myH; // Euclidean radius of the ball kernel.
201 
202  // meanFunctor.init( h, radius );
203 
204  RealPoint rOrigin = RealPoint::zero;
205  Point pOrigin = Point::zero;
206  myKernel = CountedPtr<KernelSupport>( new KernelSupport( rOrigin, eRadius ) ); // acquired
207  myDigKernel = CountedPtr<DigitalShapeKernel>( new DigitalShapeKernel() );
208  myDigKernel->attach( *myKernel );
209  myDigKernel->init( myKernel->getLowerBound() + Point::diagonal(-1), myKernel->getUpperBound() + Point::diagonal(1), myH );
210  Domain neighborhood( Point::diagonal(-1), Point::diagonal(1) );
211  unsigned int n = BasicMathFunctions::power( (unsigned int) 3, Space::dimension );
212  myKernels = std::vector< PairIterators > ( n );
213  myKernelsSet = std::vector< DigitalSet* >( n );
214  unsigned int offset = 0;
215  unsigned int middle = n / 2;
216  RealPoint shiftPoint;
217  for ( typename Domain::ConstIterator it_neigh = neighborhood.begin(),
218  it_neigh_end = neighborhood.end();
219  it_neigh != it_neigh_end;
220  ++it_neigh, ++offset )
221  {
222  /// Computation of shifting masks
223  if( offset == middle ) continue; // no shift
224  for ( Dimension k = 0; k < Space::dimension; ++k )
225  shiftPoint[ k ] = (Scalar) (*it_neigh)[ k ];
226  shiftPoint *= (Scalar) myH;
227  KernelSupport* kernelShifted = new KernelSupport( shiftPoint, eRadius );
228  EuclideanMinus* current = new EuclideanMinus( *myKernel, *kernelShifted );
229  DigitalShape digCurrent;
230  digCurrent.attach( *current );
231  digCurrent.init( myKernel->getLowerBound() + Point::diagonal(-1), myKernel->getUpperBound() + Point::diagonal(1), myH );
232 
233  myKernelsSet[ offset ] = new DigitalSet( digCurrent.getDomain() );
234  Shapes< Domain>::digitalShaper ( *(myKernelsSet[ offset ]), digCurrent );
235 
236  myKernels[ offset ].first = myKernelsSet[ offset ]->begin();
237  myKernels[ offset ].second = myKernelsSet[ offset ]->end();
238 
239  delete current;
240  delete kernelShifted;
241  }
242  /// End of computation of masks
243  myConvolver->init( pOrigin, *myDigKernel, myKernels );
244 }
245 
246 //-----------------------------------------------------------------------------
247 template <typename TKSpace, typename TPointPredicate>
248 template <typename SurfelConstIterator>
249 inline
250 typename DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::Quantity
251 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::
252 eval
253 ( SurfelConstIterator it ) const
254 {
255  CovarianceMatrix2NormalDirectionFunctor normalFunctor;
256  return normalFunctor( myConvolver->evalCovarianceMatrix( it ) );
257 }
258 
259 //-----------------------------------------------------------------------------
260 template <typename TKSpace, typename TPointPredicate>
261 template <typename OutputIterator, typename SurfelConstIterator>
262 inline
263 OutputIterator
264 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::eval
265 ( SurfelConstIterator itb,
266  SurfelConstIterator ite,
267  OutputIterator result ) const
268 {
269  CovarianceMatrix2NormalDirectionFunctor normalFunctor;
270  myConvolver->evalCovarianceMatrix( itb, ite, result, normalFunctor );
271  return result;
272 }
273 
274 //-----------------------------------------------------------------------------
275 template <typename TKSpace, typename TPointPredicate>
276 inline
277 void
278 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::selfDisplay
279 ( std::ostream & out ) const
280 {
281  out << "[IntegralInvariantNormalVectorEstimator h=" << myH
282  << " digR=" << myRadius << " eucR=" << (myH*myRadius) << " ]";
283 }
284 
285 //-----------------------------------------------------------------------------
286 template <typename TKSpace, typename TPointPredicate>
287 inline
288 bool
289 DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate>::isValid() const
290 {
291  return ( myH > 0 ) && ( myRadius > 0 ) && ( myConvolver != 0 );
292 }
293 
294 //-----------------------------------------------------------------------------
295 template <typename TKSpace, typename TPointPredicate>
296 inline
297 std::ostream&
298 DGtal::operator<<
299 ( std::ostream & out,
300  const DGtal::deprecated::IntegralInvariantNormalVectorEstimator<TKSpace, TPointPredicate> & object )
301 {
302  object.selfDisplay( out );
303  return out;
304 }