DGtal  0.9.2
IntegralInvariantVolumeEstimator.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 IntegralInvariantVolumeEstimator.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 IntegralInvariantVolumeEstimator.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, typename TVolumeFunctor>
47 inline
48 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
49 ~IntegralInvariantVolumeEstimator()
50 {
51  clear();
52 }
53 
54 //-----------------------------------------------------------------------------
55 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
56 inline
57 void
58 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
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, typename TVolumeFunctor>
68 inline
69 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
70 IntegralInvariantVolumeEstimator( VolumeFunctor fct )
71  : myFct( fct ),
72  myKernelFunctor(NumberTraits<Value>::ONE),
73  myKernels(), myKernelsSet(),
74  myKernel( 0 ), myDigKernel( 0 ),
75  myPointPredicate( 0 ), myShapeDomain( 0 ),
76  myShapePointFunctor( 0 ), myShapeSpelFunctor( 0 ),
77  myConvolver( 0 ),
78  myH( 1.0 ), myRadius( 0.0 )
79 {
80 }
81 
82 //-----------------------------------------------------------------------------
83 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
84 inline
85 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
86 IntegralInvariantVolumeEstimator
87 ( ConstAlias< KSpace > K,
88  ConstAlias< PointPredicate > aPointPredicate,
89  VolumeFunctor fct )
90  : myFct( fct ),
91  myKernelFunctor(NumberTraits<Value>::ONE),
92  myKernels(), myKernelsSet(),
93  myKernel( 0 ), myDigKernel( 0 ),
94  myPointPredicate( aPointPredicate ), myShapeDomain( 0 ),
95  myShapePointFunctor( 0 ), myShapeSpelFunctor( 0 ),
96  myConvolver( 0 ),
97  myH( 1.0 ), myRadius( 0.0 )
98 {
99  CountedConstPtrOrConstPtr<KSpace> ptrK( K );
100  myShapeDomain = CountedPtr<Domain>( new Domain( ptrK->lowerBound(), ptrK->upperBound() ) );
101  myShapePointFunctor = CountedPtr<ShapePointFunctor>( new ShapePointFunctor( *myPointPredicate, *myShapeDomain, 1, 0 ) );
102  myShapeSpelFunctor = CountedPtr<ShapeSpelFunctor>( new ShapeSpelFunctor( *myShapePointFunctor, K ) );
103  myConvolver = CountedPtr<Convolver>( new Convolver( *myShapeSpelFunctor, myKernelFunctor, K ) );
104 }
105 
106 //-----------------------------------------------------------------------------
107 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
108 inline
109 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
110 IntegralInvariantVolumeEstimator
111 ( const Self& other )
112  : myFct( other.myFct ),
113  myKernelFunctor( other.myKernelFunctor ),
114  myKernels( other.myKernels ), myKernelsSet( other.myKernelsSet ),
115  myKernel( other.myKernel ), myDigKernel( other.myDigKernel ),
116  myPointPredicate( other.myPointPredicate ), myShapeDomain( other.myShapeDomain ),
117  myShapePointFunctor( other.myShapePointFunctor ), myShapeSpelFunctor( other.myShapeSpelFunctor ),
118  myConvolver( other.myConvolver ),
119  myH( other.myH ), myRadius( other.myRadius )
120 {}
121 //-----------------------------------------------------------------------------
122 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
123 inline
124 typename DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::Self&
125 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
126 operator= ( const Self& other )
127 {
128  if ( this != &other )
129  {
130  myFct = other.myFct;
131  // myKernelFunctor = other.myKernelFunctor;
132  myKernels = other.myKernels;
133  myKernelsSet = other.myKernelsSet;
134  myKernel = other.myKernel;
135  myDigKernel = other.myDigKernel;
136  myPointPredicate = other.myPointPredicate;
137  myShapeDomain = other.myShapeDomain;
138  myShapePointFunctor = other.myShapePointFunctor;
139  myShapeSpelFunctor = other.myShapeSpelFunctor;
140  myConvolver = other.myConvolver;
141  myH = other.myH;
142  myRadius = other.myRadius;
143  }
144  return *this;
145 }
146 //-----------------------------------------------------------------------------
147 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
148 inline
149 typename DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::Scalar
150 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
151 h() const
152 {
153  return myH;
154 }
155 //-----------------------------------------------------------------------------
156 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
157 inline
158 void
159 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
160 attach
161 ( ConstAlias< KSpace > K,
162  ConstAlias<PointPredicate> aPointPredicate )
163 {
164  myPointPredicate = aPointPredicate;
165  CountedConstPtrOrConstPtr<KSpace> ptrK( K );
166  myShapeDomain = CountedPtr<Domain>( new Domain( ptrK->lowerBound(), ptrK->upperBound() ) );
167  myShapePointFunctor = CountedPtr<ShapePointFunctor>( new ShapePointFunctor( *myPointPredicate, *myShapeDomain, 1, 0 ) );
168  myShapeSpelFunctor = CountedPtr<ShapeSpelFunctor>( new ShapeSpelFunctor( *myShapePointFunctor, K ) );
169  myConvolver = CountedPtr<Convolver>( new Convolver( *myShapeSpelFunctor, myKernelFunctor, K ) );
170 }
171 //-----------------------------------------------------------------------------
172 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
173 inline
174 void
175 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
176 setParams
177 ( const double dRadius )
178 {
179  ASSERT( ( dRadius > 0.0 )
180  && "[DGtal::IntegralInvariantVolumeEstimator:setParams] Radius parameter dRadius must be positive." );
181  myRadius = dRadius;
182 }
183 
184 //-----------------------------------------------------------------------------
185 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
186 template <typename SurfelConstIterator>
187 inline
188 void
189 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
190 init
191 ( const double _h, SurfelConstIterator /* itb */, SurfelConstIterator /* ite */ )
192 {
193  ASSERT( ( _h > 0.0 )
194  && "[DGtal::IntegralInvariantVolumeEstimator:init] Gridstep parameter h must be positive." );
195  ASSERT( ( myRadius > 0.0 )
196  && "[DGtal::IntegralInvariantVolumeEstimator:init] Radius parameter dRadius must have been initialized with a call to 'setParams'." );
197  ASSERT( ( myConvolver != 0 )
198  && "[DGtal::IntegralInvariantVolumeEstimator:init] Shape of interest must have been initialized with a call to 'attach'." );
199 
200  typedef typename RealPoint::Component Scalar;
201  // Clear stuff
202  for( unsigned int i = 0; i < myKernelsSet.size(); ++i )
203  if ( myKernelsSet[ i ] != 0 ) delete myKernelsSet[ i ];
204 
205  myH = _h;
206  double eRadius = myRadius * myH; // Euclidean radius of the ball kernel.
207 
208  myFct.init( myH, eRadius );
209 
210  RealPoint rOrigin = RealPoint::zero;
211  Point pOrigin = Point::zero;
212  myKernel = CountedPtr<KernelSupport>( new KernelSupport( rOrigin, eRadius ) ); // acquired
213  myDigKernel = CountedPtr<DigitalShapeKernel>( new DigitalShapeKernel() );
214  myDigKernel->attach( *myKernel );
215  myDigKernel->init( myKernel->getLowerBound() + Point::diagonal(-1), myKernel->getUpperBound() + Point::diagonal(1), myH );
216  Domain neighborhood( Point::diagonal(-1), Point::diagonal(1) );
217  unsigned int n = functions::power( (unsigned int) 3, Space::dimension );
218  myKernels = std::vector< PairIterators > ( n );
219  myKernelsSet = std::vector< DigitalSet* >( n );
220  unsigned int offset = 0;
221  unsigned int middle = n / 2;
222  RealPoint shiftPoint;
223  for ( typename Domain::ConstIterator it_neigh = neighborhood.begin(),
224  it_neigh_end = neighborhood.end();
225  it_neigh != it_neigh_end;
226  ++it_neigh, ++offset )
227  {
228  /// Computation of shifting masks
229  if( offset == middle ) continue; // no shift
230  for ( Dimension k = 0; k < Space::dimension; ++k )
231  shiftPoint[ k ] = (Scalar) (*it_neigh)[ k ];
232  shiftPoint *= (Scalar) myH;
233  KernelSupport* kernelShifted = new KernelSupport( shiftPoint, eRadius );
234  EuclideanMinus* current = new EuclideanMinus( myKernel );
235  current->minus( kernelShifted );
236  DigitalShape digCurrent;
237  digCurrent.attach( *current );
238  digCurrent.init( myKernel->getLowerBound() + Point::diagonal(-1), myKernel->getUpperBound() + Point::diagonal(1), myH );
239 
240  myKernelsSet[ offset ] = new DigitalSet( digCurrent.getDomain() );
241  Shapes< Domain>::digitalShaper ( *(myKernelsSet[ offset ]), digCurrent );
242 
243  myKernels[ offset ].first = myKernelsSet[ offset ]->begin();
244  myKernels[ offset ].second = myKernelsSet[ offset ]->end();
245 
246  delete current;
247  delete kernelShifted;
248  }
249  /// End of computation of masks
250  myConvolver->init( pOrigin, *myDigKernel, myKernels );
251 }
252 
253 //-----------------------------------------------------------------------------
254 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
255 template <typename SurfelConstIterator>
256 inline
257 typename DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::Quantity
258 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
259 eval
260 ( SurfelConstIterator it ) const
261 {
262  return myFct( myConvolver->eval( it ) );
263 }
264 
265 //-----------------------------------------------------------------------------
266 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
267 template <typename OutputIterator, typename SurfelConstIterator>
268 inline
269 OutputIterator
270 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::eval
271 ( SurfelConstIterator itb,
272  SurfelConstIterator ite,
273  OutputIterator result ) const
274 {
275  myConvolver->eval( itb, ite, result, myFct );
276  return result;
277 }
278 
279 //-----------------------------------------------------------------------------
280 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
281 inline
282 void
283 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::selfDisplay
284 ( std::ostream & out ) const
285 {
286  out << "[IntegralInvariantVolumeEstimator h=" << myH
287  << " digR=" << myRadius << " eucR=" << (myH*myRadius) << " ]";
288 }
289 
290 //-----------------------------------------------------------------------------
291 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
292 inline
293 bool
294 DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::isValid() const
295 {
296  return ( myH > 0 ) && ( myRadius > 0 ) && ( myConvolver != 0 );
297 }
298 
299 //-----------------------------------------------------------------------------
300 template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
301 inline
302 std::ostream&
303 DGtal::operator<<
304 ( std::ostream & out,
305  const IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor> & object )
306 {
307  object.selfDisplay( out );
308  return out;
309 }