DGtalTools  1.5.beta
ATu0v1.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
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 ATu0v1.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::ATu0v1<TKSpace, TLinearAlgebra>::
44 ATu0v1( int _verbose )
45  : Base( _verbose ),
46  v1( calculus ), former_v1( calculus ),
47  L1( calculus ), alpha_Id0( calculus ),
48  l_L1( calculus ), l_1_over_4( calculus ),
49  left_V1( calculus ), l_1_over_4e( calculus )
50 {}
51 
52 template <typename TKSpace, typename TLinearAlgebra>
53 void
54 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
55 init( Clone<KSpace> aKSpace )
56 {
57  Base::init( aKSpace );
58  if ( verbose > 0 ) trace.beginBlock( "Initialize DEC specific operators" );
59  if ( verbose > 1 ) trace.info() << "primal_L1" << std::endl;
60  L1 = -1.0 * ( D0 * dual_h2 * dual_D1 * primal_h1
61  + dual_h1 * dual_D0 * primal_h2 * D1 );
62  // JOL: Does not work ! Sign problem ?
63  // L1 = -1.0 * ( D0 * AD1 + AD2 * D1 );
64  if ( verbose > 1 ) trace.info() << "v1" << std::endl;
65  v1 = KForm<Calculus, 1, PRIMAL>::ones( calculus );
66  if ( verbose > 0 ) trace.endBlock();
67 }
68 
69 //-----------------------------------------------------------------------------
70 template <typename TKSpace, typename TLinearAlgebra>
71 template <typename Image, typename Function>
72 void
73 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
74 addInput( const Image& image,
75  const Function& f,
76  bool perfect_data )
77 {
78  if ( perfect_data ) i0.push_back( PrimalForm0( calculus ) );
79  else g0.push_back( PrimalForm0( calculus ) );
80  PrimalForm0& g = perfect_data ? i0.back() : g0.back();
81  const KSpace& K = calculus.myKSpace;
82  for ( Index index = 0; index < g.myContainer.rows(); index++)
83  {
84  SCell cell = g.getSCell( index );
85  g.myContainer( index ) = f( image( K.sCoords( cell ) ) );
86  }
87 }
88 
89 //-----------------------------------------------------------------------------
90 template <typename TKSpace, typename TLinearAlgebra>
91 typename DGtal::ATu0v1<TKSpace, TLinearAlgebra>::Scalar
92 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
93 computeSNR() const
94 {
95  Scalar MSE = 0.0;
96  for ( Dimension i = 0; i < i0.size(); ++i )
97  {
98  Scalar MSEi = 0.0;
99  const PrimalForm0 u_minus_i_snr = u0[ i ] - i0[ i ];
100  for ( Index j = 0; j < u_minus_i_snr.length(); ++j )
101  MSEi += u_minus_i_snr.myContainer( j ) * u_minus_i_snr.myContainer( j );
102  MSE += MSEi / (Scalar) u_minus_i_snr.length();
103  }
104  MSE /= 3.0;
105  return 10.0 * log10(1.0 / MSE);
106 }
107 
108 //-----------------------------------------------------------------------------
109 template <typename TKSpace, typename TLinearAlgebra>
110 void
111 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
112 setAlpha( Scalar _alpha )
113 {
114  ASSERT( _alpha >= 0.0 );
115  alpha = _alpha;
116  // Building alpha_Id0
117  alpha_Id0 = _alpha * calculus.template identity<0, PRIMAL>();
118  alpha_g0.clear();
119  for ( unsigned int i = 0; i < g0.size(); i++ )
120  alpha_g0.push_back( alpha_Id0 * g0[ i ] );
121 }
122 
123 //-----------------------------------------------------------------------------
124 template <typename TKSpace, typename TLinearAlgebra>
125 void
126 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
127 setAlpha( Scalar _alpha, const PrimalForm0& m )
128 {
129  ASSERT( _alpha >= 0.0 );
130  alpha = _alpha;
131  // Building alpha_Id0
132  alpha_Id0 = _alpha * functions::dec::diagonal( m ) * calculus.template identity<0, PRIMAL>();
133  alpha_g0.clear();
134  for ( unsigned int i = 0; i < g0.size(); i++ )
135  alpha_g0.push_back( alpha_Id0 * g0[ i ] );
136 }
137 
138 //-----------------------------------------------------------------------------
139 template <typename TKSpace, typename TLinearAlgebra>
140 void
141 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
142 setLambda( Scalar _lambda )
143 {
144  ASSERT( _lambda >= 0.0 );
145  lambda = _lambda;
146  l_L1 = lambda * L1;
147  l_1_over_4 = (lambda / 4.0 ) * KForm< Calculus, 1, PRIMAL >::ones( calculus );
148 }
149 
150 //-----------------------------------------------------------------------------
151 template <typename TKSpace, typename TLinearAlgebra>
152 void
153 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
154 setEpsilon( Scalar _epsilon )
155 {
156  ASSERT( _epsilon > 0.0 );
157  epsilon = _epsilon;
158  left_V1 = epsilon * l_L1 +
159  ( lambda/(4.0*epsilon) ) * calculus.template identity<1, PRIMAL>();
160  l_1_over_4e = (1.0/epsilon) * l_1_over_4;
161 }
162 
163 //-----------------------------------------------------------------------------
164 template <typename TKSpace, typename TLinearAlgebra>
165 void
166 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
167 setUFromInput()
168 {
169  u0 = g0;
170 }
171 
172 //-----------------------------------------------------------------------------
173 template <typename TKSpace, typename TLinearAlgebra>
174 bool
175 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
176 solveU()
177 {
178  if ( verbose > 0 ) trace.beginBlock("Solving for u");
179  if ( verbose > 1 ) trace.info() << "- building matrix M : = alpha_Id0 - tA_Diag(v)^2_A" << std::endl;
180 
181  PrimalIdentity1 diag_vv = functions::dec::squaredDiagonal( v1 );
182  PrimalIdentity0 Av2A = D0.transpose() * diag_vv * D0 + alpha_Id0;
183  if ( verbose > 1 ) trace.info() << "- prefactoring matrix M" << std::endl;
184  solver_u.compute( Av2A );
185  bool ok = true;
186  for ( Dimension i = 0; i < u0.size(); ++i )
187  {
188  if ( verbose > 1 ) trace.info() << "- solving M u[" << i << "] = alpha g[" << i << "]" << std::endl;
189  u0[ i ] = solver_u.solve( alpha_g0[ i ] );
190  if ( verbose > 1 ) trace.info() << ( solver_u.isValid() ? "=> OK" : "ERROR" ) << " " << solver_u.myLinearAlgebraSolver.info() << std::endl;
191  ok = ok && solver_u.isValid();
192  }
193  if ( verbose > 0 ) trace.endBlock();
194  return ok;
195 }
196 
197 //-----------------------------------------------------------------------------
198 template <typename TKSpace, typename TLinearAlgebra>
199 bool
200 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
201 solveV()
202 {
203  former_v1 = v1;
204  if ( verbose > 0 ) trace.beginBlock("Solving for v");
205  if ( verbose > 1 ) trace.info() << "- building matrix N := l/4e Id1 + le (tA' A' + tB B) + sum Diag(Au_i)^2" << std::endl;
206  // Not the same result as below ! Seems the "good" way of doing it.
207  PrimalIdentity1 N = left_V1;
208  for ( Dimension i = 0; i < u0.size(); ++i )
209  {
210  PrimalIdentity1 U2 = functions::dec::squaredDiagonal( D0 * u0[ i ] );
211  N.myContainer += U2.myContainer;
212  }
213  // Seems the "bad" way of doing it.
214  // PrimalIdentity0 U2 = functions::dec::squaredDiagonal( u0[ 0 ] );
215  // for ( Dimension i = 1; i < u0.size(); ++i )
216  // U2.myContainer += functions::dec::squaredDiagonal( u0[ i ] ).myContainer;
217  // PrimalIdentity1 N = left_V1 + D0 * U2 * D0.transpose();
218  typedef typename PrimalIdentity1::Container Matrix;
219  const Matrix & M = N.myContainer;
220  if ( verbose > 2 )
221  for (int k = 0; k < M.outerSize(); ++k)
222  for ( typename Matrix::InnerIterator it( M, k ); it; ++it )
223  if ( ( verbose > 3 ) || ( it.row() == it.col() ) )
224  trace.info() << "[" << it.row() << "," << it.col() << "] = " << it.value() << std::endl;
225  if ( verbose > 1 ) trace.info() << "- prefactoring matrix N" << std::endl;
226  solver_v.compute( N );
227  if ( verbose > 1 ) trace.info() << "- solving N v = l/4e 1" << std::endl;
228  v1 = solver_v.solve( l_1_over_4e );
229  if ( verbose > 1 ) trace.info() << ( solver_v.isValid() ? "OK" : "ERROR" )
230  << " " << solver_v.myLinearAlgebraSolver.info()
231  << std::endl;
232  if ( verbose > 0 ) trace.endBlock();
233  return solver_v.isValid();
234 }
235 
236 //-----------------------------------------------------------------------------
237 template <typename TKSpace, typename TLinearAlgebra>
238 typename DGtal::ATu0v1<TKSpace, TLinearAlgebra>::Scalar
239 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
240 computeVariation()
241 {
242  if ( verbose > 0 ) trace.beginBlock( "Compute variation of v.");
243  delta_v_l1 = 0.0;
244  delta_v_l2 = 0.0;
245  delta_v_loo = 0.0;
246  for ( Index index = 0; index < size1(); index++)
247  {
248  delta_v_loo = std::max( delta_v_loo, std::fabs( v1.myContainer( index )
249  - former_v1.myContainer( index ) ) );
250  delta_v_l2 += ( v1.myContainer( index ) - former_v1.myContainer( index ) )
251  * ( v1.myContainer( index ) - former_v1.myContainer( index ) );
252  delta_v_l1 += fabs( v1.myContainer( index )
253  - former_v1.myContainer( index ) );
254  }
255  delta_v_l1 /= size1();
256  delta_v_l2 = sqrt( delta_v_l2 / size1() );
257  if ( verbose > 0 ) {
258  trace.info() << "Variation |v^k+1 - v^k|_oo = " << delta_v_loo << std::endl;
259  trace.info() << "Variation |v^k+1 - v^k|_2 = " << delta_v_l2 << std::endl;
260  trace.info() << "Variation |v^k+1 - v^k|_1 = " << delta_v_l1 << std::endl;
261  }
262  if ( verbose > 0 ) trace.endBlock();
263  return delta_v_loo;
264 }
265 
266 //-----------------------------------------------------------------------------
267 template <typename TKSpace, typename TLinearAlgebra>
268 typename DGtal::ATu0v1<TKSpace, TLinearAlgebra>::Scalar
269 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::
270 checkV()
271 {
272  if ( verbose > 0 ) trace.beginBlock("Checking v");
273  Scalar m1 = 1.0;
274  Scalar m2 = 0.0;
275  Scalar ma = 0.0;
276  for ( Index index = 0; index < size1(); index++)
277  {
278  Scalar val = v1.myContainer( index );
279  m1 = std::min( m1, val );
280  m2 = std::max( m2, val );
281  ma += val;
282  }
283  if ( verbose > 0 )
284  trace.info() << "1-form v: min=" << m1 << " avg=" << ( ma / size1() )
285  << " max=" << m2 << std::endl;
286  // for ( Index index = 0; index < size1(); index++)
287  // v1.myContainer( index ) = std::min( std::max(v1.myContainer( index ), 0.0) , 1.0 );
288  if ( verbose > 0 ) trace.endBlock();
289  return std::max( std::fabs( m1 ), std::fabs( m2 - 1.0 ) );
290 }
291 
292 
293 ///////////////////////////////////////////////////////////////////////////////
294 // Interface - public :
295 
296 /**
297  * Writes/Displays the object on an output stream.
298  * @param out the output stream where the object is written.
299  */
300 template <typename TKSpace, typename TLinearAlgebra>
301 inline
302 void
303 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::selfDisplay ( std::ostream & out ) const
304 {
305  out << "[ ATu0v1 #g=" << g0.size() << " dec=" << calculus << " ]";
306 }
307 
308 /**
309  * Checks the validity/consistency of the object.
310  * @return 'true' if the object is valid, 'false' otherwise.
311  */
312 template <typename TKSpace, typename TLinearAlgebra>
313 inline
314 bool
315 DGtal::ATu0v1<TKSpace, TLinearAlgebra>::isValid() const
316 {
317  return true;
318 }
319 
320 
321 
322 ///////////////////////////////////////////////////////////////////////////////
323 // Implementation of inline functions //
324 
325 template <typename TKSpace, typename TLinearAlgebra>
326 inline
327 std::ostream&
328 DGtal::operator<< ( std::ostream & out,
329  const ATu0v1<TKSpace, TLinearAlgebra> & object )
330 {
331  object.selfDisplay( out );
332  return out;
333 }
334 
335 // //
336 ///////////////////////////////////////////////////////////////////////////////
337 
338 
Companion project associated to DGtal ; Generated on Mon Aug 26 2024 21:31:30 for DGtalTools by   doxygen 1.9.1