DGtal  1.2.0
FMMPointFunctors.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 FMMPointFunctors.ih
19  * @author Tristan Roussillon (\c
20  * tristan.roussillon@liris.cnrs.fr ) Laboratoire d'InfoRmatique en
21  * Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS,
22  * France
23  *
24  *
25  * @date 2012/02/21
26  *
27  * @brief Implementation of inline methods defined in FMMPointFunctors.h
28  *
29  * This file is part of the DGtal library.
30  */
31 
32 
33 //////////////////////////////////////////////////////////////////////////////
34 #include <cstdlib>
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // Helper functions defined in the compilation unit (anonymous namespace)
39 // see operator() below
40 namespace
41 {
42  //comparator in absolute value
43  template <typename Value>
44  bool absComparator(const Value& i, const Value& j)
45  {
46  return ( std::abs(static_cast<double>(i)) < std::abs(static_cast<double>(j)) );
47  }
48 
49  //pair second member comparator in absolute value
50  template <typename Pair>
51  bool secondAbsComparator(const Pair& i, const Pair& j)
52  {
53  return absComparator( i.second, j.second );
54  }
55 
56 }
57 
58 ///////////////////////////////////////////////////////////////////////////////
59 // IMPLEMENTATION of inline methods.
60 ///////////////////////////////////////////////////////////////////////////////
61 
62 ///////////////////////////////////////////////////////////////////////////////
63 //-----------------------------------------------------------------------------
64 template <typename TImage, typename TSet>
65 inline
66 DGtal::L2FirstOrderLocalDistance<TImage,TSet>::L2FirstOrderLocalDistance
67  (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
68 {
69  ASSERT( myImgPtr );
70  ASSERT( mySetPtr );
71 }
72 
73 //-----------------------------------------------------------------------------
74 template <typename TImage, typename TSet>
75 inline
76 DGtal::L2FirstOrderLocalDistance<TImage,TSet>::L2FirstOrderLocalDistance
77  (const L2FirstOrderLocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
78 {
79  ASSERT( myImgPtr );
80  ASSERT( mySetPtr );
81 }
82 
83 //-----------------------------------------------------------------------------
84 template <typename TImage, typename TSet>
85 inline
86 DGtal::L2FirstOrderLocalDistance<TImage,TSet>&
87 DGtal::L2FirstOrderLocalDistance<TImage,TSet>::operator=
88  (const L2FirstOrderLocalDistance& other)
89 {
90  if( this != &other)
91  {
92  myImgPtr = other.myImgPtr;
93  mySetPtr = other.mySetPtr;
94  ASSERT( myImgPtr );
95  ASSERT( mySetPtr );
96  }
97  return *this;
98 }
99 
100 //-----------------------------------------------------------------------------
101 template <typename TImage, typename TSet>
102 inline
103 DGtal::L2FirstOrderLocalDistance<TImage,TSet>::~L2FirstOrderLocalDistance
104  ()
105 {
106 }
107 
108 //-----------------------------------------------------------------------------
109 template <typename TImage, typename TSet>
110 inline
111 typename DGtal::L2FirstOrderLocalDistance<TImage, TSet>::Value
112 DGtal::L2FirstOrderLocalDistance<TImage, TSet>::operator()
113  (const Point& aPoint)
114 {
115 
116  //distance values
117  Values v;
118  v.reserve(Point::dimension);
119 
120  //two 1-neighbors
121  Point neighbor1 = aPoint;
122  Point neighbor2 = aPoint;
123 
124  typename Point::Iterator it1 = neighbor1.begin();
125  typename Point::Iterator it2 = neighbor2.begin();
126  typename Point::ConstIterator it = aPoint.begin();
127  typename Point::ConstIterator itEnd = aPoint.end();
128  for ( ; it != itEnd; ++it, ++it1, ++it2)
129  {//for each dimension
130 
131  typename Point::Coordinate c = *it;
132  *it1 = (c+1);
133  *it2 = (c-1);
134 
135  //neighboring values
136  Value d = 0, d1 = 0, d2 = 0;
137  bool flag1 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor1, d1 );
138  bool flag2 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor2, d2 );
139  if ( flag1 || flag2 )
140  {
141  if ( flag1 && flag2 )
142  { //take the minimal value
143  if (std::abs(d1) < std::abs(d2))
144  d = d1;
145  else
146  d = d2;
147  } else
148  {
149  if (flag1) d = d1;
150  if (flag2) d = d2;
151  }
152 
153  v.push_back(d);
154  }
155 
156  *it1 = c;
157  *it2 = c;
158  } //end for each dimension
159 
160  //computation of the new value
161  return this->compute(v);
162 }
163 
164 
165 //-----------------------------------------------------------------------------
166 template <typename TImage, typename TSet>
167 inline
168 typename DGtal::L2FirstOrderLocalDistance<TImage, TSet>::Value
169 DGtal::L2FirstOrderLocalDistance<TImage, TSet>::compute
170 (Values& aValueList) const
171 {
172  ASSERT(aValueList.size() > 0);
173 
174  if ( aValueList.size() == 1 )
175  {
176  Value d = aValueList.back();
177  if (d >= 0) return d + 1.0;
178  else return d - 1.0;
179  }
180  else
181  {
182  //function computation
183  typename Values::iterator itMax =
184  std::max_element( aValueList.begin(), aValueList.end(), absComparator<Value> );
185  if ( gradientNorm( *itMax, aValueList ) > 1 )
186  {
187  aValueList.erase( itMax );
188  return this->compute(aValueList);
189  }
190  else
191  { //resolution
192  double a = 0;
193  double b = 0;
194  double c = -1;
195 
196  for (typename Values::iterator it = aValueList.begin();
197  it != aValueList.end(); ++it)
198  {
199  Value d = *it;
200 
201  a += 1;
202  b -= static_cast<double>(2*d);
203  c += static_cast<double>(d*d);
204  }
205 
206  //discriminant
207  double disc = b*b - 4*a*c;
208  ASSERT(disc >= 0);
209 
210  if ( b < 0 )
211  return static_cast<Value>( ( -b + std::sqrt(disc) ) / (2*a) );
212  else
213  return static_cast<Value>( ( -b - std::sqrt(disc) ) / (2*a) );
214  }
215  }
216 
217 }
218 
219 //-----------------------------------------------------------------------------
220 template <typename TImage, typename TSet>
221 inline
222 typename DGtal::L2FirstOrderLocalDistance<TImage, TSet>::Value
223 DGtal::L2FirstOrderLocalDistance<TImage, TSet>
224 ::gradientNorm(const Value& aValue, const Values& aValueList) const
225 {
226  Value sum = 0;
227  for (typename Values::const_iterator it = aValueList.begin();
228  it != aValueList.end(); ++it)
229  {
230  Value d = (aValue - *it);
231  sum += (d*d);
232  }
233  return sum;
234 }
235 
236 //-----------------------------------------------------------------------------
237 template <typename TImage, typename TSet>
238 inline
239 void
240 DGtal::L2FirstOrderLocalDistance<TImage, TSet>
241 ::selfDisplay ( std::ostream & out ) const
242 {
243  out << "L2";
244 }
245 
246 ///////////////////////////////////////////////////////////////////////////////
247 //-----------------------------------------------------------------------------
248 template <typename TImage, typename TSet>
249 inline
250 DGtal::L2SecondOrderLocalDistance<TImage,TSet>::L2SecondOrderLocalDistance
251  (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
252 {
253  ASSERT( myImgPtr );
254  ASSERT( mySetPtr );
255 }
256 
257 //-----------------------------------------------------------------------------
258 template <typename TImage, typename TSet>
259 inline
260 DGtal::L2SecondOrderLocalDistance<TImage,TSet>::L2SecondOrderLocalDistance
261  (const L2SecondOrderLocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
262 {
263  ASSERT( myImgPtr );
264  ASSERT( mySetPtr );
265 }
266 
267 //-----------------------------------------------------------------------------
268 template <typename TImage, typename TSet>
269 inline
270 DGtal::L2SecondOrderLocalDistance<TImage,TSet>&
271 DGtal::L2SecondOrderLocalDistance<TImage,TSet>::operator=
272  (const L2SecondOrderLocalDistance& other)
273 {
274  if( this != &other)
275  {
276  myImgPtr = other.myImgPtr;
277  mySetPtr = other.mySetPtr;
278  ASSERT( myImgPtr );
279  ASSERT( mySetPtr );
280  }
281  return *this;
282 }
283 
284 //-----------------------------------------------------------------------------
285 template <typename TImage, typename TSet>
286 inline
287 DGtal::L2SecondOrderLocalDistance<TImage,TSet>::~L2SecondOrderLocalDistance
288  ()
289 {
290 }
291 
292 //-----------------------------------------------------------------------------
293 template <typename TImage, typename TSet>
294 inline
295 typename DGtal::L2SecondOrderLocalDistance<TImage, TSet>::Value
296 DGtal::L2SecondOrderLocalDistance<TImage, TSet>::operator()
297  (const Point& aPoint)
298 {
299 
300  //distance values
301  List l;
302  l.reserve(Point::dimension);
303 
304  //two 1-neighbors
305  Point neighbor1 = aPoint;
306  Point neighbor2 = aPoint;
307 
308  typename Point::Iterator it1 = neighbor1.begin();
309  typename Point::Iterator it2 = neighbor2.begin();
310  typename Point::ConstIterator it = aPoint.begin();
311  typename Point::ConstIterator itEnd = aPoint.end();
312  for ( ; it != itEnd; ++it, ++it1, ++it2)
313  {//for each dimension
314 
315  typename Point::Coordinate c = *it;
316 
317  /// first-order
318  double coeff = 1.0;
319  ++(*it1);
320  --(*it2);
321  Value d = 0, d1 = 0, d2 = 0;
322  bool flag1 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor1, d1 );
323  bool flag2 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor2, d2 );
324  if ( flag1 || flag2 )
325  {
326  /// second-order
327  ++(*it1);
328  --(*it2);
329  Value d12 = 0, d22 = 0;
330  bool flag12 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor1, d12 );
331  bool flag22 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor2, d22 );
332 
333  if ( flag1 && flag2 )
334  {
335 
336  if ( flag12 && flag22 )
337  { //take the minimal value
338  d1 = getValue( d1, d12 );
339  d2 = getValue( d2, d22 );
340  if (std::abs(d1) < std::abs(d2))
341  d = d1;
342  else
343  d = d2;
344  coeff = 1.5;
345  }
346  else
347  { //like first-order accurate case
348  //take the minimal value
349  if (std::abs(d1) < std::abs(d2))
350  d = d1;
351  else
352  d = d2;
353  }
354 
355  }
356  else //if not flag1 AND flag2 both true
357  {
358  if (flag1)
359  {
360 
361  if ( flag12 )
362  {
363  d1 = getValue( d1, d12 );
364  coeff = 1.5;
365  }
366  d = d1;
367 
368  }
369 
370  if (flag2)
371  {
372 
373  if ( flag22 )
374  {
375  d2 = getValue( d2, d22 );
376  coeff = 1.5;
377  }
378  d = d2;
379  }
380  }
381 
382  l.push_back( CoeffValue( coeff, d ) );
383 
384  } //end if flag1 or flag2
385 
386  *it1 = c;
387  *it2 = c;
388  } //end for each dimension
389 
390  //computation of the new value
391  return this->compute(l);
392 }
393 
394 
395 //-----------------------------------------------------------------------------
396 template <typename TImage, typename TSet>
397 inline
398 typename DGtal::L2SecondOrderLocalDistance<TImage, TSet>::Value
399 DGtal::L2SecondOrderLocalDistance<TImage, TSet>::compute
400 (List& aList) const
401 {
402  ASSERT(aList.size() > 0);
403 
404  if ( aList.size() == 1 )
405  {
406  CoeffValue pair = aList.back();
407  if (pair.second >= 0)
408  return ( (pair.second + 1.0)/(pair.first) );
409  else
410  return ( (pair.second - 1.0)/(pair.first) );
411  }
412  else
413  { //resolution
414  double a = 0;
415  double b = 0;
416  double c = -1;
417 
418  for (typename List::iterator it = aList.begin();
419  it != aList.end(); ++it)
420  {
421  double coeff = it->first;
422  Value v = it->second;
423 
424  a += (coeff*coeff);
425  b -= 2 * coeff * static_cast<double>(v);
426  c += static_cast<double>(v*v);
427  }
428 
429  //discriminant
430  double disc = b*b - 4*a*c;
431  ASSERT(disc >= 0);
432 
433  if ( b < 0 )
434  return static_cast<Value>( ( -b + std::sqrt(disc) ) / (2*a) );
435  else
436  return static_cast<Value>( ( -b - std::sqrt(disc) ) / (2*a) );
437  }
438 }
439 
440 //-----------------------------------------------------------------------------
441 template <typename TImage, typename TSet>
442 inline
443 typename DGtal::L2SecondOrderLocalDistance<TImage, TSet>::Value
444 DGtal::L2SecondOrderLocalDistance<TImage, TSet>
445 ::getValue(const Value& aValue1, const Value& aValue2) const
446 {
447  return (2.0*aValue1 - aValue2/2.0);
448 }
449 
450 //-----------------------------------------------------------------------------
451 template <typename TImage, typename TSet>
452 inline
453 void
454 DGtal::L2SecondOrderLocalDistance<TImage, TSet>
455 ::selfDisplay ( std::ostream & out ) const
456 {
457  out << "L2";
458 }
459 
460 ///////////////////////////////////////////////////////////////////////////////
461 //-----------------------------------------------------------------------------
462 template <typename TImage, typename TSet>
463 inline
464 DGtal::LInfLocalDistance<TImage,TSet>::LInfLocalDistance
465  (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
466 {
467  ASSERT( myImgPtr );
468  ASSERT( mySetPtr );
469 }
470 
471 //-----------------------------------------------------------------------------
472 template <typename TImage, typename TSet>
473 inline
474 DGtal::LInfLocalDistance<TImage,TSet>::LInfLocalDistance
475  (const LInfLocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
476 {
477  ASSERT( myImgPtr );
478  ASSERT( mySetPtr );
479 }
480 
481 //-----------------------------------------------------------------------------
482 template <typename TImage, typename TSet>
483 inline
484 DGtal::LInfLocalDistance<TImage,TSet>&
485 DGtal::LInfLocalDistance<TImage,TSet>::operator=
486  (const LInfLocalDistance& other)
487 {
488  if( this != &other)
489  {
490  myImgPtr = other.myImgPtr;
491  mySetPtr = other.mySetPtr;
492  ASSERT( myImgPtr );
493  ASSERT( mySetPtr );
494  }
495  return *this;
496 }
497 
498 //-----------------------------------------------------------------------------
499 template <typename TImage, typename TSet>
500 inline
501 DGtal::LInfLocalDistance<TImage,TSet>::~LInfLocalDistance
502  ()
503 {
504 }
505 //-----------------------------------------------------------------------------
506 template <typename TImage, typename TSet>
507 inline
508 typename DGtal::LInfLocalDistance<TImage, TSet>::Value
509 DGtal::LInfLocalDistance<TImage, TSet>::operator()
510  (const Point& aPoint)
511 {
512 
513  //distance values
514  Values v;
515  v.reserve(Point::dimension);
516 
517  //two 1-neighbors
518  Point neighbor1 = aPoint;
519  Point neighbor2 = aPoint;
520 
521  typename Point::Iterator it1 = neighbor1.begin();
522  typename Point::Iterator it2 = neighbor2.begin();
523  typename Point::ConstIterator it = aPoint.begin();
524  typename Point::ConstIterator itEnd = aPoint.end();
525  for ( ; it != itEnd; ++it, ++it1, ++it2)
526  {//for each dimension
527 
528  typename Point::Coordinate c = *it;
529  *it1 = (c+1);
530  *it2 = (c-1);
531 
532  //neighboring values
533  Value d = 0, d1 = 0, d2 = 0;
534  bool flag1 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor1, d1 );
535  bool flag2 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor2, d2 );
536  if ( flag1 || flag2 )
537  {
538  if ( flag1 && flag2 )
539  { //take the minimal value
540  if (std::abs(d1) < std::abs(d2))
541  d = d1;
542  else
543  d = d2;
544  } else
545  {
546  if (flag1) d = d1;
547  if (flag2) d = d2;
548  }
549 
550  v.push_back(d);
551 
552  }
553 
554  *it1 = c;
555  *it2 = c;
556  } //end for each dimension
557 
558  //computation of the new value
559  return this->compute(v);
560 
561 }
562 
563 //-----------------------------------------------------------------------------
564 template <typename TImage, typename TSet>
565 inline
566 typename DGtal::LInfLocalDistance<TImage, TSet>::Value
567 DGtal::LInfLocalDistance<TImage, TSet>::compute
568 (Values& aValueList) const
569 {
570 
571  ASSERT(aValueList.size() > 0);
572 
573  if ( aValueList.size() == 1 )
574  {
575  Value d = aValueList.back();
576  if (d >= 0) ++d;
577  else --d;
578  return d;
579  }
580  else
581  { //max element
582  typename Values::iterator it =
583  std::max_element( aValueList.begin(), aValueList.end(), absComparator<Value> );
584  return *it;
585  }
586 }
587 
588 //-----------------------------------------------------------------------------
589 template <typename TImage, typename TSet>
590 inline
591 void
592 DGtal::LInfLocalDistance<TImage, TSet>
593 ::selfDisplay ( std::ostream & out ) const
594 {
595  out << "LInf";
596 }
597 
598 ///////////////////////////////////////////////////////////////////////////////
599 //-----------------------------------------------------------------------------
600 template <typename TImage, typename TSet>
601 inline
602 DGtal::L1LocalDistance<TImage,TSet>::L1LocalDistance
603  (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
604 {
605  ASSERT( myImgPtr );
606  ASSERT( mySetPtr );
607 }
608 
609 //-----------------------------------------------------------------------------
610 template <typename TImage, typename TSet>
611 inline
612 DGtal::L1LocalDistance<TImage,TSet>::L1LocalDistance
613  (const L1LocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
614 {
615  ASSERT( myImgPtr );
616  ASSERT( mySetPtr );
617 }
618 
619 //-----------------------------------------------------------------------------
620 template <typename TImage, typename TSet>
621 inline
622 DGtal::L1LocalDistance<TImage,TSet>&
623 DGtal::L1LocalDistance<TImage,TSet>::operator=
624  (const L1LocalDistance& other)
625 {
626  if( this != &other)
627  {
628  myImgPtr = other.myImgPtr;
629  mySetPtr = other.mySetPtr;
630  ASSERT( myImgPtr );
631  ASSERT( mySetPtr );
632  }
633  return *this;
634 }
635 
636 //-----------------------------------------------------------------------------
637 template <typename TImage, typename TSet>
638 inline
639 DGtal::L1LocalDistance<TImage,TSet>::~L1LocalDistance
640  ()
641 {
642 }
643 //-----------------------------------------------------------------------------
644 template <typename TImage, typename TSet>
645 inline
646 typename DGtal::L1LocalDistance<TImage, TSet>::Value
647 DGtal::L1LocalDistance<TImage, TSet>::operator()
648  (const Point& aPoint)
649 {
650 
651  //distance values
652  Values v;
653  v.reserve(2*Point::dimension);
654 
655  //two 1-neighbors
656  Point neighbor1 = aPoint;
657  Point neighbor2 = aPoint;
658 
659  typename Point::Iterator it1 = neighbor1.begin();
660  typename Point::Iterator it2 = neighbor2.begin();
661  typename Point::ConstIterator it = aPoint.begin();
662  typename Point::ConstIterator itEnd = aPoint.end();
663  for ( ; it != itEnd; ++it, ++it1, ++it2)
664  {//for each dimension
665 
666  typename Point::Coordinate c = *it;
667  *it1 = (c+1);
668  *it2 = (c-1);
669 
670  //neighboring values
671  Value d1 = 0, d2 = 0;
672  bool flag1 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor1, d1 );
673  bool flag2 = findAndGetValue( *myImgPtr, *mySetPtr, neighbor2, d2 );
674  if (flag1) v.push_back( d1 );
675  if (flag2) v.push_back( d2 );
676 
677  *it1 = c;
678  *it2 = c;
679  } //end for each dimension
680 
681  //computation of the new value
682  return this->compute(v);
683 
684 }
685 
686 //-----------------------------------------------------------------------------
687 template <typename TImage, typename TSet>
688 inline
689 typename DGtal::L1LocalDistance<TImage, TSet>::Value
690 DGtal::L1LocalDistance<TImage, TSet>::compute
691 (Values& aValueList) const
692 {
693  ASSERT(aValueList.size() > 0);
694 
695  //min (in absolute values)
696  typename Values::iterator it =
697  std::min_element( aValueList.begin(), aValueList.end(), absComparator<Value> );
698  Value vmin = *it;
699 
700  //sign
701  if (vmin >= 0)
702  return vmin + 1;
703  else
704  return vmin - 1;
705 }
706 
707 //-----------------------------------------------------------------------------
708 template <typename TImage, typename TSet>
709 inline
710 void
711 DGtal::L1LocalDistance<TImage, TSet>
712 ::selfDisplay ( std::ostream & out ) const
713 {
714  out << "L1";
715 }
716 // //
717 ///////////////////////////////////////////////////////////////////////////////
718 
719 ///////////////////////////////////////////////////////////////////////////////
720 ///////////////////////////////////////////////////////////////////////////////
721 // Helper classes defined in the compilation unit (anonymous namespace)
722 // see operator() below
723 namespace
724 {
725  //
726  template<bool complementToOne>
727  struct ValueBetween0And1
728  {
729  template<typename Value>
730  static Value get(const Value& v)
731  {
732  ASSERT( (v>=0)&&(v<=1) );
733  return v;
734  }
735  };
736  //specialization
737  template< >
738  struct ValueBetween0And1<true>
739  {
740  template<typename Value>
741  static Value get(const Value& v)
742  {
743  ASSERT( (v>=0)&&(v<=1) );
744  return (1 - v);
745  }
746  };
747 }
748 
749 ///////////////////////////////////////////////////////////////////////////////
750 //-----------------------------------------------------------------------------
751 template <typename TKSpace, typename TMap, bool isIndirect>
752 inline
753 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
754 ::L2FirstOrderLocalDistanceFromCells
755  ( ConstAlias<KSpace> aK, Map& aMap): myKSpace(&aK), myMap(&aMap)
756 {
757 }
758 
759 //-----------------------------------------------------------------------------
760 template <typename TKSpace, typename TMap, bool isIndirect>
761 inline
762 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
763 ::L2FirstOrderLocalDistanceFromCells
764  (const L2FirstOrderLocalDistanceFromCells& other)
765 : myKSpace(other.myKSpace), myMap(other.myMap)
766 {
767 }
768 
769 //-----------------------------------------------------------------------------
770 template <typename TKSpace, typename TMap, bool isIndirect>
771 inline
772 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>&
773 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
774 ::operator=
775  (const L2FirstOrderLocalDistanceFromCells& other)
776 {
777  if( this != &other)
778  {
779  myMap = other.myMap;
780  myKSpace = other.myKSpace;
781  }
782  return *this;
783 }
784 
785 //-----------------------------------------------------------------------------
786 template <typename TKSpace, typename TMap, bool isIndirect>
787 inline
788 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
789 ::~L2FirstOrderLocalDistanceFromCells
790  ()
791 {
792 }
793 
794 //-----------------------------------------------------------------------------
795 template <typename TKSpace, typename TMap, bool isIndirect>
796 inline
797 typename DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>::Value
798 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>::operator()
799  (const Point& aPoint)
800 {
801 
802  //distance values
803  Values v;
804  v.reserve(Point::dimension);
805 
806  Cell spel = myKSpace->uSpel( aPoint );
807  for ( typename KSpace::DirIterator q = myKSpace->uDirs( spel );
808  (q != 0); ++q )
809  { //for each dimension
810  const DGtal::Dimension dir = *q;
811 
812  /// for the direct orientation
813  Cell surfel1 = myKSpace->uIncident( spel, dir, true );
814  ASSERT( myKSpace->uDim( surfel1 ) == (KSpace::dimension - 1) );
815 
816  /// for the indirect orientation
817  Cell surfel2 = myKSpace->uIncident( spel, dir, false );
818  ASSERT( myKSpace->uDim( surfel2 ) == (KSpace::dimension - 1) );
819 
820  //neighboring values
821  Value d = 0;
822  typename Map::iterator it1 = myMap->find( surfel1 );
823  typename Map::iterator it2 = myMap->find( surfel2 );
824  bool flag1 = ( it1 != myMap->end() );
825  bool flag2 = ( it2 != myMap->end() );
826  if ( flag1 || flag2 )
827  {
828  if ( flag1 && flag2 )
829  { //take the minimal value
830  ASSERT( (it1->second >= 0)&&(it1->second <= 1) );
831  ASSERT( (it2->second >= 0)&&(it2->second <= 1) );
832  if (it1->second < it2->second)
833  d = ValueBetween0And1<isIndirect>
834  ::get( it1->second );
835  else
836  d = ValueBetween0And1<isIndirect>
837  ::get( it2->second );
838  } else
839  {
840  if (flag1)
841  {
842  ASSERT( (it1->second >= 0)&&(it1->second <= 1) );
843  d = ValueBetween0And1<isIndirect>
844  ::get( it1->second );
845  }
846  if (flag2)
847  {
848  ASSERT( (it2->second >= 0)&&(it2->second <= 1) );
849  d = ValueBetween0And1<isIndirect>
850  ::get( it2->second );
851  }
852  }
853 
854  if (d == 0) return 0;
855 
856  v.push_back(d);
857  }
858 
859  } //end for each dimension
860 
861  //computation of the new value
862  return this->compute(v);
863 }
864 
865 
866 //-----------------------------------------------------------------------------
867 template <typename TKSpace, typename TMap, bool isIndirect>
868 inline
869 typename DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>::Value
870 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>::compute
871 (Values& aValueList) const
872 {
873  ASSERT(aValueList.size() > 0);
874 
875  if ( aValueList.size() == 1 )
876  {
877  return aValueList.back();
878  }
879  else
880  { //resolution
881  double a = 0;
882 
883  for (typename Values::iterator it = aValueList.begin();
884  it != aValueList.end(); ++it)
885  {
886  Value d = *it;
887  a += ( 1.0 / static_cast<double>( d*d ) );
888  }
889 
890  return static_cast<Value>( std::sqrt(a) / a );
891  }
892 }
893 
894 //-----------------------------------------------------------------------------
895 template <typename TKSpace, typename TMap, bool isIndirect>
896 inline
897 void
898 DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap, isIndirect>
899 ::selfDisplay ( std::ostream & out ) const
900 {
901  out << "L2";
902 }
903 
904 ///////////////////////////////////////////////////////////////////////////////
905 //-----------------------------------------------------------------------------
906 template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
907 inline
908 DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::SpeedExtrapolator
909 (const DistanceImage& aDistImg, const Set& aSet, SpeedFunctor& aSpeedFunc)
910  : myDistImgPtr(&aDistImg), mySetPtr(&aSet), mySpeedFuncPtr(&aSpeedFunc)
911 {
912  ASSERT( myDistImgPtr );
913  ASSERT( mySetPtr );
914  ASSERT( mySpeedFuncPtr );
915 }
916 
917 //-----------------------------------------------------------------------------
918 template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
919 inline
920 DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::SpeedExtrapolator
921  (const SpeedExtrapolator& other)
922  : myDistImgPtr(other.myDistImgPtr), mySetPtr(other.mySetPtr), mySpeedFuncPtr(other.mySpeedFuncPtr)
923 {
924  ASSERT( myDistImgPtr );
925  ASSERT( mySetPtr );
926  ASSERT( mySpeedFuncPtr );
927 }
928 
929 //-----------------------------------------------------------------------------
930 template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
931 inline
932 DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>&
933 DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::operator=
934  (const SpeedExtrapolator& other)
935 {
936  if( this != &other)
937  {
938  myDistImgPtr = other.myDistImgPtr;
939  mySetPtr = other.mySetPtr;
940  mySpeedFuncPtr = other.mySpeedFuncPtr;
941  ASSERT( myDistImgPtr );
942  ASSERT( mySetPtr );
943  ASSERT( mySpeedFuncPtr );
944  }
945  return *this;
946 }
947 
948 //-----------------------------------------------------------------------------
949 template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
950 inline
951 DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::~SpeedExtrapolator
952  ()
953 {
954 }
955 
956 //-----------------------------------------------------------------------------
957 template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
958 inline
959 typename DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::Value
960 DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::operator()
961  (const Point& aPoint)
962 {
963 
964  //speed values
965  Value num = 0.0, den = 0;
966 
967  //two 1-neighbors for one dimension
968  Point neighbor1 = aPoint;
969  Point neighbor2 = aPoint;
970 
971  typename Point::Iterator it1 = neighbor1.begin();
972  typename Point::Iterator it2 = neighbor2.begin();
973  typename Point::ConstIterator it = aPoint.begin();
974  typename Point::ConstIterator itEnd = aPoint.end();
975  for ( ; it != itEnd; ++it, ++it1, ++it2)
976  {//for each dimension
977 
978  typename Point::Coordinate c = *it;
979 
980  //distance value
981  DistanceValue d = 0;
982  bool flag = findAndGetValue( *myDistImgPtr, *mySetPtr, aPoint, d );
983  ASSERT( flag ); boost::ignore_unused_variable_warning( flag );
984 
985  //neighbors
986  *it1 = (c+1);
987  *it2 = (c-1);
988  //neighboring speed value
989  Value s = 0;
990  //neighboring distance values
991  DistanceValue d0 = 0, d1 = 0, d2 = 0;
992  bool flag1 = findAndGetValue( *myDistImgPtr, *mySetPtr, neighbor1, d1 );
993  bool flag2 = findAndGetValue( *myDistImgPtr, *mySetPtr, neighbor2, d2 );
994  if ( flag1 || flag2 )
995  {
996  if ( flag1 && flag2 )
997  { //take the minimal value
998  if (std::abs(d1) < std::abs(d2))
999  {
1000  d0 = d1;
1001  s = (*mySpeedFuncPtr)( neighbor1 );
1002  }
1003  else
1004  {
1005  d0 = d2;
1006  s = (*mySpeedFuncPtr)( neighbor2 );
1007  }
1008  } else
1009  {
1010  if (flag1)
1011  {
1012  d0 = d1;
1013  s = (*mySpeedFuncPtr)( neighbor1 );
1014  }
1015  if (flag2)
1016  {
1017  d0 = d2;
1018  s = (*mySpeedFuncPtr)( neighbor2 );
1019  }
1020  }
1021 
1022  Value diff = static_cast<Value>(d - d0);
1023  num += ( s * diff );
1024  den += diff;
1025  }
1026 
1027  *it1 = c;
1028  *it2 = c;
1029  } //end for each dimension
1030 
1031  //computation of the new value
1032  ASSERT(den != 0);
1033  return (num/den);
1034 }
Generated on Tue Jun 1 2021 17:05:17 for DGtal by doxygen 1.9.1