DGtal 1.4.0
Loading...
Searching...
No Matches
AngleLinearMinimizer.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 AngleLinearMinimizer.ih
19
20 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
21 * Laboratory of Mathematics (CNRS, UMR 5807), University of Savoie, France
22 *
23 * @author (backported from ImaGene by) Bertrand Kerautret (\c kerautre@loria.fr )
24 * LORIA (CNRS, UMR 7503), University of Nancy, France
25 *
26 * @date 2011/08/31
27 *
28 * Implementation of inline methods defined in AngleLinearMinimizer.h
29 *
30 * This file is part of the DGtal library.
31 */
32
33///////////////////////////////////////////////////////////////////////////////
34// IMPLEMENTATION of inline methods.
35///////////////////////////////////////////////////////////////////////////////
36
37//////////////////////////////////////////////////////////////////////////////
38#include <cstdlib>
39//////////////////////////////////////////////////////////////////////////////
40
41
42
43///////////////////////////////////////////////////////////////////////////////
44// Implementation of inline methods //
45
46
47/**
48 * Sum of all the absolute displacements of the last optimisation step.
49 */
50inline
51double
52DGtal::AngleLinearMinimizer::sum() const
53{
54 return mySum;
55}
56
57/**
58 * Max of all the absolute displacements of the last optimisation step.
59 */
60inline
61double
62DGtal::AngleLinearMinimizer::max() const
63{
64 return myMax;
65}
66
67
68/**
69 * Default constructor. Does nothing.
70 */
71inline
72DGtal::AngleLinearMinimizerByRelaxation::AngleLinearMinimizerByRelaxation()
73{}
74
75
76/**
77 * Destructor. Does nothing.
78 */
79inline
80DGtal::AngleLinearMinimizerByRelaxation::~AngleLinearMinimizerByRelaxation()
81{}
82
83
84/**
85 * Default constructor. Does nothing.
86 */
87inline
88DGtal::AngleLinearMinimizerByGradientDescent::AngleLinearMinimizerByGradientDescent
89( double aStep )
90 : myStep( aStep )
91{}
92
93
94/**
95 * Destructor. Does nothing.
96 */
97inline
98DGtal::AngleLinearMinimizerByGradientDescent::~AngleLinearMinimizerByGradientDescent()
99{}
100
101
102/**
103 * Default constructor. Does nothing.
104 */
105inline
106DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::AngleLinearMinimizerByAdaptiveStepGradientDescent
107( double aStep )
108 : myStep( aStep )
109{}
110
111/**
112 * Destructor. Does nothing.
113 */
114inline
115DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::~AngleLinearMinimizerByAdaptiveStepGradientDescent()
116{}
117
118
119
120/**
121 * @return a reference on the information structure of the [i]th value.
122 */
123inline
124DGtal::AngleLinearMinimizer::ValueInfo &
125DGtal::AngleLinearMinimizer::rw( unsigned int i )
126{
127 ASSERT( ( myValues != 0 ) && ( i < maxSize() ) );
128 return myValues[ i ];
129}
130
131
132/**
133 * @return a const reference on the information structure of the [i]th value.
134 */
135inline
136const DGtal::AngleLinearMinimizer::ValueInfo &
137DGtal::AngleLinearMinimizer::ro( unsigned int i ) const
138{
139 ASSERT( ( myValues != 0 ) && ( i < maxSize() ) );
140 return myValues[ i ];
141}
142
143
144/**
145 * @return the maximum number of values stored in the object.
146 */
147inline
148unsigned int
149DGtal::AngleLinearMinimizer::maxSize() const
150{
151 return myMaxSize;
152}
153
154
155/**
156 * @return the number of values stored in the object.
157 */
158inline
159unsigned int
160DGtal::AngleLinearMinimizer::size() const
161{
162 return mySize;
163}
164
165
166/**
167 * Specifies the exact number of valid values.
168 * @param nb any number below 'maxSize()'.
169 */
170inline
171void
172DGtal::AngleLinearMinimizer::setSize( unsigned int nb )
173{
174 ASSERT( nb <= maxSize() );
175 mySize = nb;
176}
177
178
179/**
180 * Specifies if the curve is open or not.
181 * @param isCurveOpen when 'true' the curve is open and the last
182 * value does not depend on the first one, otherwise the curve is
183 * closed and the last value is linked to the first one.
184 */
185inline
186void
187DGtal::AngleLinearMinimizer::setIsCurveOpen( bool isCurveOpen)
188{
189 myIsCurveOpen = isCurveOpen;
190}
191
192inline
193DGtal::AngleLinearMinimizer::~AngleLinearMinimizer()
194{
195 reset();
196}
197
198
199inline
200DGtal::AngleLinearMinimizer::AngleLinearMinimizer()
201 : myValues( 0 )
202{
203}
204
205
206inline
207void
208DGtal::AngleLinearMinimizer::reset()
209{
210 if ( myValues != 0 )
211 {
212 delete[] myValues;
213 myValues = 0;
214 }
215 mySum = 0.0;
216 myMax = 0.0;
217}
218
219
220inline
221void
222DGtal::AngleLinearMinimizer::init( unsigned int nbMax )
223{
224 reset();
225 myValues = new ValueInfo[ nbMax ];
226 myMaxSize = nbMax;
227 mySize = nbMax;
228 myIsCurveOpen = false;
229}
230
231
232
233
234///////////////////////////////////////////////////////////////////////////////
235// ------------------------- Optimization services --------------------------
236
237inline
238double
239DGtal::AngleLinearMinimizer::getEnergy( unsigned int i1, unsigned int i2 ) const
240{
241 ModuloComputer<int> mc( size() );
242 double E = 0.0;
243 for ( unsigned int i = mc.next( i1 ); i != i2; )
244 {
245 unsigned int inext = mc.next( i );
246 const ValueInfo & vi = this->ro( i );
247 const ValueInfo & viprev = this->ro( mc.previous( i ) );
248 double dev = AngleComputer::deviation( vi.value, viprev.value );
249 E += (dev * dev) / viprev.distToNext;
250 i = inext;
251 }
252 return E;
253}
254
255
256
257inline
258double
259DGtal::AngleLinearMinimizer::getFormerEnergy( unsigned int i1, unsigned int i2 ) const
260{
261 ModuloComputer<int> mc( size() );
262 double E = 0.0;
263 for ( unsigned int i = mc.next( i1 ); i != i2; )
264 {
265 unsigned int inext = mc.next( i );
266 const ValueInfo & vi = this->ro( i );
267 const ValueInfo & viprev = this->ro( mc.previous( i ) );
268 double dev = AngleComputer::deviation( vi.oldValue, viprev.oldValue );
269 E += (dev * dev) / viprev.distToNext;
270 i = inext;
271 }
272 return E;
273}
274
275
276
277
278inline
279std::vector<double>
280DGtal::AngleLinearMinimizer::getGradient() const
281{
282 ModuloComputer<int> mc( size() );
283
284 std::vector<double> grad( size() );
285 for ( unsigned int i = 0; i < size(); i++ )
286 {
287 unsigned int iprev = mc.previous( i );
288 unsigned int inext = mc.next( i );
289 const ValueInfo & vi = this->ro( i );
290 double val = vi.value;
291 if ( myIsCurveOpen && ( i == ( size() - 1 ) ) )
292 { // free extremity to the front/right.
293 const ValueInfo & viprev = this->ro( iprev );
294 double valp = viprev.value;
295 grad[ i ] = 2.0 * AngleComputer::deviation( val, valp ) / viprev.distToNext;
296 }
297 else if ( myIsCurveOpen && ( i == 0 ) )
298 { // free extremity to the back/left.
299 const ValueInfo & vinext = this->ro( inext );
300 double valn = vinext.value;
301 grad[ i ] = -2.0 * AngleComputer::deviation( valn, val ) / vi.distToNext;
302 }
303 else
304 { // standard case.
305 const ValueInfo & viprev = this->ro( iprev );
306 const ValueInfo & vinext = this->ro( inext );
307 double valp = viprev.value;
308 double valn = vinext.value;
309 grad[ i ] = 2.0*( AngleComputer::deviation( val, valp ) / viprev.distToNext
310 - AngleComputer::deviation( valn, val ) / vi.distToNext ) ;
311 }
312 }
313 return grad;
314}
315
316
317
318
319inline
320std::vector<double>
321DGtal::AngleLinearMinimizer::getFormerGradient() const
322{
323 ModuloComputer< int> mc( size() );
324
325 std::vector<double> grad( size() );
326 for ( unsigned int i = 0; i < size(); i++ )
327 {
328 unsigned int iprev = mc.previous( i );
329 unsigned int inext = mc.next( i );
330 const ValueInfo & vi = this->ro( i );
331 double val = vi.oldValue;
332 if ( myIsCurveOpen && ( i == ( size() - 1 ) ) )
333 { // free extremity to the front/right.
334 const ValueInfo & viprev = this->ro( iprev );
335 double valp = viprev.oldValue;
336 grad[ i ] = 2.0 * AngleComputer::deviation( val, valp ) / viprev.distToNext;
337 }
338 else if ( myIsCurveOpen && ( i == 0 ) )
339 { // free extremity to the back/left.
340 const ValueInfo & vinext = this->ro( inext );
341 double valn = vinext.oldValue;
342 grad[ i ] = -2.0 * AngleComputer::deviation( valn, val ) / vi.distToNext;
343 }
344 else
345 { // standard case.
346 const ValueInfo & viprev = this->ro( iprev );
347 const ValueInfo & vinext = this->ro( inext );
348 double valp = viprev.oldValue;
349 double valn = vinext.oldValue;
350 grad[ i ] = 2.0*( AngleComputer::deviation( val, valp ) / viprev.distToNext
351 - AngleComputer::deviation( valn, val ) / vi.distToNext ) ;
352 }
353 }
354 return grad;
355}
356
357
358inline
359double
360DGtal::AngleLinearMinimizer::optimize()
361{
362 return optimize( 0, 0 );
363}
364
365
366inline
367double
368DGtal::AngleLinearMinimizer::optimize( unsigned int i1, unsigned int i2 )
369{
370 ASSERT( size() > 2 );
371 ModuloComputer< int> mc( size() );
372
373 unsigned int i = i1;
374 // (I) first pass to work on old values.
375 do
376 {
377 ValueInfo & vi = this->rw( i );
378 vi.oldValue = vi.value;
379 // go to next.
380 i = mc.next( i );
381 }
382 while ( i != i2 );
383 this->oneStep( i1, i2 );
384
385 mySum = 0.0;
386 myMax = 0.0;
387 i = i1;
388 do
389 {
390 const ValueInfo & vi = this->ro( i );
391 double diff = fabs( AngleComputer::deviation( vi.value, vi.oldValue ) );
392 if ( diff > myMax )
393 myMax = diff;
394 mySum += diff;
395 i = mc.next( i );
396 }
397 while ( i != i2 );
398
399 return mySum;
400}
401
402
403inline
404double
405DGtal::AngleLinearMinimizer::lastDelta() const
406{
407 return max();
408}
409
410
411inline
412void
413DGtal::AngleLinearMinimizer::oneStep( unsigned int i1, unsigned int i2 )
414{
415 ModuloComputer< int> mc( size() );
416
417 double mid;
418 unsigned int i = i1;
419 unsigned int iprev = mc.previous( i );
420 do
421 {
422 unsigned int inext = mc.next( i );
423 ValueInfo & vi = this->rw( i );
424 if ( myIsCurveOpen && ( i == ( size() - 1 ) ) )
425 { // free extremity to the front/right.
426 const ValueInfo & viprev = this->ro( iprev );
427 mid = viprev.oldValue; // - viprev.delta;
428 }
429 else if ( myIsCurveOpen && ( i == 0 ) )
430 { // free extremity to the back/left.
431 const ValueInfo & vinext = this->ro( inext );
432 mid = vinext.oldValue; // + vi.delta;
433 }
434 else
435 { // standard case.
436 const ValueInfo & viprev = this->ro( iprev );
437 const ValueInfo & vinext = this->ro( inext );
438 double valp = viprev.oldValue; // - viprev.delta;
439 double valn = vinext.oldValue; // + vi.delta;
440
441 // old
442 double y = AngleComputer::deviation( valn, valp );
443 mid = ( viprev.distToNext * y )
444 / ( vi.distToNext + viprev.distToNext );
445 mid = AngleComputer::cast( mid + valp );
446
447 }
448 if ( AngleComputer::less( mid, vi.min ) ) mid = vi.min;
449 if ( AngleComputer::less( vi.max, mid ) ) mid = vi.max;
450 double mvt = AngleComputer::deviation( mid, vi.oldValue );
451 vi.value = AngleComputer::cast( vi.oldValue + 0.5 * mvt );
452 // go to next.
453 iprev = i;
454 i = inext;
455 }
456 while ( i != i2 );
457
458
459}
460
461
462
463inline
464void
465DGtal::AngleLinearMinimizerByRelaxation::oneStep( unsigned int i1, unsigned int i2 )
466{
467 ModuloComputer< int> mc( size() );
468
469 double mid;
470 unsigned int i = i1;
471 unsigned int iprev = mc.previous( i );
472 do
473 {
474 unsigned int inext = mc.next( i );
475 ValueInfo & vi = this->rw( i );
476 // vi.oldValue = vi.value;
477 if ( myIsCurveOpen && ( i == ( size() - 1 ) ) )
478 { // free extremity to the front/right.
479 const ValueInfo & viprev = this->ro( iprev );
480 mid = viprev.value; // - viprev.delta;
481 }
482 else if ( myIsCurveOpen && ( i == 0 ) )
483 { // free extremity to the back/left.
484 const ValueInfo & vinext = this->ro( inext );
485 mid = vinext.oldValue; // + vi.delta;
486 }
487 else
488 { // standard case.
489 const ValueInfo & viprev = this->ro( iprev );
490 const ValueInfo & vinext = this->ro( inext );
491 double valp = viprev.value; // - viprev.delta;
492 double valn = vinext.value;
493
494 // old
495 double y = AngleComputer::deviation( valn, valp );
496 mid = ( viprev.distToNext * y )
497 / ( vi.distToNext + viprev.distToNext );
498 mid = AngleComputer::cast( mid + valp );
499 }
500 if ( AngleComputer::less( mid, vi.min ) ) mid = vi.min;
501 if ( AngleComputer::less( vi.max, mid ) ) mid = vi.max;
502 vi.value = mid;
503 iprev = i;
504 i = inext;
505 }
506 while ( i != i2 );
507}
508
509
510
511
512inline
513double
514DGtal::AngleLinearMinimizerByRelaxation::lastDelta() const
515{
516 return max();
517}
518
519
520
521inline
522void
523DGtal::AngleLinearMinimizerByGradientDescent::oneStep( unsigned int i1, unsigned int i2 )
524{
525 ModuloComputer< int> mc( size() );
526
527 std::vector<double> grad ( getFormerGradient() );
528 double mid;
529 unsigned int i = i1;
530 do
531 {
532 unsigned int inext = mc.next( i );
533 ValueInfo & vi = this->rw( i );
534 double val = vi.oldValue;
535 mid = AngleComputer::cast( val - myStep*grad[ i ] );
536 if ( AngleComputer::less( mid, vi.min ) ) mid = vi.min;
537 if ( AngleComputer::less( vi.max, mid ) ) mid = vi.max;
538 vi.value = mid;
539 // go to next.
540 i = inext;
541 }
542 while ( i != i2 );
543}
544
545
546
547inline
548double
549DGtal::AngleLinearMinimizerByGradientDescent::lastDelta() const
550{
551 std::vector<double> grad ( getFormerGradient() );
552 double ninf = 0.0;
553 for ( unsigned int i = 0; i < size(); i++ )
554 {
555 const ValueInfo & vi = this->ro( i );
556 if ( vi.value != vi.oldValue )
557 {
558 double n = fabs( grad[ i ] );
559 if ( n > ninf ) ninf = n;
560 }
561 }
562 return ninf;
563}
564
565
566
567
568inline
569void
570DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::oneStep( unsigned int i1, unsigned int i2 )
571{
572 ModuloComputer< int> mc( size() );
573 std::vector<double> grad ( getFormerGradient() );
574
575 double mid;
576 unsigned int i = i1;
577 do
578 {
579 unsigned int inext = mc.next( i );
580 ValueInfo & vi = this->rw( i );
581 double val = vi.oldValue;
582 mid = AngleComputer::cast( val - myStep*grad[ i ] );
583 if ( AngleComputer::less( mid, vi.min ) ) mid = vi.min;
584 if ( AngleComputer::less( vi.max, mid ) ) mid = vi.max;
585 vi.value = mid;
586 // go to next.
587 i = inext;
588 }
589 while ( i != i2 );
590
591 double E1 = getFormerEnergy( i1, i2 );
592 double E2 = getEnergy( i1, i2 );
593 if ( E1 <= E2 )
594 {
595 myStep /= 4.0;
596 }
597 else
598 {
599 /* doubling step. */
600 myStep *= 2.0;
601 }
602}
603
604
605
606inline
607double
608DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::lastDelta() const
609{
610 std::vector<double> grad ( getFormerGradient() );
611 double ninf = 0.0;
612 for ( unsigned int i = 0; i < size(); i++ )
613 {
614 const ValueInfo & vi = this->ro( i );
615 if ( vi.value != vi.oldValue )
616 {
617 double n = fabs( grad[ i ] );
618 if ( n > ninf ) ninf = n;
619 }
620 }
621 return ninf;
622}
623
624
625
626
627///////////////////////////////////////////////////////////////////////////////
628// Interface - public :
629
630/**
631 * Writes/Displays the object on an output stream.
632 * @param aStream the output stream where the object is written.
633 */
634inline
635void
636DGtal::AngleLinearMinimizer::selfDisplay ( std::ostream & aStream ) const
637{
638 aStream << "[AngleLinearMinimizer::standard 0.5]";
639}
640
641
642/**
643 * Writes/Displays the object on an output stream.
644 * @param aStream the output stream where the object is written.
645 */
646inline
647void
648DGtal::AngleLinearMinimizerByRelaxation::selfDisplay( std::ostream & aStream ) const
649{
650 aStream << "[LinearMinimizer::relaxation]";
651}
652
653/**
654 * Writes/Displays the object on an output stream.
655 * @param aStream the output stream where the object is written.
656 */
657inline
658void
659DGtal::AngleLinearMinimizerByGradientDescent::selfDisplay( std::ostream& aStream ) const
660{
661 aStream << "[LinearMinimizer::gradient descent " << myStep << "]";
662}
663
664/**
665 * Writes/Displays the object on an output stream.
666 * @param aStream the output stream where the object is written.
667 */
668inline
669void
670DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::selfDisplay( std::ostream& aStream ) const
671{
672 aStream << "[LinearMinimizer::gradient descent with adaptive step " << myStep << "]";
673}
674
675/**
676 * Checks the validity/consistency of the object.
677 * @return 'true' if the object is valid, 'false' otherwise.
678 */
679inline
680bool
681DGtal::AngleLinearMinimizer::isValid() const
682{
683 return true;
684}
685
686
687
688///////////////////////////////////////////////////////////////////////////////
689// Internals - private :
690
691// //
692///////////////////////////////////////////////////////////////////////////////
693
694
695
696
697///////////////////////////////////////////////////////////////////////////////
698// Implementation of inline functions and external operators //
699
700inline
701std::string
702DGtal::AngleLinearMinimizer::className() const
703{
704 return "AngleLinearMinimizer";
705}
706
707/**
708 * Overloads 'operator<<' for displaying objects of class 'AngleLinearMinimizer'.
709 * @param out the output stream where the object is written.
710 * @param object the object of class 'AngleLinearMinimizer' to write.
711 * @return the output stream after the writing.
712 */
713inline
714std::ostream&
715DGtal::operator<< ( std::ostream & out,
716 const AngleLinearMinimizer & object )
717{
718 object.selfDisplay ( out );
719 return out;
720}
721
722// //
723///////////////////////////////////////////////////////////////////////////////
724
725