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.
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.
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/>.
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
24 * Implementation of inline methods defined in FrechetShortcut.h
26 * This file is part of the DGtal library.
29 ///////////////////////////////////////////////////////////////////////////////
30 // IMPLEMENTATION of inline methods.
31 ///////////////////////////////////////////////////////////////////////////////
33 //////////////////////////////////////////////////////////////////////////////
35 //////////////////////////////////////////////////////////////////////////////
38 ///////////////////////////////////////////////////////////////
40 ////////////////////////////////////////////////////////////////
42 //#define PRECISION 0.00001
45 //creation of a backPath
46 // Default constructor
47 template <typename TIterator, typename TInteger>
49 DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath()
56 //creation of a backPath
57 template <typename TIterator, typename TInteger>
59 DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath(const FrechetShortcut<TIterator,TInteger> *s,int q): myS(s),myQuad(q),myFlag(false)
63 template <typename TIterator, typename TInteger>
65 DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath(const Backpath & other): myS(other.myS),myQuad(other.myQuad),myFlag(other.myFlag),myOcculters(other.myOcculters),myForbiddenIntervals(other.myForbiddenIntervals)
70 template <typename TIterator, typename TInteger>
72 typename DGtal::FrechetShortcut<TIterator,TInteger>::Backpath& DGtal::FrechetShortcut<TIterator,TInteger>::Backpath::operator=(const Backpath & other)
75 myQuad = other.myQuad;
76 myFlag = other.myFlag;
77 myOcculters = other.myOcculters;
78 myForbiddenIntervals = IntervalSet(other.myForbiddenIntervals);
83 template <typename TIterator, typename TInteger>
85 void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::reset()
89 myForbiddenIntervals.clear();
93 //destruction of a backPath
94 template <typename TIterator, typename TInteger>
96 DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::~Backpath()
99 template <typename TIterator, typename TInteger>
101 void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::updateBackPathFirstQuad(int d, const ConstIterator& it)
127 // update the list of active occulters
128 template <typename TIterator, typename TInteger>
130 void DGtal::FrechetShortcut<TIterator,TInteger>::Backpath::updateOcculters()
133 // The potential new occulter is the last-but-one point
134 Point p = Point(*(myIt-1));
138 Vector dir = Tools::chainCode2Vect(myQuad);
139 Vector dir_ortho = Tools::chainCode2Vect((myQuad+6)%8);
144 u2 = Tools::chainCode2Vect((myQuad+1)%8);;
147 double angle_max=M_PI_4;
151 IntegerComputer<TInteger> ic;
153 if(myOcculters.size()==0)
161 typename occulter_list::iterator iter=myOcculters.begin();
163 iter = myOcculters.begin();
164 for(typename occulter_list::size_type i=0; i < myOcculters.size() && ok ; ++i)
166 pi = Point(*(iter->first));
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)
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)
180 iter = myOcculters.erase(iter);
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)
189 double alpha = Tools::angleVectVect(v,dir_ortho);
191 if(alpha >= iter->second.angle_min && alpha <=
192 iter->second.angle_max)
194 // p is a new occulter
198 // pi's angle_min is updated
199 iter->second.angle_min = alpha;
203 if(alpha > iter->second.angle_max)
205 //pi is not an occulter anymore
206 iter = myOcculters.erase(iter);
211 // if alpha < iter->second.angle_min, pi does not
212 // change, p may be an occulter -> do nothing
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]
217 double alpha = Tools::angleVectVect(v,dir_ortho);
218 alpha = M_PI - alpha;
220 if(alpha >= iter->second.angle_min && alpha <=
221 iter->second.angle_max)
226 // pi's angle_max is updated
227 iter->second.angle_max = alpha;
231 if(alpha < iter->second.angle_min)
233 //pi is not an occulter anymore
234 iter = myOcculters.erase(iter);
241 // if(alpha > iter->second.angle_max), pi does not
242 // change, p may be an occulter -> do nothing
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));
262 // update the set of intervals
263 template <typename TIterator, typename TInteger>
265 void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::updateIntervals()
267 Point p = Point(*myIt);
272 IntegerComputer<TInteger> ic;
274 for(typename occulter_list::iterator iter =
275 myOcculters.begin(); iter!=myOcculters.end() ;++iter)
277 pi = Point(*(iter->first));
281 dir = Tools::chainCode2Vect(myQuad);
282 dir1 = Tools::chainCode2Vect((myQuad+1)%8);
284 if(ic.dotProduct(v,dir)<0 || ic.dotProduct(v,dir1)<0)
286 if(v.norm()>=myS->myError/sqrt(2.0F))
288 if(ic.crossProduct(dir,v)<=0)
293 double angle_v = Tools::angleVectVect(v,dir);
295 double tmp = acos((double) myS->myError/(sqrt(2.0F)*v.norm()));
296 double angle1 = -tmp+angle_v;
297 double angle2 = tmp+angle_v;
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);
314 // update the length of the longest backpath on a curve part
315 template <typename TIterator, typename TInteger>
317 void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::addPositivePoint()
319 // if we were on a monotone backpath, the point is an end of backpath
320 // otherwise, do nothing
329 /************************************************************/
332 template <typename TIterator, typename TInteger>
334 void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::addNegativePoint()
337 // if we were on a monotone backpath, do nothing, the backpath
339 // otherwise it is the beggining of a new monotone backpath,
340 // possibly a locally maximal occulting point
342 //trace.info() << "add negative point" << std::endl;
357 ///////////////////////////////////////////////////////////////////////////
358 // End of class backpath
359 ////////////////////////////////////////////////////////////////////////////
361 //////////////////////////////////////////////////////////////////////////
363 ////////////////////////////////////////////////////////////////////////////
368 template <typename TIterator, typename TInteger>
370 DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone()
379 template <typename TIterator, typename TInteger>
381 DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone(double angle0, double angle1)
384 // angle0 and angle1 are ordered in direct orientation such that the
385 // angle made by the two directions is lower than PI.
388 // case angle0-angle1 = PI -> infinite cone
389 if(fabs(fabs(angle0-angle1)-M_PI) < PRECISION)
391 // the orientation is supposed to be ok (depends on the points involved)
396 if(fabs(angle0-angle1)<M_PI)
411 // the cone includes the direction of angle=0
426 template <typename TIterator, typename TInteger>
428 DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone(double x, double y, double x0, double y0, double x1, double y1)
430 double angle0 = Tools::computeAngle(x, y, x0, y0);
431 double angle1 = Tools::computeAngle(x, y, x1, y1);
433 assert(angle0 != -1 && angle1 != -1);
435 *this = Cone(angle0,angle1);
439 template <typename TIterator, typename TInteger>
441 bool DGtal::FrechetShortcut<TIterator,TInteger>::Cone::isEmpty() const
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.
453 template <typename TIterator, typename TInteger>
455 typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone& DGtal::FrechetShortcut<TIterator,TInteger>::Cone::operator=(const Cone& c)
463 // // Computes the symmetrical cone
464 template <typename TIterator, typename TInteger>
466 typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone DGtal::FrechetShortcut<TIterator,TInteger>::Cone::symmetricalCone()
468 Cone cnew(myMin+M_PI,myMax+M_PI);
472 // Computes the intersection between the self Cone and another one.
473 template <typename TIterator, typename TInteger>
475 void DGtal::FrechetShortcut<TIterator,TInteger>::Cone::intersectCones(Cone c)
479 // computes the intersection between the self cone and one half of another cone
480 res = intersectConesSimple(c);
482 // if they are disjoint, try the intersection with the other half
485 Cone sym = c.symmetricalCone();
486 res = intersectConesSimple(sym);
493 //intersection of the self cone with another cone: considers only one half of the cone
494 template <typename TIterator, typename TInteger>
496 typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone DGtal::FrechetShortcut<TIterator,TInteger>::Cone::intersectConesSimple(Cone c)
500 // if the cone is infinite, the new cone is c
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,
512 // first possibility: the cones are disjoint
513 if(!Tools::isBetween(myMin, c.myMin, c.myMax, 2*M_PI) && !Tools::isBetween(myMax, c.myMin,
515 res = Cone(-1,-1); // empty cone: both angles are set to -1
517 // or the new cone includes the old one, nothing changes, the cone remains the same.
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))
526 res = Cone(c.myMin, myMax);
528 res = Cone(myMin,c.myMax);
534 template <typename TIterator, typename TInteger>
536 void DGtal::FrechetShortcut<TIterator,TInteger>::Cone::selfDisplay ( std::ostream & out) const
538 out << "[Cone]" << std::endl;
540 out << "Infinite" << std::endl;
542 out << "[Cone min = " << myMin << " max = " << myMax << "]" << std::endl;
543 out << "[End Cone]" << std::endl;
547 /////////////////////////////////////////////////////////////////////////////
548 /// FrechetShortcut class
549 /////////////////////////////////////////////////////////////////////////////
551 ///////////////////////////////////////////////////////////////////////////////
552 // Implementation of inline methods //
554 template <typename TIterator, typename TInteger>
556 DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut()
565 myBackpath.push_back(b);
571 template <typename TIterator, typename TInteger>
573 DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut(double error)
581 //backpath b(i,error);
583 myBackpath.push_back(b);
590 template <typename TIterator, typename TInteger>
592 void DGtal::FrechetShortcut<TIterator,TInteger>::init(const ConstIterator& it)
601 template <typename TIterator, typename TInteger>
603 DGtal::FrechetShortcut<TIterator,TInteger> DGtal::FrechetShortcut<TIterator,TInteger>::getSelf()
605 FrechetShortcut<TIterator,TInteger> other = FrechetShortcut(myError);
610 template <typename TIterator, typename TInteger>
612 DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut (const FrechetShortcut<TIterator,TInteger> & other ) : myError(other.myError), myBackpath(other.myBackpath), myCone(other.myCone), myBegin(other.myBegin), myEnd(other.myEnd){
618 template <typename TIterator, typename TInteger>
620 DGtal::FrechetShortcut<TIterator,TInteger> & DGtal::FrechetShortcut<TIterator,TInteger>::operator=(const FrechetShortcut<TIterator,TInteger> & other)
625 myError = other.myError;
626 myBackpath = other.myBackpath;
627 myCone = other.myCone;
628 myBegin = other.myBegin;
634 template <typename TIterator, typename TInteger>
636 DGtal::FrechetShortcut<std::reverse_iterator<TIterator>,TInteger>
637 DGtal::FrechetShortcut<TIterator,TInteger>
640 return Reverse(myError);
643 template <typename TIterator, typename TInteger>
646 DGtal::FrechetShortcut<TIterator,TInteger>::operator==(
647 const DGtal::FrechetShortcut<TIterator,TInteger>& other) const {
648 return ((myBegin == other.myBegin) && (myEnd == other.myEnd) && (myError == other.myError));
653 template <typename TIterator, typename TInteger>
656 DGtal::FrechetShortcut<TIterator,TInteger>::operator!=(
657 const DGtal::FrechetShortcut<TIterator,TInteger>& other) const {
658 return (!(*this == other));
662 template <typename TIterator, typename TInteger>
665 DGtal::FrechetShortcut<TIterator,TInteger>::extendFront()
667 bool flag = (updateWidth() && updateBackpath());
676 template <typename TIterator, typename TInteger>
679 DGtal::FrechetShortcut<TIterator,TInteger>::isExtendableFront()
682 return (testUpdateWidth() && testUpdateBackpath());
686 template <typename TIterator, typename TInteger>
688 typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone
689 DGtal::FrechetShortcut<TIterator,TInteger>::computeNewCone()
693 Point firstP = Point(*myBegin);
694 Point newP = Point(*(myEnd+1));
701 // compute the tangent points defined by the first point and the
702 // circle C(newP,error)
705 bool intersect = Tools::circleTangentPoints(firstP[0],firstP[1], newP[0], newP[1], myError/(sqrt(2.0F)), &x0, &y0,
710 // define a cone according to the new tangent points
712 // case where there is one single tangent point
713 if(fabs(x0-x1) < PRECISION && fabs(y0-y1) < PRECISION)
715 double angle = Tools::computeAngle(firstP[0],firstP[1],newP[0],newP[1]);
717 // the cone is reduced to a line
718 c = Cone(angle,angle);
721 c = Cone(firstP[0],firstP[1],x0,y0,x1,y1);
723 newCone.intersectCones(c);
731 // Test if the new direction belongs to the new cone, but does not
733 template <typename TIterator, typename TInteger>
735 bool DGtal::FrechetShortcut<TIterator,TInteger>::testUpdateWidth()
737 Cone c = computeNewCone();
739 Point firstP = Point(*myBegin);
740 Point newP = Point(*(myEnd+1));
747 double angle = Tools::computeAngle(firstP[0], firstP[1], newP[0], newP[1]);
749 return Tools::isBetween(angle,c.myMin,c.myMax,2*M_PI);
757 template <typename TIterator, typename TInteger>
759 bool DGtal::FrechetShortcut<TIterator,TInteger>::testUpdateBackpath()
761 // Save the current value of the backpath
762 std::vector <typename DGtal::FrechetShortcut<TIterator,TInteger>::Backpath> BackpathSave;
764 for(unsigned int i=0;i<8;i++)
766 Backpath b(myBackpath[i]);
767 BackpathSave.push_back(b);
770 // Check whether the next point could be added or not with respect to the backpath
771 bool flag = updateBackpath();
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]);
782 // Same as testUpdateWidth() but myCone is modified.
783 template <typename TIterator, typename TInteger>
785 bool DGtal::FrechetShortcut<TIterator,TInteger>::updateWidth()
787 Cone c = computeNewCone();
791 Point firstP = Point(*myBegin);
792 Point newP = Point(*(myEnd+1));
801 double angle = Tools::computeAngle(firstP[0], firstP[1], newP[0],
804 flag = Tools::isBetween(angle,c.myMin,c.myMax,2*M_PI);
812 template <typename TIterator, typename TInteger>
814 bool DGtal::FrechetShortcut<TIterator,TInteger>::updateBackpath()
816 Point prevP = Point(*myEnd);
817 Point P = Point(*(myEnd+1));
819 int d = Tools::computeChainCode(prevP,P);
821 for(unsigned int j=0;j<8;j++)
822 myBackpath[j].updateBackPathFirstQuad(Tools::rot(d,j),myEnd+1);
825 return isBackpathOk();
829 template <typename TIterator, typename TInteger>
831 bool DGtal::FrechetShortcut<TIterator,TInteger>::isBackpathOk()
833 // compute the quadrant of the direction of P(i,j)
835 Point firstP = Point(*myBegin);
836 Point P = Point(*(myEnd+1));
838 int q = Tools::computeOctant(firstP,P);
841 // to handle non simple curves (a point is visited twice)
845 // compute the direction vector pipj
847 v[0] = P[0]-firstP[0];
848 v[1] = P[1]-firstP[1];
850 // compute the angle between the direction vector and the elementary
851 // direction (defined by the quadrant)
852 Point dir_elem = Tools::chainCode2Vect(q);
854 double angle = Tools::angleVectVect(v,dir_elem);
856 boost::icl::interval_set<double> intervals = myBackpath[q].myForbiddenIntervals;
858 if(boost::icl::contains(intervals,angle))
866 template <typename TIterator, typename TInteger>
868 void DGtal::FrechetShortcut<TIterator,TInteger>::resetBackpath()
870 for(unsigned int i=0;i<8;i++)
872 myBackpath[i].reset();
876 template <typename TIterator, typename TInteger>
878 void DGtal::FrechetShortcut<TIterator,TInteger>::resetCone()
881 myCone.myMax = 2*M_PI; // default cone is the whole space
886 template <typename TIterator, typename TInteger>
889 DGtal::FrechetShortcut<TIterator,TInteger>::begin() const {
893 template <typename TIterator, typename TInteger>
896 DGtal::FrechetShortcut<TIterator,TInteger>::end() const {
897 ConstIterator i(myEnd); ++i;
903 template <typename TIterator, typename TInteger>
906 DGtal::FrechetShortcut<TIterator,TInteger>::className() const
908 return "FrechetShortcut";
914 * Writes/Displays the object on an output stream.
915 * @param out the output stream where the object is written.
918 template <typename TIterator, typename TInteger>
921 DGtal::FrechetShortcut<TIterator,TInteger>::selfDisplay ( std::ostream & out) const
924 out << "[FrechetShortcut]" << std::endl;
925 out << "(Begin, End)=";
926 out << "("<< Point(*myBegin) << ", " << Point(*myEnd) << ")\n";
927 out << "[End FrechetShortcut]" << std::endl;
931 // Implementation of inline functions //
933 template <typename TIterator, typename TInteger>
936 DGtal::operator<< ( std::ostream & out,
937 const DGtal::FrechetShortcut<TIterator,TInteger> & object )
939 object.selfDisplay( out );
947 ///////////////////////////////////////////////////////////////////////////////