DGtal  1.5.beta
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  */
50 inline
51 double
52 DGtal::AngleLinearMinimizer::sum() const
53 {
54  return mySum;
55 }
56 
57 /**
58  * Max of all the absolute displacements of the last optimisation step.
59  */
60 inline
61 double
62 DGtal::AngleLinearMinimizer::max() const
63 {
64  return myMax;
65 }
66 
67 
68 /**
69  * Default constructor. Does nothing.
70  */
71 inline
72 DGtal::AngleLinearMinimizerByRelaxation::AngleLinearMinimizerByRelaxation()
73 {}
74 
75 
76 /**
77  * Destructor. Does nothing.
78  */
79 inline
80 DGtal::AngleLinearMinimizerByRelaxation::~AngleLinearMinimizerByRelaxation()
81 {}
82 
83 
84 /**
85  * Default constructor. Does nothing.
86  */
87 inline
88 DGtal::AngleLinearMinimizerByGradientDescent::AngleLinearMinimizerByGradientDescent
89 ( double aStep )
90  : myStep( aStep )
91 {}
92 
93 
94 /**
95  * Destructor. Does nothing.
96  */
97 inline
98 DGtal::AngleLinearMinimizerByGradientDescent::~AngleLinearMinimizerByGradientDescent()
99 {}
100 
101 
102 /**
103  * Default constructor. Does nothing.
104  */
105 inline
106 DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::AngleLinearMinimizerByAdaptiveStepGradientDescent
107 ( double aStep )
108  : myStep( aStep )
109 {}
110 
111 /**
112  * Destructor. Does nothing.
113  */
114 inline
115 DGtal::AngleLinearMinimizerByAdaptiveStepGradientDescent::~AngleLinearMinimizerByAdaptiveStepGradientDescent()
116 {}
117 
118 
119 
120 /**
121  * @return a reference on the information structure of the [i]th value.
122  */
123 inline
124 DGtal::AngleLinearMinimizer::ValueInfo &
125 DGtal::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  */
135 inline
136 const DGtal::AngleLinearMinimizer::ValueInfo &
137 DGtal::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  */
147 inline
148 unsigned int
149 DGtal::AngleLinearMinimizer::maxSize() const
150 {
151  return myMaxSize;
152 }
153 
154 
155 /**
156  * @return the number of values stored in the object.
157  */
158 inline
159 unsigned int
160 DGtal::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  */
170 inline
171 void
172 DGtal::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  */
185 inline
186 void
187 DGtal::AngleLinearMinimizer::setIsCurveOpen( bool isCurveOpen)
188 {
189  myIsCurveOpen = isCurveOpen;
190 }
191 
192 inline
193 DGtal::AngleLinearMinimizer::~AngleLinearMinimizer()
194 {
195  reset();
196 }
197 
198 
199 inline
200 DGtal::AngleLinearMinimizer::AngleLinearMinimizer()
201  : myValues( 0 )
202 {
203 }
204 
205 
206 inline
207 void
208 DGtal::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 
220 inline
221 void
222 DGtal::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 
237 inline
238 double
239 DGtal::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 
257 inline
258 double
259 DGtal::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 
278 inline
279 std::vector<double>
280 DGtal::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 
319 inline
320 std::vector<double>
321 DGtal::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 
358 inline
359 double
360 DGtal::AngleLinearMinimizer::optimize()
361 {
362  return optimize( 0, 0 );
363 }
364 
365 
366 inline
367 double
368 DGtal::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 
403 inline
404 double
405 DGtal::AngleLinearMinimizer::lastDelta() const
406 {
407  return max();
408 }
409 
410 
411 inline
412 void
413 DGtal::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 
463 inline
464 void
465 DGtal::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 
512 inline
513 double
514 DGtal::AngleLinearMinimizerByRelaxation::lastDelta() const
515 {
516  return max();
517 }
518 
519 
520 
521 inline
522 void
523 DGtal::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 
547 inline
548 double
549 DGtal::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 
568 inline
569 void
570 DGtal::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 
606 inline
607 double
608 DGtal::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  */
634 inline
635 void
636 DGtal::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  */
646 inline
647 void
648 DGtal::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  */
657 inline
658 void
659 DGtal::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  */
668 inline
669 void
670 DGtal::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  */
679 inline
680 bool
681 DGtal::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 
700 inline
701 std::string
702 DGtal::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  */
713 inline
714 std::ostream&
715 DGtal::operator<< ( std::ostream & out,
716  const AngleLinearMinimizer & object )
717 {
718  object.selfDisplay ( out );
719  return out;
720 }
721 
722 // //
723 ///////////////////////////////////////////////////////////////////////////////
724 
725