DGtal 1.4.0
Loading...
Searching...
No Matches
FrechetShortcut.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 FrechetShortcut.ih
19 * @author Isabelle Sivignon (\c isabelle.sivignon@gipsa-lab.grenoble-inp.fr )
20 * gipsa-lab Grenoble Images Parole Signal Automatique (CNRS, UMR 5216), CNRS, France
21 *
22 * @date 2012/02/24
23 *
24 * Implementation of inline methods defined in FrechetShortcut.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29///////////////////////////////////////////////////////////////////////////////
30// IMPLEMENTATION of inline methods.
31///////////////////////////////////////////////////////////////////////////////
32
33//////////////////////////////////////////////////////////////////////////////
34#include <cstdlib>
35//////////////////////////////////////////////////////////////////////////////
36
37
38///////////////////////////////////////////////////////////////
39// Class backpath
40////////////////////////////////////////////////////////////////
41
42//#define PRECISION 0.00001
43
44
45//creation of a backPath
46// Default constructor
47template <typename TIterator, typename TInteger>
48inline
49DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath()
50{
51 myQuad = 0;
52 myFlag = false;
53}
54
55
56//creation of a backPath
57template <typename TIterator, typename TInteger>
58inline
59DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath(const FrechetShortcut<TIterator,TInteger> *s,int q): myS(s),myQuad(q),myFlag(false)
60{
61}
62
63template <typename TIterator, typename TInteger>
64inline
65DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath(const Backpath & other): myS(other.myS),myQuad(other.myQuad),myFlag(other.myFlag),myOcculters(other.myOcculters),myForbiddenIntervals(other.myForbiddenIntervals)
66{
67}
68
69
70template <typename TIterator, typename TInteger>
71inline
72typename DGtal::FrechetShortcut<TIterator,TInteger>::Backpath& DGtal::FrechetShortcut<TIterator,TInteger>::Backpath::operator=(const Backpath & other)
73{
74 myS = other.myS;
75 myQuad = other.myQuad;
76 myFlag = other.myFlag;
77 myOcculters = other.myOcculters;
78 myForbiddenIntervals = IntervalSet(other.myForbiddenIntervals);
79 return *this;
80}
81
82
83template <typename TIterator, typename TInteger>
84inline
85void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::reset()
86{
87 myFlag = false;
88 myOcculters.clear();
89 myForbiddenIntervals.clear();
90}
91
92
93//destruction of a backPath
94template <typename TIterator, typename TInteger>
95inline
96DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::~Backpath()
97{ }
98
99template <typename TIterator, typename TInteger>
100inline
101void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::updateBackPathFirstQuad(int d, const ConstIterator& it)
102{
103
104 myIt = it;
105
106 switch(d)
107 {
108 case 0:
109 case 1:
110 case 2:
111 case 7:
112 {
113 addPositivePoint();
114 break;
115 }
116 case 3:
117 case 4:
118 case 5:
119 case 6:
120 {
121 addNegativePoint();
122 break;
123 }
124 }
125}
126
127// update the list of active occulters
128template <typename TIterator, typename TInteger>
129inline
130void DGtal::FrechetShortcut<TIterator,TInteger>::Backpath::updateOcculters()
131{
132
133 // The potential new occulter is the last-but-one point
134 Point p = Point(*(myIt-1));
135
136
137 Point pi,v;
138 Vector dir = Tools::chainCode2Vect(myQuad);
139 Vector dir_ortho = Tools::chainCode2Vect((myQuad+6)%8);
140
141
142 Point u1,u2;
143 u1 = dir;
144 u2 = Tools::chainCode2Vect((myQuad+1)%8);;
145
146 double angle_min=0;
147 double angle_max=M_PI_4;
148 bool occ = false;
149 bool ok = true;
150
151 IntegerComputer<TInteger> ic;
152
153 if(myOcculters.size()==0)
154 {
155 occ =true;
156 angle_min=0;
157 angle_max=M_PI_4;
158 }
159 else
160 {
161 typename occulter_list::iterator iter=myOcculters.begin();
162
163 iter = myOcculters.begin();
164 for(typename occulter_list::size_type i=0; i < myOcculters.size() && ok ; ++i)
165 {
166 pi = Point(*(iter->first));
167 v = p-pi;
168
169 // pi is after p for all directions -> p is not an occulter
170 if(ic.dotProduct(v,u1) < 0 && ic.dotProduct(v,u2) <0)
171 {
172 ok = false;
173 occ = false;
174 }
175 else
176 // p is after pi for all directions -> pi is not an occulter
177 // anymore, p is a new occulter.
178 if(ic.dotProduct(v,u1) > 0 && ic.dotProduct(v,u2) > 0)
179 {
180 iter = myOcculters.erase(iter);
181 occ = true;
182 angle_min = 0;
183 angle_max = M_PI_4;
184 }
185 else
186 // p is after pi on [0,alpha], before pi on [alpha,pi/4]
187 if(ic.dotProduct(v,u1) > 0 && ic.dotProduct(v,u2) <= 0)
188 {
189 double alpha = Tools::angleVectVect(v,dir_ortho);
190
191 if(alpha >= iter->second.angle_min && alpha <=
192 iter->second.angle_max)
193 {
194 // p is a new occulter
195 occ = true;
196 angle_min = 0;
197 angle_max = alpha;
198 // pi's angle_min is updated
199 iter->second.angle_min = alpha;
200 iter++;
201 }
202 else
203 if(alpha > iter->second.angle_max)
204 {
205 //pi is not an occulter anymore
206 iter = myOcculters.erase(iter);
207 occ=true;
208 angle_min = 0;
209 angle_max = M_PI_4;
210 }
211 // if alpha < iter->second.angle_min, pi does not
212 // change, p may be an occulter -> do nothing
213 }
214 else // scalar_product(v,u1) < 0 && scalar_product(v,u2) > 0
215 // p is after pi on [alpha,pi/4], before pi on [0,alpha]
216 {
217 double alpha = Tools::angleVectVect(v,dir_ortho);
218 alpha = M_PI - alpha;
219
220 if(alpha >= iter->second.angle_min && alpha <=
221 iter->second.angle_max)
222 {
223 occ = true;
224 angle_min = alpha;
225 angle_max = M_PI_4;
226 // pi's angle_max is updated
227 iter->second.angle_max = alpha;
228 iter++;
229 }
230 else
231 if(alpha < iter->second.angle_min)
232 {
233 //pi is not an occulter anymore
234 iter = myOcculters.erase(iter);
235 occ=true;
236 angle_min = 0;
237 angle_max = M_PI_4;
238 }
239 else
240 iter++;
241 // if(alpha > iter->second.angle_max), pi does not
242 // change, p may be an occulter -> do nothing
243
244 }
245 }
246 }
247
248 if(occ)
249 {
250 occulter_attributes new_occ;
251 new_occ.angle_min = angle_min;
252 new_occ.angle_max = angle_max;
253 myOcculters.insert(myOcculters.end(),std::pair<const ConstIterator,occulter_attributes>(myIt-1,new_occ));
254
255 }
256
257
258}
259
260
261
262// update the set of intervals
263template <typename TIterator, typename TInteger>
264inline
265void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::updateIntervals()
266{
267 Point p = Point(*myIt);
268
269 Point pi,v;
270 Vector dir,dir1;
271
272 IntegerComputer<TInteger> ic;
273
274 for(typename occulter_list::iterator iter =
275 myOcculters.begin(); iter!=myOcculters.end() ;++iter)
276 {
277 pi = Point(*(iter->first));
278
279 v = p-pi;
280
281 dir = Tools::chainCode2Vect(myQuad);
282 dir1 = Tools::chainCode2Vect((myQuad+1)%8);
283
284 if(ic.dotProduct(v,dir)<0 || ic.dotProduct(v,dir1)<0)
285 {
286 if(v.norm()>=myS->myError/sqrt(2.0F))
287 {
288 if(ic.crossProduct(dir,v)<=0)
289 {
290 v[0] = -v[0];
291 v[1] = -v[1];
292 }
293 double angle_v = Tools::angleVectVect(v,dir);
294
295 double tmp = acos((double) myS->myError/(sqrt(2.0F)*v.norm()));
296 double angle1 = -tmp+angle_v;
297 double angle2 = tmp+angle_v;
298 if(angle1 < 0)
299 angle1 = 0;
300 if(angle2 > M_PI_4)
301 angle2 = M_PI_4;
302
303 // Define a new interval of forbidden angles and insert it in the list.
304 boost::icl::interval<double>::type s = boost::icl::interval<double>::closed(angle1,angle2);
305 myForbiddenIntervals.insert(s);
306
307 }
308 }
309 }
310
311}
312
313
314// update the length of the longest backpath on a curve part
315template <typename TIterator, typename TInteger>
316inline
317void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::addPositivePoint()
318{
319 // if we were on a monotone backpath, the point is an end of backpath
320 // otherwise, do nothing
321 if(myFlag)
322 {
323 myFlag=false;
324 }
325}
326
327
328
329/************************************************************/
330
331
332template <typename TIterator, typename TInteger>
333inline
334void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::addNegativePoint()
335{
336
337 // if we were on a monotone backpath, do nothing, the backpath
338 // continues
339 // otherwise it is the beggining of a new monotone backpath,
340 // possibly a locally maximal occulting point
341
342 //trace.info() << "add negative point" << std::endl;
343
344 if(!myFlag)
345 {
346 myFlag=true;
347 updateOcculters();
348 updateIntervals();
349 }
350 else
351 {
352 updateIntervals();
353 }
354
355
356}
357///////////////////////////////////////////////////////////////////////////
358// End of class backpath
359////////////////////////////////////////////////////////////////////////////
360
361//////////////////////////////////////////////////////////////////////////
362// Class cone
363////////////////////////////////////////////////////////////////////////////
364
365
366//creation of a cone
367
368template <typename TIterator, typename TInteger>
369inline
370DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone()
371{
372 myInf = true;
373 myMin = 0;
374 myMax = 2*M_PI;
375}
376
377
378
379template <typename TIterator, typename TInteger>
380inline
381DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone(double angle0, double angle1)
382{
383
384 // angle0 and angle1 are ordered in direct orientation such that the
385 // angle made by the two directions is lower than PI.
386
387
388 // case angle0-angle1 = PI -> infinite cone
389 if(fabs(fabs(angle0-angle1)-M_PI) < PRECISION)
390 {
391 // the orientation is supposed to be ok (depends on the points involved)
392 myMin = angle0;
393 myMax = angle1;
394 }
395 else
396 if(fabs(angle0-angle1)<M_PI)
397 {
398 if(angle0-angle1>0)
399 {
400 myMin = angle1;
401 myMax = angle0;
402 }
403 else
404 {
405 myMin = angle0;
406 myMax = angle1;
407 }
408 }
409 else
410 {
411 // the cone includes the direction of angle=0
412 if(angle0>angle1)
413 {
414 myMin = angle0;
415 myMax = angle1;
416 }
417 else
418 {
419 myMin = angle1;
420 myMax = angle0;
421 }
422 }
423 myInf = false;
424}
425
426template <typename TIterator, typename TInteger>
427inline
428DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone(double x, double y, double x0, double y0, double x1, double y1)
429{
430 double angle0 = Tools::computeAngle(x, y, x0, y0);
431 double angle1 = Tools::computeAngle(x, y, x1, y1);
432
433 assert(angle0 != -1 && angle1 != -1);
434
435 *this = Cone(angle0,angle1);
436 myInf = false;
437}
438
439template <typename TIterator, typename TInteger>
440inline
441bool DGtal::FrechetShortcut<TIterator,TInteger>::Cone::isEmpty() const
442{
443 if(myInf)
444 return false;
445 else
446 // Fix 05/2024 to enable error = 0: a cone may be defined by two values myMin=myMax --> check for empty cone by setting myMin=myMax= -1 instead
447 if(myMin==-1) // and then myMax = -1 too: way to represent the empty intersection of two cones.
448 return true;
449 else
450 return false;
451}
452
453template <typename TIterator, typename TInteger>
454inline
455typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone& DGtal::FrechetShortcut<TIterator,TInteger>::Cone::operator=(const Cone& c)
456{
457 myMin =c.myMin;
458 myMax=c.myMax;
459 myInf = c.myInf;
460 return *this;
461}
462
463// // Computes the symmetrical cone
464template <typename TIterator, typename TInteger>
465inline
466typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone DGtal::FrechetShortcut<TIterator,TInteger>::Cone::symmetricalCone()
467{
468 Cone cnew(myMin+M_PI,myMax+M_PI);
469 return cnew;
470}
471
472// Computes the intersection between the self Cone and another one.
473template <typename TIterator, typename TInteger>
474inline
475void DGtal::FrechetShortcut<TIterator,TInteger>::Cone::intersectCones(Cone c)
476{
477 Cone res;
478
479 // computes the intersection between the self cone and one half of another cone
480 res = intersectConesSimple(c);
481
482 // if they are disjoint, try the intersection with the other half
483 if(res.isEmpty())
484 {
485 Cone sym = c.symmetricalCone();
486 res = intersectConesSimple(sym);
487 }
488
489 *this = res;
490}
491
492
493//intersection of the self cone with another cone: considers only one half of the cone
494template <typename TIterator, typename TInteger>
495inline
496typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone DGtal::FrechetShortcut<TIterator,TInteger>::Cone::intersectConesSimple(Cone c)
497{
498 Cone res;
499
500 // if the cone is infinite, the new cone is c
501 if(myInf)
502 {
503 res = c;
504 res.myInf = false;
505 }
506 else
507 // the directions of the new cone are not included in the old one
508 if(!Tools::isBetween(c.myMin, myMin, myMax, 2*M_PI) && !Tools::isBetween(c.myMax, myMin,
509 myMax,
510 2*M_PI))
511 {
512 // first possibility: the cones are disjoint
513 if(!Tools::isBetween(myMin, c.myMin, c.myMax, 2*M_PI) && !Tools::isBetween(myMax, c.myMin,
514 c.myMax, 2*M_PI))
515 res = Cone(-1,-1); // empty cone: both angles are set to -1
516 else
517 // or the new cone includes the old one, nothing changes, the cone remains the same.
518 res = *this;
519 }
520 else
521 // the old cone is "cut" by the new one
522 if(Tools::isBetween(c.myMin, myMin, myMax, 2*M_PI))
523 if(Tools::isBetween(c.myMax, myMin, myMax, 2*M_PI))
524 res = c;
525 else
526 res = Cone(c.myMin, myMax);
527 else
528 res = Cone(myMin,c.myMax);
529
530
531 return res;
532}
533
534template <typename TIterator, typename TInteger>
535inline
536void DGtal::FrechetShortcut<TIterator,TInteger>::Cone::selfDisplay ( std::ostream & out) const
537{
538 out << "[Cone]" << std::endl;
539 if(myInf)
540 out << "Infinite" << std::endl;
541 else
542 out << "[Cone min = " << myMin << " max = " << myMax << "]" << std::endl;
543 out << "[End Cone]" << std::endl;
544}
545
546
547/////////////////////////////////////////////////////////////////////////////
548/// FrechetShortcut class
549/////////////////////////////////////////////////////////////////////////////
550
551///////////////////////////////////////////////////////////////////////////////
552// Implementation of inline methods //
553
554template <typename TIterator, typename TInteger>
555inline
556DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut()
557{
558 myError = 0;
559 myCone = Cone();
560
561 for(int i=0;i<8;i++)
562 {
563 //backpath b(i,0);
564 Backpath b(this,i);
565 myBackpath.push_back(b);
566 }
567
568}
569
570
571template <typename TIterator, typename TInteger>
572inline
573DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut(double error)
574{
575
576 myError = error;
577 myCone = Cone();
578
579 for(int i=0;i<8;i++)
580 {
581 //backpath b(i,error);
582 Backpath b(this,i);
583 myBackpath.push_back(b);
584 }
585}
586
587
588
589
590template <typename TIterator, typename TInteger>
591inline
592void DGtal::FrechetShortcut<TIterator,TInteger>::init(const ConstIterator& it)
593{
594 myBegin = it;
595 myEnd = it;
596 resetCone();
597 resetBackpath();
598}
599
600
601template <typename TIterator, typename TInteger>
602inline
603DGtal::FrechetShortcut<TIterator,TInteger> DGtal::FrechetShortcut<TIterator,TInteger>::getSelf()
604{
605 FrechetShortcut<TIterator,TInteger> other = FrechetShortcut(myError);
606 return other;
607}
608
609
610template <typename TIterator, typename TInteger>
611inline
612DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut (const FrechetShortcut<TIterator,TInteger> & other ) : myError(other.myError), myBackpath(other.myBackpath), myCone(other.myCone), myBegin(other.myBegin), myEnd(other.myEnd){
613 resetBackpath();
614 resetCone();
615
616}
617
618template <typename TIterator, typename TInteger>
619inline
620DGtal::FrechetShortcut<TIterator,TInteger> & DGtal::FrechetShortcut<TIterator,TInteger>::operator=(const FrechetShortcut<TIterator,TInteger> & other)
621{
622
623 if(this != &other)
624 {
625 myError = other.myError;
626 myBackpath = other.myBackpath;
627 myCone = other.myCone;
628 myBegin = other.myBegin;
629 myEnd = other.myEnd;
630 }
631 return *this;
632}
633
634template <typename TIterator, typename TInteger>
635inline
636DGtal::FrechetShortcut<std::reverse_iterator<TIterator>,TInteger>
637DGtal::FrechetShortcut<TIterator,TInteger>
638::getReverse() const
639{
640 return Reverse(myError);
641}
642
643template <typename TIterator, typename TInteger>
644inline
645bool
646DGtal::FrechetShortcut<TIterator,TInteger>::operator==(
647 const DGtal::FrechetShortcut<TIterator,TInteger>& other) const {
648 return ((myBegin == other.myBegin) && (myEnd == other.myEnd) && (myError == other.myError));
649}
650
651
652
653template <typename TIterator, typename TInteger>
654inline
655bool
656DGtal::FrechetShortcut<TIterator,TInteger>::operator!=(
657 const DGtal::FrechetShortcut<TIterator,TInteger>& other) const {
658 return (!(*this == other));
659}
660
661
662template <typename TIterator, typename TInteger>
663inline
664bool
665DGtal::FrechetShortcut<TIterator,TInteger>::extendFront()
666{
667 bool flag = (updateWidth() && updateBackpath());
668
669 if(flag)
670 ++myEnd;
671
672 return flag;
673}
674
675
676template <typename TIterator, typename TInteger>
677inline
678bool
679DGtal::FrechetShortcut<TIterator,TInteger>::isExtendableFront()
680{
681
682 return (testUpdateWidth() && testUpdateBackpath());
683
684}
685
686template <typename TIterator, typename TInteger>
687inline
688typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone
689DGtal::FrechetShortcut<TIterator,TInteger>::computeNewCone()
690{
691 double x0, y0,x1,y1;
692
693 Point firstP = Point(*myBegin);
694 Point newP = Point(*(myEnd+1));
695
696 Cone newCone=myCone;
697
698 if(firstP == newP)
699 return newCone;
700
701 // compute the tangent points defined by the first point and the
702 // circle C(newP,error)
703
704
705 bool intersect = Tools::circleTangentPoints(firstP[0],firstP[1], newP[0], newP[1], myError/(sqrt(2.0F)), &x0, &y0,
706 &x1, &y1);
707
708 if(intersect)
709 {
710 // define a cone according to the new tangent points
711 Cone c;
712 // case where there is one single tangent point
713 if(fabs(x0-x1) < PRECISION && fabs(y0-y1) < PRECISION)
714 {
715 double angle = Tools::computeAngle(firstP[0],firstP[1],newP[0],newP[1]);
716
717 // the cone is reduced to a line
718 c = Cone(angle,angle);
719 }
720 else
721 c = Cone(firstP[0],firstP[1],x0,y0,x1,y1);
722
723 newCone.intersectCones(c);
724 }
725
726
727 return newCone;
728
729}
730
731// Test if the new direction belongs to the new cone, but does not
732// modify myCone
733template <typename TIterator, typename TInteger>
734inline
735bool DGtal::FrechetShortcut<TIterator,TInteger>::testUpdateWidth()
736{
737 Cone c = computeNewCone();
738
739 Point firstP = Point(*myBegin);
740 Point newP = Point(*(myEnd+1));
741
742 if(!(c.isEmpty()))
743 if(c.myInf)
744 return true;
745 else
746 {
747 double angle = Tools::computeAngle(firstP[0], firstP[1], newP[0], newP[1]);
748 assert(angle != -1);
749 return Tools::isBetween(angle,c.myMin,c.myMax,2*M_PI);
750 }
751 else
752 return false;
753
754}
755
756
757template <typename TIterator, typename TInteger>
758inline
759bool DGtal::FrechetShortcut<TIterator,TInteger>::testUpdateBackpath()
760{
761 // Save the current value of the backpath
762 std::vector <typename DGtal::FrechetShortcut<TIterator,TInteger>::Backpath> BackpathSave;
763
764 for(unsigned int i=0;i<8;i++)
765 {
766 Backpath b(myBackpath[i]);
767 BackpathSave.push_back(b);
768 }
769
770 // Check whether the next point could be added or not with respect to the backpath
771 bool flag = updateBackpath();
772
773 // Copy back the values of backpath before the test.
774 for(unsigned int i=0;i<8;i++)
775 myBackpath[i] = Backpath(BackpathSave[i]);
776
777 return flag;
778
779}
780
781
782// Same as testUpdateWidth() but myCone is modified.
783template <typename TIterator, typename TInteger>
784inline
785bool DGtal::FrechetShortcut<TIterator,TInteger>::updateWidth()
786{
787 Cone c = computeNewCone();
788
789 myCone = c;
790
791 Point firstP = Point(*myBegin);
792 Point newP = Point(*(myEnd+1));
793
794 bool flag = true;
795
796 if(!(c.isEmpty()))
797 if(c.myInf)
798 flag = true;
799 else
800 {
801 double angle = Tools::computeAngle(firstP[0], firstP[1], newP[0],
802 newP[1]);
803 assert(angle != -1);
804 flag = Tools::isBetween(angle,c.myMin,c.myMax,2*M_PI);
805 }
806 else
807 flag = false;
808
809 return flag;
810}
811
812template <typename TIterator, typename TInteger>
813inline
814bool DGtal::FrechetShortcut<TIterator,TInteger>::updateBackpath()
815{
816 Point prevP = Point(*myEnd);
817 Point P = Point(*(myEnd+1));
818
819 int d = Tools::computeChainCode(prevP,P);
820
821 for(unsigned int j=0;j<8;j++)
822 myBackpath[j].updateBackPathFirstQuad(Tools::rot(d,j),myEnd+1);
823
824
825 return isBackpathOk();
826
827}
828
829template <typename TIterator, typename TInteger>
830inline
831bool DGtal::FrechetShortcut<TIterator,TInteger>::isBackpathOk()
832{
833 // compute the quadrant of the direction of P(i,j)
834
835 Point firstP = Point(*myBegin);
836 Point P = Point(*(myEnd+1));
837
838 int q = Tools::computeOctant(firstP,P);
839
840
841 // to handle non simple curves (a point is visited twice)
842 if(firstP==P)
843 return true;
844
845 // compute the direction vector pipj
846 Point v;
847 v[0] = P[0]-firstP[0];
848 v[1] = P[1]-firstP[1];
849
850 // compute the angle between the direction vector and the elementary
851 // direction (defined by the quadrant)
852 Point dir_elem = Tools::chainCode2Vect(q);
853
854 double angle = Tools::angleVectVect(v,dir_elem);
855
856 boost::icl::interval_set<double> intervals = myBackpath[q].myForbiddenIntervals;
857
858 if(boost::icl::contains(intervals,angle))
859 return false;
860
861 return true;
862
863}
864
865
866template <typename TIterator, typename TInteger>
867inline
868void DGtal::FrechetShortcut<TIterator,TInteger>::resetBackpath()
869{
870 for(unsigned int i=0;i<8;i++)
871 {
872 myBackpath[i].reset();
873 }
874}
875
876template <typename TIterator, typename TInteger>
877inline
878void DGtal::FrechetShortcut<TIterator,TInteger>::resetCone()
879{
880 myCone.myMin = 0;
881 myCone.myMax = 2*M_PI; // default cone is the whole space
882 myCone.myInf = true;
883}
884
885
886template <typename TIterator, typename TInteger>
887inline
888TIterator
889DGtal::FrechetShortcut<TIterator,TInteger>::begin() const {
890 return myBegin;
891}
892
893template <typename TIterator, typename TInteger>
894inline
895TIterator
896DGtal::FrechetShortcut<TIterator,TInteger>::end() const {
897 ConstIterator i(myEnd); ++i;
898 return i;
899}
900
901
902
903template <typename TIterator, typename TInteger>
904inline
905std::string
906DGtal::FrechetShortcut<TIterator,TInteger>::className() const
907{
908 return "FrechetShortcut";
909}
910
911
912
913/**
914 * Writes/Displays the object on an output stream.
915 * @param out the output stream where the object is written.
916 */
917
918template <typename TIterator, typename TInteger>
919inline
920void
921DGtal::FrechetShortcut<TIterator,TInteger>::selfDisplay ( std::ostream & out) const
922{
923
924 out << "[FrechetShortcut]" << std::endl;
925 out << "(Begin, End)=";
926 out << "("<< Point(*myBegin) << ", " << Point(*myEnd) << ")\n";
927 out << "[End FrechetShortcut]" << std::endl;
928
929}
930
931// Implementation of inline functions //
932
933template <typename TIterator, typename TInteger>
934inline
935std::ostream&
936DGtal::operator<< ( std::ostream & out,
937 const DGtal::FrechetShortcut<TIterator,TInteger> & object )
938{
939 object.selfDisplay( out );
940 return out;
941}
942
943
944
945
946// //
947///////////////////////////////////////////////////////////////////////////////
948
949