DGtalTools 2.0.0
Loading...
Searching...
No Matches
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 ------------------------------
42template <typename TKSpace, typename TLinearAlgebra>
43DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
44ATu2v0( 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
54template <typename TKSpace, typename TLinearAlgebra>
55void
56DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
57init( 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//-----------------------------------------------------------------------------
75template <typename TKSpace, typename TLinearAlgebra>
76template <typename Image, typename Function>
77void
78DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
79addInput( 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//-----------------------------------------------------------------------------
95template <typename TKSpace, typename TLinearAlgebra>
96typename DGtal::ATu2v0<TKSpace, TLinearAlgebra>::Scalar
97DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
98computeSNR() 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//-----------------------------------------------------------------------------
114template <typename TKSpace, typename TLinearAlgebra>
115void
116DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
117setMetricAverage( bool average )
118{
119 metric_average = average;
120}
121
122//-----------------------------------------------------------------------------
123template <typename TKSpace, typename TLinearAlgebra>
124void
125DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
126setAlpha( 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//-----------------------------------------------------------------------------
138template <typename TKSpace, typename TLinearAlgebra>
139void
140DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
141setAlpha( 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//-----------------------------------------------------------------------------
153template <typename TKSpace, typename TLinearAlgebra>
154void
155DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
156setLambda( 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//-----------------------------------------------------------------------------
171template <typename TKSpace, typename TLinearAlgebra>
172void
173DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
174setEpsilon( 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//-----------------------------------------------------------------------------
191template <typename TKSpace, typename TLinearAlgebra>
192void
193DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
194setUFromInput()
195{
196 u2 = g2;
197}
198
199//-----------------------------------------------------------------------------
200template <typename TKSpace, typename TLinearAlgebra>
201void
202DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
203setUFromInputAndMask()
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//-----------------------------------------------------------------------------
218template <typename TKSpace, typename TLinearAlgebra>
219bool
220DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
221solveU()
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//-----------------------------------------------------------------------------
245template <typename TKSpace, typename TLinearAlgebra>
246bool
247DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
248solveV()
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//-----------------------------------------------------------------------------
280template <typename TKSpace, typename TLinearAlgebra>
281typename DGtal::ATu2v0<TKSpace, TLinearAlgebra>::Scalar
282DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
283computeVariation()
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//-----------------------------------------------------------------------------
310template <typename TKSpace, typename TLinearAlgebra>
311typename DGtal::ATu2v0<TKSpace, TLinearAlgebra>::Scalar
312DGtal::ATu2v0<TKSpace, TLinearAlgebra>::
313checkV()
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 */
343template <typename TKSpace, typename TLinearAlgebra>
344inline
345void
346DGtal::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 */
355template <typename TKSpace, typename TLinearAlgebra>
356inline
357bool
358DGtal::ATu2v0<TKSpace, TLinearAlgebra>::isValid() const
359{
360 return true;
361}
362
363
364
365///////////////////////////////////////////////////////////////////////////////
366// Implementation of inline functions //
367
368template <typename TKSpace, typename TLinearAlgebra>
369inline
370std::ostream&
371DGtal::operator<< ( std::ostream & out,
372 const ATu2v0<TKSpace, TLinearAlgebra> & object )
373{
374 object.selfDisplay( out );
375 return out;
376}
377
378// //
379///////////////////////////////////////////////////////////////////////////////
380
381