DGtal 1.3.0
Loading...
Searching...
No Matches
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
40namespace
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//-----------------------------------------------------------------------------
64template <typename TImage, typename TSet>
65inline
66DGtal::L2FirstOrderLocalDistance<TImage,TSet>::L2FirstOrderLocalDistance
67 (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
68{
69 ASSERT( myImgPtr );
70 ASSERT( mySetPtr );
71}
72
73//-----------------------------------------------------------------------------
74template <typename TImage, typename TSet>
75inline
76DGtal::L2FirstOrderLocalDistance<TImage,TSet>::L2FirstOrderLocalDistance
77 (const L2FirstOrderLocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
78{
79 ASSERT( myImgPtr );
80 ASSERT( mySetPtr );
81}
82
83//-----------------------------------------------------------------------------
84template <typename TImage, typename TSet>
85inline
86DGtal::L2FirstOrderLocalDistance<TImage,TSet>&
87DGtal::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//-----------------------------------------------------------------------------
101template <typename TImage, typename TSet>
102inline
103DGtal::L2FirstOrderLocalDistance<TImage,TSet>::~L2FirstOrderLocalDistance
104 ()
105{
106}
107
108//-----------------------------------------------------------------------------
109template <typename TImage, typename TSet>
110inline
111typename DGtal::L2FirstOrderLocalDistance<TImage, TSet>::Value
112DGtal::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//-----------------------------------------------------------------------------
166template <typename TImage, typename TSet>
167inline
168typename DGtal::L2FirstOrderLocalDistance<TImage, TSet>::Value
169DGtal::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//-----------------------------------------------------------------------------
220template <typename TImage, typename TSet>
221inline
222typename DGtal::L2FirstOrderLocalDistance<TImage, TSet>::Value
223DGtal::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//-----------------------------------------------------------------------------
237template <typename TImage, typename TSet>
238inline
239void
240DGtal::L2FirstOrderLocalDistance<TImage, TSet>
241::selfDisplay ( std::ostream & out ) const
242{
243 out << "L2";
244}
245
246///////////////////////////////////////////////////////////////////////////////
247//-----------------------------------------------------------------------------
248template <typename TImage, typename TSet>
249inline
250DGtal::L2SecondOrderLocalDistance<TImage,TSet>::L2SecondOrderLocalDistance
251 (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
252{
253 ASSERT( myImgPtr );
254 ASSERT( mySetPtr );
255}
256
257//-----------------------------------------------------------------------------
258template <typename TImage, typename TSet>
259inline
260DGtal::L2SecondOrderLocalDistance<TImage,TSet>::L2SecondOrderLocalDistance
261 (const L2SecondOrderLocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
262{
263 ASSERT( myImgPtr );
264 ASSERT( mySetPtr );
265}
266
267//-----------------------------------------------------------------------------
268template <typename TImage, typename TSet>
269inline
270DGtal::L2SecondOrderLocalDistance<TImage,TSet>&
271DGtal::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//-----------------------------------------------------------------------------
285template <typename TImage, typename TSet>
286inline
287DGtal::L2SecondOrderLocalDistance<TImage,TSet>::~L2SecondOrderLocalDistance
288 ()
289{
290}
291
292//-----------------------------------------------------------------------------
293template <typename TImage, typename TSet>
294inline
295typename DGtal::L2SecondOrderLocalDistance<TImage, TSet>::Value
296DGtal::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//-----------------------------------------------------------------------------
396template <typename TImage, typename TSet>
397inline
398typename DGtal::L2SecondOrderLocalDistance<TImage, TSet>::Value
399DGtal::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//-----------------------------------------------------------------------------
441template <typename TImage, typename TSet>
442inline
443typename DGtal::L2SecondOrderLocalDistance<TImage, TSet>::Value
444DGtal::L2SecondOrderLocalDistance<TImage, TSet>
445::getValue(const Value& aValue1, const Value& aValue2) const
446{
447 return (2.0*aValue1 - aValue2/2.0);
448}
449
450//-----------------------------------------------------------------------------
451template <typename TImage, typename TSet>
452inline
453void
454DGtal::L2SecondOrderLocalDistance<TImage, TSet>
455::selfDisplay ( std::ostream & out ) const
456{
457 out << "L2";
458}
459
460///////////////////////////////////////////////////////////////////////////////
461//-----------------------------------------------------------------------------
462template <typename TImage, typename TSet>
463inline
464DGtal::LInfLocalDistance<TImage,TSet>::LInfLocalDistance
465 (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
466{
467 ASSERT( myImgPtr );
468 ASSERT( mySetPtr );
469}
470
471//-----------------------------------------------------------------------------
472template <typename TImage, typename TSet>
473inline
474DGtal::LInfLocalDistance<TImage,TSet>::LInfLocalDistance
475 (const LInfLocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
476{
477 ASSERT( myImgPtr );
478 ASSERT( mySetPtr );
479}
480
481//-----------------------------------------------------------------------------
482template <typename TImage, typename TSet>
483inline
484DGtal::LInfLocalDistance<TImage,TSet>&
485DGtal::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//-----------------------------------------------------------------------------
499template <typename TImage, typename TSet>
500inline
501DGtal::LInfLocalDistance<TImage,TSet>::~LInfLocalDistance
502 ()
503{
504}
505//-----------------------------------------------------------------------------
506template <typename TImage, typename TSet>
507inline
508typename DGtal::LInfLocalDistance<TImage, TSet>::Value
509DGtal::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//-----------------------------------------------------------------------------
564template <typename TImage, typename TSet>
565inline
566typename DGtal::LInfLocalDistance<TImage, TSet>::Value
567DGtal::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//-----------------------------------------------------------------------------
589template <typename TImage, typename TSet>
590inline
591void
592DGtal::LInfLocalDistance<TImage, TSet>
593::selfDisplay ( std::ostream & out ) const
594{
595 out << "LInf";
596}
597
598///////////////////////////////////////////////////////////////////////////////
599//-----------------------------------------------------------------------------
600template <typename TImage, typename TSet>
601inline
602DGtal::L1LocalDistance<TImage,TSet>::L1LocalDistance
603 (Image& aImg, TSet& aSet): myImgPtr(&aImg), mySetPtr(&aSet)
604{
605 ASSERT( myImgPtr );
606 ASSERT( mySetPtr );
607}
608
609//-----------------------------------------------------------------------------
610template <typename TImage, typename TSet>
611inline
612DGtal::L1LocalDistance<TImage,TSet>::L1LocalDistance
613 (const L1LocalDistance& other): myImgPtr(other.myImgPtr), mySetPtr(other.mySetPtr)
614{
615 ASSERT( myImgPtr );
616 ASSERT( mySetPtr );
617}
618
619//-----------------------------------------------------------------------------
620template <typename TImage, typename TSet>
621inline
622DGtal::L1LocalDistance<TImage,TSet>&
623DGtal::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//-----------------------------------------------------------------------------
637template <typename TImage, typename TSet>
638inline
639DGtal::L1LocalDistance<TImage,TSet>::~L1LocalDistance
640 ()
641{
642}
643//-----------------------------------------------------------------------------
644template <typename TImage, typename TSet>
645inline
646typename DGtal::L1LocalDistance<TImage, TSet>::Value
647DGtal::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//-----------------------------------------------------------------------------
687template <typename TImage, typename TSet>
688inline
689typename DGtal::L1LocalDistance<TImage, TSet>::Value
690DGtal::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//-----------------------------------------------------------------------------
708template <typename TImage, typename TSet>
709inline
710void
711DGtal::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
723namespace
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//-----------------------------------------------------------------------------
751template <typename TKSpace, typename TMap, bool isIndirect>
752inline
753DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
754::L2FirstOrderLocalDistanceFromCells
755 ( ConstAlias<KSpace> aK, Map& aMap): myKSpace(&aK), myMap(&aMap)
756{
757}
758
759//-----------------------------------------------------------------------------
760template <typename TKSpace, typename TMap, bool isIndirect>
761inline
762DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
763::L2FirstOrderLocalDistanceFromCells
764 (const L2FirstOrderLocalDistanceFromCells& other)
765: myKSpace(other.myKSpace), myMap(other.myMap)
766{
767}
768
769//-----------------------------------------------------------------------------
770template <typename TKSpace, typename TMap, bool isIndirect>
771inline
772DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>&
773DGtal::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//-----------------------------------------------------------------------------
786template <typename TKSpace, typename TMap, bool isIndirect>
787inline
788DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>
789::~L2FirstOrderLocalDistanceFromCells
790 ()
791{
792}
793
794//-----------------------------------------------------------------------------
795template <typename TKSpace, typename TMap, bool isIndirect>
796inline
797typename DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>::Value
798DGtal::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//-----------------------------------------------------------------------------
867template <typename TKSpace, typename TMap, bool isIndirect>
868inline
869typename DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap,isIndirect>::Value
870DGtal::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//-----------------------------------------------------------------------------
895template <typename TKSpace, typename TMap, bool isIndirect>
896inline
897void
898DGtal::L2FirstOrderLocalDistanceFromCells<TKSpace,TMap, isIndirect>
899::selfDisplay ( std::ostream & out ) const
900{
901 out << "L2";
902}
903
904///////////////////////////////////////////////////////////////////////////////
905//-----------------------------------------------------------------------------
906template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
907inline
908DGtal::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//-----------------------------------------------------------------------------
918template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
919inline
920DGtal::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//-----------------------------------------------------------------------------
930template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
931inline
932DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>&
933DGtal::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//-----------------------------------------------------------------------------
949template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
950inline
951DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::~SpeedExtrapolator
952 ()
953{
954}
955
956//-----------------------------------------------------------------------------
957template <typename TDistanceImage, typename TSet, typename TSpeedFunctor>
958inline
959typename DGtal::SpeedExtrapolator<TDistanceImage,TSet,TSpeedFunctor>::Value
960DGtal::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 Sun Nov 27 2022 15:15:42 for DGtal by doxygen 1.9.5