DGtal 2.1.1
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
46template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
47inline
48DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
49~IntegralInvariantVolumeEstimator()
50{
51 clear();
52}
53
54//-----------------------------------------------------------------------------
55template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
56inline
57void
58DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
59clear()
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//-----------------------------------------------------------------------------
67template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
68inline
69DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
70IntegralInvariantVolumeEstimator( 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//-----------------------------------------------------------------------------
83template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
84inline
85DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
86IntegralInvariantVolumeEstimator
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//-----------------------------------------------------------------------------
107template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
108inline
109DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
110IntegralInvariantVolumeEstimator
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//-----------------------------------------------------------------------------
122template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
123inline
124typename DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::Self&
125DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
126operator= ( 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//-----------------------------------------------------------------------------
147template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
148inline
149typename DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::Scalar
150DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
151h() const
152{
153 return myH;
154}
155//-----------------------------------------------------------------------------
156template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
157inline
158void
159DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
160attach
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//-----------------------------------------------------------------------------
172template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
173inline
174void
175DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
176setParams
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//-----------------------------------------------------------------------------
185template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
186template <typename SurfelConstIterator>
187inline
188void
189DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
190init
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 // Clear stuff
201 for( unsigned int i = 0; i < myKernelsSet.size(); ++i )
202 if ( myKernelsSet[ i ] != 0 ) delete myKernelsSet[ i ];
203
204 myH = _h;
205 double eRadius = myRadius * myH; // Euclidean radius of the ball kernel.
206
207 myFct.init( myH, eRadius );
208
209 RealPoint rOrigin = RealPoint::zero;
210 Point pOrigin = Point::zero;
211 myKernel = CountedPtr<KernelSupport>( new KernelSupport( rOrigin, eRadius ) ); // acquired
212 myDigKernel = CountedPtr<DigitalShapeKernel>( new DigitalShapeKernel() );
213 myDigKernel->attach( *myKernel );
214 myDigKernel->init( myKernel->getLowerBound() + Point::diagonal(-1), myKernel->getUpperBound() + Point::diagonal(1), myH );
215 Domain neighborhood( Point::diagonal(-1), Point::diagonal(1) );
216 unsigned int n = functions::power( (unsigned int) 3, Space::dimension );
217 myKernels = std::vector< PairIterators > ( n );
218 myKernelsSet = std::vector< DigitalSet* >( n );
219 unsigned int offset = 0;
220 unsigned int middle = n / 2;
221 const Domain kernelDomain = myDigKernel->getDomain();
222 for ( typename Domain::ConstIterator it_neigh = neighborhood.begin(),
223 it_neigh_end = neighborhood.end();
224 it_neigh != it_neigh_end;
225 ++it_neigh, ++offset )
226 {
227 /// Computation of shifting masks
228 if( offset == middle ) continue; // no shift
229
230 myKernelsSet[ offset ] = new DigitalSet( kernelDomain );
231 for ( typename Domain::ConstIterator it = kernelDomain.begin(),
232 itEnd = kernelDomain.end();
233 it != itEnd; ++it )
234 {
235 if ( myDigKernel->operator()( *it ) )
236 {
237 const Point shiftedPoint = *it - *it_neigh;
238 const bool isInsideShifted = kernelDomain.isInside( shiftedPoint )
239 && myDigKernel->operator()( shiftedPoint );
240 if ( ! isInsideShifted )
241 myKernelsSet[ offset ]->insert( *it );
242 }
243 }
244
245 myKernels[ offset ].first = myKernelsSet[ offset ]->begin();
246 myKernels[ offset ].second = myKernelsSet[ offset ]->end();
247 }
248 /// End of computation of masks
249 myConvolver->init( pOrigin, *myDigKernel, myKernels );
250}
251
252//-----------------------------------------------------------------------------
253template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
254template <typename SurfelConstIterator>
255inline
256typename DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::Quantity
257DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::
258eval
259( SurfelConstIterator it ) const
260{
261 return myFct( myConvolver->eval( it ) );
262}
263
264//-----------------------------------------------------------------------------
265template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
266template <typename OutputIterator, typename SurfelConstIterator>
267inline
268OutputIterator
269DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::eval
270( SurfelConstIterator itb,
271 SurfelConstIterator ite,
272 OutputIterator result ) const
273{
274 myConvolver->eval( itb, ite, result, myFct );
275 return result;
276}
277
278//-----------------------------------------------------------------------------
279template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
280inline
281void
282DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::selfDisplay
283( std::ostream & out ) const
284{
285 out << "[IntegralInvariantVolumeEstimator h=" << myH
286 << " digR=" << myRadius << " eucR=" << (myH*myRadius) << " ]";
287}
288
289//-----------------------------------------------------------------------------
290template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
291inline
292bool
293DGtal::IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor>::isValid() const
294{
295 return ( myH > 0 ) && ( myRadius > 0 ) && ( myConvolver != 0 );
296}
297
298//-----------------------------------------------------------------------------
299template <typename TKSpace, typename TPointPredicate, typename TVolumeFunctor>
300inline
301std::ostream&
302DGtal::operator<<
303( std::ostream & out,
304 const IntegralInvariantVolumeEstimator<TKSpace, TPointPredicate, TVolumeFunctor> & object )
305{
306 object.selfDisplay( out );
307 return out;
308}