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