DGtalTools 2.0.0
Loading...
Searching...
No Matches
ATVu2v0.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 ATVu2v0.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 * @author Marion Foare (\c marion.foare@univ-savoie.fr )
22 * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
23 *
24 * @date 2016/10/12
25 *
26 * Implementation of inline methods defined in ATVu2v0.h
27 *
28 * This file is part of the DGtal library.
29 */
30
31
32//////////////////////////////////////////////////////////////////////////////
33#include <cstdlib>
34//////////////////////////////////////////////////////////////////////////////
35
36///////////////////////////////////////////////////////////////////////////////
37// IMPLEMENTATION of inline methods.
38///////////////////////////////////////////////////////////////////////////////
39
40///////////////////////////////////////////////////////////////////////////////
41// ----------------------- Standard services ------------------------------
42template <typename TKSpace, typename TLinearAlgebra>
43DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
44ATVu2v0( int _verbose )
45 : Base( _verbose ),
46 M01( calculus ), M12( calculus ),
47 primal_AD2( calculus ), primal_L0( calculus ),
48 v0( calculus ), former_v0( calculus ),
49 alpha_Id2( calculus ), left_V0( calculus ),
50 l_1( calculus ), l_over_e_1( calculus ),
51 l_over_e_Id0( calculus ),
52 metric_average( false )
53{}
54
55template <typename TKSpace, typename TLinearAlgebra>
56void
57DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
58init( Clone<KSpace> aKSpace )
59{
60 Base::init( aKSpace );
61 if ( verbose > 0 ) trace.beginBlock( "Initialize DEC specific operators" );
62 if ( verbose > 1 ) trace.info() << "M01" << std::endl;
63 M01 = calculus.template derivative<0, PRIMAL>();
64 M01.myContainer = .5 * M01.myContainer.cwiseAbs();
65 if ( verbose > 1 ) trace.info() << "M12" << std::endl;
66 M12 = calculus.template derivative<1, PRIMAL>();
67 M12.myContainer = .25 * M12.myContainer.cwiseAbs();
68 if ( verbose > 1 ) trace.info() << "primal_AD2" << std::endl;
69 primal_AD2 = calculus.template antiderivative<2,PRIMAL>();
70 if ( verbose > 1 ) trace.info() << "primal_L0" << std::endl;
71 primal_L0 = calculus.template laplace<PRIMAL>();
72 if ( verbose > 1 ) trace.info() << "v0" << std::endl;
73 v0 = KForm<Calculus, 0, PRIMAL>::ones( calculus );
74 if ( verbose > 0 ) trace.endBlock();
75}
76
77//-----------------------------------------------------------------------------
78template <typename TKSpace, typename TLinearAlgebra>
79template <typename Image, typename Function>
80void
81DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
82addInput( const Image& image,
83 const Function& f,
84 bool perfect_data )
85{
86 if ( perfect_data ) i2.push_back( PrimalForm2( calculus ) );
87 else g2.push_back( PrimalForm2( calculus ) );
88 PrimalForm2& g = perfect_data ? i2.back() : g2.back();
89 const KSpace& K = calculus.myKSpace;
90 for ( Index index = 0; index < g.myContainer.rows(); index++)
91 {
92 SCell cell = g.getSCell( index );
93 g.myContainer( index ) = f( image( K.sCoords( cell ) ) );
94 }
95}
96
97//-----------------------------------------------------------------------------
98template <typename TKSpace, typename TLinearAlgebra>
99typename DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::Scalar
100DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
101computeSNR() const
102{
103 Scalar MSE = 0.0;
104 for ( Dimension i = 0; i < i2.size(); ++i )
105 {
106 Scalar MSEi = 0.0;
107 const PrimalForm2 u_minus_i_snr = u2[ i ] - i2[ i ];
108 for ( Index j = 0; j < u_minus_i_snr.length(); ++j )
109 MSEi += u_minus_i_snr.myContainer( j ) * u_minus_i_snr.myContainer( j );
110 MSE += MSEi / (Scalar) u_minus_i_snr.length();
111 }
112 MSE /= 3.0;
113 return 10.0 * log10(1.0 / MSE);
114}
115
116//-----------------------------------------------------------------------------
117template <typename TKSpace, typename TLinearAlgebra>
118void
119DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
120setMetricAverage( bool average )
121{
122 metric_average = average;
123}
124
125//-----------------------------------------------------------------------------
126template <typename TKSpace, typename TLinearAlgebra>
127void
128DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
129setAlpha( Scalar _alpha )
130{
131 ASSERT( _alpha >= 0.0 );
132 alpha = _alpha;
133 // Building alpha_Id2
134 alpha_Id2 = _alpha * calculus.template identity<2, PRIMAL>();
135 alpha_g2.clear();
136 for ( unsigned int i = 0; i < g2.size(); i++ )
137 alpha_g2.push_back( alpha_Id2 * g2[ i ] );
138}
139
140//-----------------------------------------------------------------------------
141template <typename TKSpace, typename TLinearAlgebra>
142void
143DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
144setAlpha( Scalar _alpha, const PrimalForm2& m )
145{
146 ASSERT( _alpha >= 0.0 );
147 alpha = _alpha;
148 // Building alpha_Id0
149 alpha_Id2 = _alpha * functions::dec::diagonal( m ) * calculus.template identity<2, PRIMAL>();
150 alpha_g2.clear();
151 for ( unsigned int i = 0; i < g2.size(); i++ )
152 alpha_g2.push_back( alpha_Id2 * g2[ i ] );
153}
154
155//-----------------------------------------------------------------------------
156template <typename TKSpace, typename TLinearAlgebra>
157void
158DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
159setLambda( Scalar _lambda )
160{
161 ASSERT( _lambda >= 0.0 );
162 lambda = _lambda;
163 if ( metric_average )
164 l_1 = lambda
165 * M01.transpose() * M12.transpose()
166 * M12 * M01
167 * KForm< Calculus, 0, PRIMAL >::ones( calculus );
168 else
169 l_1 = lambda * KForm< Calculus, 0, PRIMAL >::ones( calculus );
170}
171
172//-----------------------------------------------------------------------------
173template <typename TKSpace, typename TLinearAlgebra>
174void
175DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
176setEpsilon( Scalar _epsilon )
177{
178 ASSERT( _epsilon > 0.0 );
179 epsilon = _epsilon;
180 if ( metric_average )
181 left_V0 = ( lambda/epsilon )
182 * M01.transpose() * M12.transpose()
183 * M12 * M01
184 + (lambda*epsilon*epsilon*epsilon) * primal_L0.transpose() * primal_L0;
185 else
186 left_V0 = ( lambda/epsilon )
187 * calculus.template identity<0, PRIMAL>()
188 + (lambda*epsilon*epsilon*epsilon) * primal_L0.transpose() * primal_L0;
189 l_over_e_1 = (1.0/epsilon) * l_1;
190}
191
192//-----------------------------------------------------------------------------
193template <typename TKSpace, typename TLinearAlgebra>
194void
195DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
196setUFromInput()
197{
198 u2 = g2;
199}
200
201//-----------------------------------------------------------------------------
202template <typename TKSpace, typename TLinearAlgebra>
203void
204DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
205setUFromInputAndMask()
206{
207 setUFromInput();
208 for ( int i = 0; i < u2.size(); ++i )
209 {
210 const Index nb = u2[ i ].length();
211 for ( Index index = 0; index < nb; index++)
212 {
213 if ( alpha_g2[ i ].myContainer( index ) == 0.0 )
214 u2[ i ].myContainer( index ) = ((double) rand()) / (double) RAND_MAX;
215 }
216 }
217}
218
219//-----------------------------------------------------------------------------
220template <typename TKSpace, typename TLinearAlgebra>
221bool
222DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
223solveU()
224{
225 if ( verbose > 0 ) trace.beginBlock("Solving for u");
226 if ( verbose > 1 ) trace.info() << "- building matrix M : = alpha_Id2 - tB'_Diag(M01 v)^2_B'" << std::endl;
227
228 const PrimalIdentity1 diag_Mv_squared = functions::dec::squaredDiagonal( M01 * v0 );
229 // JOL: clarify sign below.
230 const PrimalIdentity2 M = alpha_Id2
231 + primal_AD2.transpose() * diag_Mv_squared * primal_AD2;
232 if ( verbose > 1 ) trace.info() << "- prefactoring matrix M" << std::endl;
233 solver_u.compute( M );
234 bool ok = true;
235 for ( Dimension i = 0; i < u2.size(); ++i )
236 {
237 if ( verbose > 1 ) trace.info() << "- solving M u[" << i << "] = alpha g[" << i << "]" << std::endl;
238 u2[ i ] = solver_u.solve( alpha_g2[ i ] );
239 if ( verbose > 1 ) trace.info() << ( solver_u.isValid() ? "=> OK" : "ERROR" ) << " " << solver_u.myLinearAlgebraSolver.info() << std::endl;
240 ok = ok && solver_u.isValid();
241 }
242 if ( verbose > 0 ) trace.endBlock();
243 return ok;
244}
245
246//-----------------------------------------------------------------------------
247template <typename TKSpace, typename TLinearAlgebra>
248bool
249DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
250solveV()
251{
252 former_v0 = v0;
253 if ( verbose > 0 ) trace.beginBlock("Solving for v");
254 if ( verbose > 1 ) trace.info() << "- building matrix N := l/e Id0 + le^3 L^t L + M01^t sum Diag(B' u_i)^2 M01" << std::endl;
255
256 PrimalForm1 squared_norm_D_u2 = PrimalForm1::zeros(calculus);
257
258 PrimalIdentity1 U2 = functions::dec::squaredDiagonal( primal_AD2 * u2[ 0 ] );
259 for ( Dimension i = 1; i < u2.size(); ++i )
260 U2.myContainer += functions::dec::squaredDiagonal( primal_AD2 * u2[ i ] ).myContainer;
261 const PrimalIdentity0 N = left_V0 + M01.transpose() * U2 * M01;
262
263 typedef typename PrimalIdentity0::Container Matrix;
264 const Matrix & M = N.myContainer;
265 if ( verbose > 2 )
266 for (int k = 0; k < M.outerSize(); ++k)
267 for ( typename Matrix::InnerIterator it( M, k ); it; ++it )
268 if ( ( verbose > 3 ) || ( it.row() == it.col() ) )
269 trace.info() << "[" << it.row() << "," << it.col() << "] = " << it.value() << std::endl;
270 if ( verbose > 1 ) trace.info() << "- prefactoring matrix N" << std::endl;
271 solver_v.compute( N );
272 if ( verbose > 1 ) trace.info() << "- solving N v = l/e 1" << std::endl;
273 v0 = solver_v.solve( l_over_e_1 );
274 if ( verbose > 1 ) trace.info() << ( solver_v.isValid() ? "OK" : "ERROR" )
275 << " " << solver_v.myLinearAlgebraSolver.info()
276 << std::endl;
277 if ( verbose > 0 ) trace.endBlock();
278 return solver_v.isValid();
279}
280
281//-----------------------------------------------------------------------------
282template <typename TKSpace, typename TLinearAlgebra>
283typename DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::Scalar
284DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
285computeVariation()
286{
287 if ( verbose > 0 ) trace.beginBlock( "Compute variation of v.");
288 delta_v_l1 = 0.0;
289 delta_v_l2 = 0.0;
290 delta_v_loo = 0.0;
291 for ( Index index = 0; index < size0(); index++)
292 {
293 delta_v_loo = std::max( delta_v_loo,
294 std::fabs( v0.myContainer( index )
295 - former_v0.myContainer( index ) ) );
296 delta_v_l2 += ( v0.myContainer( index ) - former_v0.myContainer( index ) )
297 * ( v0.myContainer( index ) - former_v0.myContainer( index ) );
298 delta_v_l1 += fabs( v0.myContainer( index )
299 - former_v0.myContainer( index ) );
300 }
301 delta_v_l1 /= size0();
302 delta_v_l2 = sqrt( delta_v_l2 / size0() );
303 if ( verbose > 0 )
304 {
305 trace.info() << "Variation |v^k+1 - v^k|_oo = " << delta_v_loo << std::endl;
306 trace.info() << "Variation |v^k+1 - v^k|_2 = " << delta_v_l2 << std::endl;
307 trace.info() << "Variation |v^k+1 - v^k|_1 = " << delta_v_l1 << std::endl;
308 }
309 if ( verbose > 0 ) trace.endBlock();
310 return delta_v_loo;
311}
312
313//-----------------------------------------------------------------------------
314template <typename TKSpace, typename TLinearAlgebra>
315typename DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::Scalar
316DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
317checkV()
318{
319 if ( verbose > 0 ) trace.beginBlock("Checking v");
320 Scalar m1 = 1.0;
321 Scalar m2 = 0.0;
322 Scalar ma = 0.0;
323 for ( Index index = 0; index < size0(); index++)
324 {
325 Scalar val = v0.myContainer( index );
326 m1 = std::min( m1, val );
327 m2 = std::max( m2, val );
328 ma += val;
329 }
330 if ( verbose > 0 )
331 trace.info() << "1-form v: min=" << m1 << " avg=" << ( ma / size0() )
332 << " max=" << m2 << std::endl;
333 if ( verbose > 0 ) trace.endBlock();
334 return std::max( std::fabs( m1 ), std::fabs( m2 - 1.0 ) );
335}
336
337
338///////////////////////////////////////////////////////////////////////////////
339// Interface - public :
340
341/**
342 * Writes/Displays the object on an output stream.
343 * @param out the output stream where the object is written.
344 */
345template <typename TKSpace, typename TLinearAlgebra>
346inline
347void
348DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::selfDisplay ( std::ostream & out ) const
349{
350 out << "[ ATVu2v0 #g=" << g2.size() << " dec=" << calculus << " ]";
351}
352
353/**
354 * Checks the validity/consistency of the object.
355 * @return 'true' if the object is valid, 'false' otherwise.
356 */
357template <typename TKSpace, typename TLinearAlgebra>
358inline
359bool
360DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::isValid() const
361{
362 return true;
363}
364
365
366
367///////////////////////////////////////////////////////////////////////////////
368// Implementation of inline functions //
369
370template <typename TKSpace, typename TLinearAlgebra>
371inline
372std::ostream&
373DGtal::operator<< ( std::ostream & out,
374 const ATVu2v0<TKSpace, TLinearAlgebra> & object )
375{
376 object.selfDisplay( out );
377 return out;
378}
379
380// //
381///////////////////////////////////////////////////////////////////////////////
382
383