DGtalTools  1.5.beta
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 ------------------------------
42 template <typename TKSpace, typename TLinearAlgebra>
43 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
44 ATVu2v0( 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 
55 template <typename TKSpace, typename TLinearAlgebra>
56 void
57 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
58 init( 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 //-----------------------------------------------------------------------------
78 template <typename TKSpace, typename TLinearAlgebra>
79 template <typename Image, typename Function>
80 void
81 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
82 addInput( 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 //-----------------------------------------------------------------------------
98 template <typename TKSpace, typename TLinearAlgebra>
99 typename DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::Scalar
100 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
101 computeSNR() 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 //-----------------------------------------------------------------------------
117 template <typename TKSpace, typename TLinearAlgebra>
118 void
119 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
120 setMetricAverage( bool average )
121 {
122  metric_average = average;
123 }
124 
125 //-----------------------------------------------------------------------------
126 template <typename TKSpace, typename TLinearAlgebra>
127 void
128 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
129 setAlpha( 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 //-----------------------------------------------------------------------------
141 template <typename TKSpace, typename TLinearAlgebra>
142 void
143 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
144 setAlpha( 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 //-----------------------------------------------------------------------------
156 template <typename TKSpace, typename TLinearAlgebra>
157 void
158 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
159 setLambda( 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 //-----------------------------------------------------------------------------
173 template <typename TKSpace, typename TLinearAlgebra>
174 void
175 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
176 setEpsilon( 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 //-----------------------------------------------------------------------------
193 template <typename TKSpace, typename TLinearAlgebra>
194 void
195 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
196 setUFromInput()
197 {
198  u2 = g2;
199 }
200 
201 //-----------------------------------------------------------------------------
202 template <typename TKSpace, typename TLinearAlgebra>
203 void
204 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
205 setUFromInputAndMask()
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 //-----------------------------------------------------------------------------
220 template <typename TKSpace, typename TLinearAlgebra>
221 bool
222 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
223 solveU()
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 //-----------------------------------------------------------------------------
247 template <typename TKSpace, typename TLinearAlgebra>
248 bool
249 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
250 solveV()
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 //-----------------------------------------------------------------------------
282 template <typename TKSpace, typename TLinearAlgebra>
283 typename DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::Scalar
284 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
285 computeVariation()
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 //-----------------------------------------------------------------------------
314 template <typename TKSpace, typename TLinearAlgebra>
315 typename DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::Scalar
316 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::
317 checkV()
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  */
345 template <typename TKSpace, typename TLinearAlgebra>
346 inline
347 void
348 DGtal::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  */
357 template <typename TKSpace, typename TLinearAlgebra>
358 inline
359 bool
360 DGtal::ATVu2v0<TKSpace, TLinearAlgebra>::isValid() const
361 {
362  return true;
363 }
364 
365 
366 
367 ///////////////////////////////////////////////////////////////////////////////
368 // Implementation of inline functions //
369 
370 template <typename TKSpace, typename TLinearAlgebra>
371 inline
372 std::ostream&
373 DGtal::operator<< ( std::ostream & out,
374  const ATVu2v0<TKSpace, TLinearAlgebra> & object )
375 {
376  object.selfDisplay( out );
377  return out;
378 }
379 
380 // //
381 ///////////////////////////////////////////////////////////////////////////////
382 
383 
Companion project associated to DGtal ; Generated on Mon Aug 26 2024 21:31:30 for DGtalTools by   doxygen 1.9.1