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
47template <typename TIterator, typename TInteger>
49DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath()
56//creation of a backPath
57template <typename TIterator, typename TInteger>
59DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath(const FrechetShortcut<TIterator,TInteger> *s,int q): myS(s),myQuad(q),myFlag(false)
63template <typename TIterator, typename TInteger>
65DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::Backpath(const Backpath & other): myS(other.myS),myQuad(other.myQuad),myFlag(other.myFlag),myOcculters(other.myOcculters),myForbiddenIntervals(other.myForbiddenIntervals)
70template <typename TIterator, typename TInteger>
72typename 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);
83template <typename TIterator, typename TInteger>
85void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::reset()
89 myForbiddenIntervals.clear();
93//destruction of a backPath
94template <typename TIterator, typename TInteger>
96DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::~Backpath()
99template <typename TIterator, typename TInteger>
101void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::updateBackPathFirstQuad(int d, const ConstIterator& it)
127// update the list of active occulters
128template <typename TIterator, typename TInteger>
130void 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;
150 IntegerComputer<TInteger> ic;
152 if(myOcculters.size()==0)
160 typename occulter_list::iterator iter, next;
162 iter = myOcculters.begin();
163 while(iter!=myOcculters.end() && ok)
165 pi = Point(*(iter->first));
170 // pi is after p for all directions -> p is not an occulter
171 if(ic.dotProduct(v,u1) < 0 && ic.dotProduct(v,u2) <0)
177 // p is after pi for all directions -> pi is not an occulter
178 // anymore, p is a new occulter.
179 if(ic.dotProduct(v,u1) > 0 && ic.dotProduct(v,u2) > 0)
181 myOcculters.erase(iter);
187 // p is after pi on [0,alpha], before pi on [alpha,pi/4]
188 if(ic.dotProduct(v,u1) > 0 && ic.dotProduct(v,u2) <= 0)
190 double alpha = Tools::angleVectVect(v,dir_ortho);
192 if(alpha >= iter->second.angle_min && alpha <=
193 iter->second.angle_max)
195 // p is a new occulter
199 // pi's angle_min is updated
200 iter->second.angle_min = alpha;
203 if(alpha > iter->second.angle_max)
205 //pi is not an occulter anymore
206 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;
230 if(alpha < iter->second.angle_min)
232 //pi is not an occulter anymore
233 myOcculters.erase(iter);
238 // if(alpha > iter->second.angle_max), pi does not
239 // change, p may be an occulter -> do nothing
248 occulter_attributes new_occ;
249 new_occ.angle_min = angle_min;
250 new_occ.angle_max = angle_max;
251 myOcculters.insert(myOcculters.end(),std::pair<const ConstIterator,occulter_attributes>(myIt-1,new_occ));
260// update the set of intervals
261template <typename TIterator, typename TInteger>
263void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::updateIntervals()
265 Point p = Point(*myIt);
270 IntegerComputer<TInteger> ic;
272 for(typename occulter_list::iterator iter =
273 myOcculters.begin(); iter!=myOcculters.end() ;++iter)
275 pi = Point(*(iter->first));
279 dir = Tools::chainCode2Vect(myQuad);
280 dir1 = Tools::chainCode2Vect((myQuad+1)%8);
282 if(ic.dotProduct(v,dir)<0 || ic.dotProduct(v,dir1)<0)
284 if(v.norm()>=myS->myError/sqrt(2.0F))
286 if(ic.crossProduct(dir,v)<=0)
291 double angle_v = Tools::angleVectVect(v,dir);
293 double tmp = acos((double) myS->myError/(sqrt(2.0F)*v.norm()));
294 double angle1 = -tmp+angle_v;
295 double angle2 = tmp+angle_v;
301 // Define a new interval of forbidden angles and insert it in the list.
302 boost::icl::interval<double>::type s = boost::icl::interval<double>::closed(angle1,angle2);
303 myForbiddenIntervals.insert(s);
312// update the length of the longest backpath on a curve part
313template <typename TIterator, typename TInteger>
315void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::addPositivePoint()
317 // if we were on a monotone backpath, the point is an end of backpath
318 // otherwise, do nothing
327/************************************************************/
330template <typename TIterator, typename TInteger>
332void DGtal::FrechetShortcut<TIterator, TInteger>::Backpath::addNegativePoint()
335 // if we were on a monotone backpath, do nothing, the backpath
337 // otherwise it is the beggining of a new monotone backpath,
338 // possibly a locally maximal occulting point
340 //trace.info() << "add negative point" << std::endl;
355///////////////////////////////////////////////////////////////////////////
356// End of class backpath
357////////////////////////////////////////////////////////////////////////////
359//////////////////////////////////////////////////////////////////////////
361////////////////////////////////////////////////////////////////////////////
366template <typename TIterator, typename TInteger>
368DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone()
377template <typename TIterator, typename TInteger>
379DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone(double angle0, double angle1)
382 // angle0 and angle1 are ordered in direct orientation such that the
383 // angle made by the two directions is lower than PI.
386 // case angle0-angle1 = PI -> infinite cone
387 if(fabs(fabs(angle0-angle1)-M_PI) < PRECISION)
389 // the orientation is supposed to be ok (depends on the points involved)
394 if(fabs(angle0-angle1)<M_PI)
409 // the cone includes the direction of angle=0
424template <typename TIterator, typename TInteger>
426DGtal::FrechetShortcut<TIterator,TInteger>::Cone::Cone(double x, double y, double x0, double y0, double x1, double y1)
428 double angle0 = Tools::computeAngle(x, y, x0, y0);
429 double angle1 = Tools::computeAngle(x, y, x1, y1);
431 assert(angle0 != -1 && angle1 != -1);
433 *this = Cone(angle0,angle1);
437template <typename TIterator, typename TInteger>
439bool DGtal::FrechetShortcut<TIterator,TInteger>::Cone::isEmpty() const
450template <typename TIterator, typename TInteger>
452typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone& DGtal::FrechetShortcut<TIterator,TInteger>::Cone::operator=(const Cone& c)
460// // Computes the symmetrical cone
461template <typename TIterator, typename TInteger>
463typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone DGtal::FrechetShortcut<TIterator,TInteger>::Cone::symmetricalCone()
465 Cone cnew(myMin+M_PI,myMax+M_PI);
469// Computes the intersection between the self Cone and another one.
470template <typename TIterator, typename TInteger>
472void DGtal::FrechetShortcut<TIterator,TInteger>::Cone::intersectCones(Cone c)
476 // computes the intersection between the self cone and one half of another cone
477 res = intersectConesSimple(c);
479 // if they are disjoint, try the intersection with the other half
482 Cone sym = c.symmetricalCone();
483 res = intersectConesSimple(sym);
490//intersection of the self cone with another cone: considers only one half of the cone
491template <typename TIterator, typename TInteger>
493typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone DGtal::FrechetShortcut<TIterator,TInteger>::Cone::intersectConesSimple(Cone c)
497 // if the cone is infinite, the new cone is c
504 // the directions of the new cone are not included in the old one
505 if(!Tools::isBetween(c.myMin, myMin, myMax, 2*M_PI) && !Tools::isBetween(c.myMax, myMin,
509 // first possibility: the cones are disjoint
510 if(!Tools::isBetween(myMin, c.myMin, c.myMax, 2*M_PI) && !Tools::isBetween(myMax, c.myMin,
514 // or the new cone includes the old one, nothing changes, the cone remains the same.
518 // the old cone is "cut" by the new one
519 if(Tools::isBetween(c.myMin, myMin, myMax, 2*M_PI))
520 if(Tools::isBetween(c.myMax, myMin, myMax, 2*M_PI))
523 res = Cone(c.myMin, myMax);
525 res = Cone(myMin,c.myMax);
530template <typename TIterator, typename TInteger>
532void DGtal::FrechetShortcut<TIterator,TInteger>::Cone::selfDisplay ( std::ostream & out) const
534 out << "[Cone]" << std::endl;
536 out << "Infinite" << std::endl;
538 out << "[Cone min = " << myMin << " max = " << myMax << "]" << std::endl;
539 out << "[End Cone]" << std::endl;
543/////////////////////////////////////////////////////////////////////////////
544/// FrechetShortcut class
545/////////////////////////////////////////////////////////////////////////////
547///////////////////////////////////////////////////////////////////////////////
548// Implementation of inline methods //
550template <typename TIterator, typename TInteger>
552DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut()
561 myBackpath.push_back(b);
567template <typename TIterator, typename TInteger>
569DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut(double error)
576 //backpath b(i,error);
578 myBackpath.push_back(b);
585template <typename TIterator, typename TInteger>
587void DGtal::FrechetShortcut<TIterator,TInteger>::init(const ConstIterator& it)
596template <typename TIterator, typename TInteger>
598DGtal::FrechetShortcut<TIterator,TInteger> DGtal::FrechetShortcut<TIterator,TInteger>::getSelf()
600 FrechetShortcut<TIterator,TInteger> other = FrechetShortcut(myError);
605template <typename TIterator, typename TInteger>
607DGtal::FrechetShortcut<TIterator,TInteger>::FrechetShortcut (const FrechetShortcut<TIterator,TInteger> & other ) : myError(other.myError), myBackpath(other.myBackpath), myCone(other.myCone), myBegin(other.myBegin), myEnd(other.myEnd){
613template <typename TIterator, typename TInteger>
615DGtal::FrechetShortcut<TIterator,TInteger> & DGtal::FrechetShortcut<TIterator,TInteger>::operator=(const FrechetShortcut<TIterator,TInteger> & other)
620 myError = other.myError;
621 myBackpath = other.myBackpath;
622 myCone = other.myCone;
623 myBegin = other.myBegin;
629template <typename TIterator, typename TInteger>
631DGtal::FrechetShortcut<std::reverse_iterator<TIterator>,TInteger>
632DGtal::FrechetShortcut<TIterator,TInteger>
635 return Reverse(myError);
638template <typename TIterator, typename TInteger>
641DGtal::FrechetShortcut<TIterator,TInteger>::operator==(
642 const DGtal::FrechetShortcut<TIterator,TInteger>& other) const {
643 return ((myBegin == other.myBegin) && (myEnd == other.myEnd) && (myError == other.myError));
648template <typename TIterator, typename TInteger>
651DGtal::FrechetShortcut<TIterator,TInteger>::operator!=(
652 const DGtal::FrechetShortcut<TIterator,TInteger>& other) const {
653 return (!(*this == other));
657template <typename TIterator, typename TInteger>
660DGtal::FrechetShortcut<TIterator,TInteger>::extendFront()
662 bool flag = (updateWidth() && updateBackpath());
671template <typename TIterator, typename TInteger>
674DGtal::FrechetShortcut<TIterator,TInteger>::isExtendableFront()
677 return (testUpdateWidth() && testUpdateBackpath());
681template <typename TIterator, typename TInteger>
683typename DGtal::FrechetShortcut<TIterator,TInteger>::Cone
684DGtal::FrechetShortcut<TIterator,TInteger>::computeNewCone()
688 Point firstP = Point(*myBegin);
689 Point newP = Point(*(myEnd+1));
694 // compute the tangent points defined by the first point and the
695 // circle C(newP,error)
696 bool intersect = Tools::circleTangentPoints(firstP[0],firstP[1], newP[0], newP[1], myError/(sqrt(2.0F)), &x0, &y0,
701 // define a cone according to the new tangent points
703 // case where there is one single tangent point
704 if(fabs(x0-x1) < PRECISION && fabs(y0-y1) < PRECISION)
706 double angle = Tools::computeAngle(firstP[0],firstP[1],newP[0],newP[1]);
708 double angle0 = angle - M_PI_2;
710 angle0 = angle0+2*M_PI;
711 double angle1 = angle + M_PI_2;
713 angle1 = angle1-2*M_PI;
714 c = Cone(angle0,angle1);
717 c = Cone(firstP[0],firstP[1],x0,y0,x1,y1);
719 newCone.intersectCones(c);
726// Test if the new direction belongs to the new cone, but does not
728template <typename TIterator, typename TInteger>
730bool DGtal::FrechetShortcut<TIterator,TInteger>::testUpdateWidth()
732 Cone c = computeNewCone();
734 Point firstP = Point(*myBegin);
735 Point newP = Point(*(myEnd+1));
742 double angle = Tools::computeAngle(firstP[0], firstP[1], newP[0], newP[1]);
744 return Tools::isBetween(angle,c.myMin,c.myMax,2*M_PI);
752template <typename TIterator, typename TInteger>
754bool DGtal::FrechetShortcut<TIterator,TInteger>::testUpdateBackpath()
756 // Save the current value of the backpath
757 std::vector <typename DGtal::FrechetShortcut<TIterator,TInteger>::Backpath> BackpathSave;
759 for(unsigned int i=0;i<8;i++)
761 Backpath b(myBackpath[i]);
762 BackpathSave.push_back(b);
765 // Check whether the next point could be added or not with respect to the backpath
766 bool flag = updateBackpath();
768 // Copy back the values of backpath before the test.
769 for(unsigned int i=0;i<8;i++)
770 myBackpath[i] = Backpath(BackpathSave[i]);
777// Same as testUpdateWidth() but myCone is modified.
778template <typename TIterator, typename TInteger>
780bool DGtal::FrechetShortcut<TIterator,TInteger>::updateWidth()
782 Cone c = computeNewCone();
786 Point firstP = Point(*myBegin);
787 Point newP = Point(*(myEnd+1));
796 double angle = Tools::computeAngle(firstP[0], firstP[1], newP[0],
799 flag = Tools::isBetween(angle,c.myMin,c.myMax,2*M_PI);
807template <typename TIterator, typename TInteger>
809bool DGtal::FrechetShortcut<TIterator,TInteger>::updateBackpath()
811 Point prevP = Point(*myEnd);
812 Point P = Point(*(myEnd+1));
814 int d = Tools::computeChainCode(prevP,P);
816 for(unsigned int j=0;j<8;j++)
817 myBackpath[j].updateBackPathFirstQuad(Tools::rot(d,j),myEnd+1);
820 return isBackpathOk();
824template <typename TIterator, typename TInteger>
826bool DGtal::FrechetShortcut<TIterator,TInteger>::isBackpathOk()
828 // compute the quadrant of the direction of P(i,j)
830 Point firstP = Point(*myBegin);
831 Point P = Point(*(myEnd+1));
833 int q = Tools::computeOctant(firstP,P);
836 // to handle non simple curves (a point is visited twice)
840 // compute the direction vector pipj
842 v[0] = P[0]-firstP[0];
843 v[1] = P[1]-firstP[1];
845 // compute the angle between the direction vector and the elementary
846 // direction (defined by the quadrant)
847 Point dir_elem = Tools::chainCode2Vect(q);
849 double angle = Tools::angleVectVect(v,dir_elem);
851 boost::icl::interval_set<double> intervals = myBackpath[q].myForbiddenIntervals;
853 if(boost::icl::contains(intervals,angle))
861template <typename TIterator, typename TInteger>
863void DGtal::FrechetShortcut<TIterator,TInteger>::resetBackpath()
865 for(unsigned int i=0;i<8;i++)
867 myBackpath[i].reset();
871template <typename TIterator, typename TInteger>
873void DGtal::FrechetShortcut<TIterator,TInteger>::resetCone()
881template <typename TIterator, typename TInteger>
884DGtal::FrechetShortcut<TIterator,TInteger>::begin() const {
888template <typename TIterator, typename TInteger>
891DGtal::FrechetShortcut<TIterator,TInteger>::end() const {
892 ConstIterator i(myEnd); ++i;
898template <typename TIterator, typename TInteger>
901DGtal::FrechetShortcut<TIterator,TInteger>::className() const
903 return "FrechetShortcut";
909 * Writes/Displays the object on an output stream.
910 * @param out the output stream where the object is written.
913template <typename TIterator, typename TInteger>
916DGtal::FrechetShortcut<TIterator,TInteger>::selfDisplay ( std::ostream & out) const
919 out << "[FrechetShortcut]" << std::endl;
920 out << "(Begin, End)=";
921 out << "("<< Point(*myBegin) << ", " << Point(*myEnd) << ")\n";
922 out << "[End FrechetShortcut]" << std::endl;
926// Implementation of inline functions //
928template <typename TIterator, typename TInteger>
931DGtal::operator<< ( std::ostream & out,
932 const DGtal::FrechetShortcut<TIterator,TInteger> & object )
934 object.selfDisplay( out );
942///////////////////////////////////////////////////////////////////////////////