DGtal 1.3.0
Loading...
Searching...
No Matches
IntegralInvariantCovarianceEstimator.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 IntegralInvariantCovarianceEstimator.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 IntegralInvariantCovarianceEstimator.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 TCovarianceMatrixFunctor>
47inline
48DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
49~IntegralInvariantCovarianceEstimator()
50{
51 clear();
52}
53
54//-----------------------------------------------------------------------------
55template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
56inline
57void
58DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
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 TCovarianceMatrixFunctor>
68inline
69DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
70IntegralInvariantCovarianceEstimator( CovarianceMatrixFunctor 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 TCovarianceMatrixFunctor>
84inline
85DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
86IntegralInvariantCovarianceEstimator
87( ConstAlias< KSpace > K,
88 ConstAlias< PointPredicate > aPointPredicate,
89 CovarianceMatrixFunctor 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 TCovarianceMatrixFunctor>
108inline
109DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
110IntegralInvariantCovarianceEstimator
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 TCovarianceMatrixFunctor>
123inline
124typename DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::Self&
125DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
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 TCovarianceMatrixFunctor>
148inline
149typename DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::Scalar
150DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
151h() const
152{
153 return myH;
154}
155//-----------------------------------------------------------------------------
156template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
157inline
158void
159DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
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 TCovarianceMatrixFunctor>
173inline
174void
175DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
176setParams
177( const double dRadius )
178{
179 ASSERT( ( dRadius > 0.0 )
180 && "[DGtal::IntegralInvariantCovarianceEstimator:setParams] Radius parameter dRadius must be positive." );
181 myRadius = dRadius;
182}
183
184//-----------------------------------------------------------------------------
185template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
186template <typename SurfelConstIterator>
187inline
188void
189DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
190init
191( const double _h, SurfelConstIterator /* itb */, SurfelConstIterator /* ite */ )
192{
193 ASSERT( ( _h > 0.0 )
194 && "[DGtal::IntegralInvariantCovarianceEstimator:init] Gridstep parameter h must be positive." );
195 ASSERT( ( myRadius > 0.0 )
196 && "[DGtal::IntegralInvariantCovarianceEstimator:init] Radius parameter dRadius must have been initialized with a call to 'setParams'." );
197 ASSERT( ( myConvolver != 0 )
198 && "[DGtal::IntegralInvariantCovarianceEstimator:init] Shape of interest must have been initialized with a call to 'attach'." );
199
200 typedef typename RealPoint::Component ScalarC;
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 ] = (ScalarC) (*it_neigh)[ k ];
232 shiftPoint *= (ScalarC) 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//-----------------------------------------------------------------------------
254template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
255template <typename SurfelConstIterator>
256inline
257typename DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::Quantity
258DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::
259eval
260( SurfelConstIterator it ) const
261{
262 return myFct( myConvolver->evalCovarianceMatrix( it ) );
263}
264
265//-----------------------------------------------------------------------------
266template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
267template <typename OutputIterator, typename SurfelConstIterator>
268inline
269OutputIterator
270DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::eval
271( SurfelConstIterator itb,
272 SurfelConstIterator ite,
273 OutputIterator result ) const
274{
275 myConvolver->evalCovarianceMatrix( itb, ite, result, myFct );
276 return result;
277}
278
279//-----------------------------------------------------------------------------
280template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
281inline
282void
283DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::selfDisplay
284( std::ostream & out ) const
285{
286 out << "[IntegralInvariantCovarianceEstimator h=" << myH
287 << " digR=" << myRadius << " eucR=" << (myH*myRadius) << " ]";
288}
289
290//-----------------------------------------------------------------------------
291template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
292inline
293bool
294DGtal::IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor>::isValid() const
295{
296 return ( myH > 0 ) && ( myRadius > 0 ) && ( myConvolver != 0 );
297}
298
299//-----------------------------------------------------------------------------
300template <typename TKSpace, typename TPointPredicate, typename TCovarianceMatrixFunctor>
301inline
302std::ostream&
303DGtal::operator<<
304( std::ostream & out,
305 const IntegralInvariantCovarianceEstimator<TKSpace, TPointPredicate, TCovarianceMatrixFunctor> & object )
306{
307 object.selfDisplay( out );
308 return out;
309}