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 ArithmeticalDSS.ih
19 * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
20 * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
24 * Implementation of inline methods defined in ArithmeticalDSS.h
26 * This file is part of the DGtal library.
30//////////////////////////////////////////////////////////////////////////////
32//////////////////////////////////////////////////////////////////////////////
34///////////////////////////////////////////////////////////////////////////////
35// IMPLEMENTATION of inline methods.
36///////////////////////////////////////////////////////////////////////////////
38///////////////////////////////////////////////////////////////////////////////
39// ----------------------- Standard services ------------------------------
41template <typename TCoordinate, typename TInteger, unsigned short adjacency>
43DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
44::ArithmeticalDSS(const Point& aPoint)
46myF(aPoint), myL(aPoint),
47 myUf(aPoint), myUl(aPoint),
48 myLf(aPoint), myLl(aPoint),
49 myDSL( NumberTraits<TCoordinate>::ZERO,
50 NumberTraits<TCoordinate>::ZERO,
51 NumberTraits<TInteger>::ZERO,
52 NumberTraits<TInteger>::ZERO,
53 std::make_pair( Vector(NumberTraits<TCoordinate>::ZERO, NumberTraits<TCoordinate>::ZERO),
54 Vector(NumberTraits<TCoordinate>::ZERO, NumberTraits<TCoordinate>::ZERO) ),
55 Vector(NumberTraits<TCoordinate>::ZERO, NumberTraits<TCoordinate>::ZERO) )
59//-----------------------------------------------------------------------------
60template <typename TCoordinate, typename TInteger, unsigned short adjacency>
62DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
63::ArithmeticalDSS(const Coordinate& aA, const Coordinate& aB,
64 const Integer& aLowerBound, const Integer& aUpperBound,
65 const Point& aF, const Point& aL,
66 const Point& aUf, const Point& aUl,
67 const Point& aLf, const Point& aLl,
68 const Steps& aSteps, const Vector& aShift)
71 myUf(aUf), myUl(aUl), myLf(aLf), myLl(aLl),
72 myDSL( aA, aB, aLowerBound, aUpperBound, aSteps, aShift )
76//-----------------------------------------------------------------------------
77template <typename TCoordinate, typename TInteger, unsigned short adjacency>
79DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
80::ArithmeticalDSS(const DSL& aDSL,
81 const Point& aF, const Point& aL,
82 const Point& aUf, const Point& aUl,
83 const Point& aLf, const Point& aLl)
86 myUf(aUf), myUl(aUl), myLf(aLf), myLl(aLl),
91//-----------------------------------------------------------------------------
92template <typename TCoordinate, typename TInteger, unsigned short adjacency>
94DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
95::ArithmeticalDSS(const Coordinate& aA, const Coordinate& aB,
96 const Point& aF, const Point& aL,
97 const Point& aUf, const Point& aUl,
98 const Point& aLf, const Point& aLl)
101 myUf(aUf), myUl(aUl), myLf(aLf), myLl(aLl),
102 myDSL( aA, aB, DSL::remainder( aA, aB, aUf ) )
107//-----------------------------------------------------------------------------
108template <typename TCoordinate, typename TInteger, unsigned short adjacency>
110DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
111::ArithmeticalDSS(const Point& aF, const Point& aL,
112 const bool& areOnTheUpperLine)
113 : myDSL( NumberTraits<Coordinate>::ZERO, NumberTraits<Coordinate>::ZERO, NumberTraits<Integer>::ZERO )
115 typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
117 if (areOnTheUpperLine)
118 *this = Factory::createPattern(aF, aL);
120 *this = Factory::createReversedPattern(aF, aL);
123//-----------------------------------------------------------------------------
124template <typename TCoordinate, typename TInteger, unsigned short adjacency>
126DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
127::ArithmeticalDSS(const DSL& aDSL, const Point& aF, const Point& aL)
128 : myDSL( NumberTraits<Coordinate>::ZERO, NumberTraits<Coordinate>::ZERO, NumberTraits<Integer>::ZERO )
130 typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
131 *this = Factory::createSubsegment(aDSL, aF, aL);
134//-----------------------------------------------------------------------------
135template <typename TCoordinate, typename TInteger, unsigned short adjacency>
137DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
138::ArithmeticalDSS(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aDSS,
139 const Point& aF, const Point& aL)
140 : myDSL( NumberTraits<Coordinate>::ZERO, NumberTraits<Coordinate>::ZERO, NumberTraits<Integer>::ZERO )
142 typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
143 *this = Factory::createSubsegment(aDSS, aF, aL);
146//-----------------------------------------------------------------------------
147template <typename TCoordinate, typename TInteger, unsigned short adjacency>
148template <typename Iterator>
150DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
151::ArithmeticalDSS(const Iterator& aItb, const Iterator& aIte)
152 : myDSL( NumberTraits<Coordinate>::ZERO, NumberTraits<Coordinate>::ZERO, NumberTraits<Integer>::ZERO )
157 //construction from the begin iterator
159 *this = ArithmeticalDSS(p);
161 for (++it; it != aIte; ++it)
166 throw InputException();
170//-----------------------------------------------------------------------------
171template <typename TCoordinate, typename TInteger, unsigned short adjacency>
173DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
174::ArithmeticalDSS(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther)
176 myF(aOther.myF), myL(aOther.myL),
177 myUf(aOther.myUf), myUl(aOther.myUl), myLf(aOther.myLf), myLl(aOther.myLl),
181//-----------------------------------------------------------------------------
182template <typename TCoordinate, typename TInteger, unsigned short adjacency>
184DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>&
185DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
186::operator=(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther)
188 if ( this != &aOther )
196 myDSL = aOther.myDSL;
201//-----------------------------------------------------------------------------
202template <typename TCoordinate, typename TInteger, unsigned short adjacency>
204DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
205DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
208 return ArithmeticalDSS(-myDSL.myA, -myDSL.myB,
209 -myDSL.myUpperBound, -myDSL.myLowerBound,
210 myL, myF, myLl, myLf, myUl, myUf,
211 std::make_pair(-myDSL.mySteps.first, -myDSL.mySteps.second),
215//-----------------------------------------------------------------------------
216template <typename TCoordinate, typename TInteger, unsigned short adjacency>
219DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
220::equalsTo(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
222 return ( (myUf == aOther.myUf) &&
223 (myUl == aOther.myUl) &&
224 (myLf == aOther.myLf) &&
225 (myLl == aOther.myLl) &&
226 (myF == aOther.myF) &&
227 (myL == aOther.myL) &&
228 (myDSL.equalsTo(aOther.myDSL)) );
231//-----------------------------------------------------------------------------
232template <typename TCoordinate, typename TInteger, unsigned short adjacency>
235DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
236::operator==(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
238 return ( equalsTo(aOther) || equalsTo(aOther.negate()) );
241//-----------------------------------------------------------------------------
242template <typename TCoordinate, typename TInteger, unsigned short adjacency>
245DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
246::operator!=(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
248 return !( operator==(aOther) );
251//-----------------------------------------------------------------------------
252template <typename TCoordinate, typename TInteger, unsigned short adjacency>
254DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::~ArithmeticalDSS()
258//-----------------------------------------------------------------------------
259template <typename TCoordinate, typename TInteger, unsigned short adjacency>
262DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::isValid() const
264 return functions::checkAll(*this);
267//-----------------------------------------------------------------------------
268template <typename TCoordinate, typename TInteger, unsigned short adjacency>
270const typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::DSL&
271DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::dsl() const
276//-----------------------------------------------------------------------------
277template <typename TCoordinate, typename TInteger, unsigned short adjacency>
279typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Coordinate
280DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::a() const
285//-----------------------------------------------------------------------------
286template <typename TCoordinate, typename TInteger, unsigned short adjacency>
288typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Coordinate
289DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::b() const
294//-----------------------------------------------------------------------------
295template <typename TCoordinate, typename TInteger, unsigned short adjacency>
297typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
298DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::mu() const
303//-----------------------------------------------------------------------------
304template <typename TCoordinate, typename TInteger, unsigned short adjacency>
306typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
307DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::omega() const
309 return myDSL.omega();
312//-----------------------------------------------------------------------------
313template <typename TCoordinate, typename TInteger, unsigned short adjacency>
315typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Vector
316DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::shift() const
318 return myDSL.shift();
321//-----------------------------------------------------------------------------
322template <typename TCoordinate, typename TInteger, unsigned short adjacency>
324typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Steps
325DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::steps() const
327 return myDSL.steps();
330//-----------------------------------------------------------------------------
331template <typename TCoordinate, typename TInteger, unsigned short adjacency>
333typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
334DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::back() const
339//-----------------------------------------------------------------------------
340template <typename TCoordinate, typename TInteger, unsigned short adjacency>
342typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
343DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::front() const
348//-----------------------------------------------------------------------------
349template <typename TCoordinate, typename TInteger, unsigned short adjacency>
351typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
352DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Uf() const
357//-----------------------------------------------------------------------------
358template <typename TCoordinate, typename TInteger, unsigned short adjacency>
360typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
361DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Ul() const
366//-----------------------------------------------------------------------------
367template <typename TCoordinate, typename TInteger, unsigned short adjacency>
369typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
370DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Lf() const
375//-----------------------------------------------------------------------------
376template <typename TCoordinate, typename TInteger, unsigned short adjacency>
378typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
379DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Ll() const
384//-----------------------------------------------------------------------------
385template <typename TCoordinate, typename TInteger, unsigned short adjacency>
387typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
388DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::remainder(const Point& aPoint) const
390 return myDSL.remainder(aPoint);
393//-----------------------------------------------------------------------------
394template <typename TCoordinate, typename TInteger, unsigned short adjacency>
396typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
397DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::orthogonalPosition(const Point& aPoint) const
399 return myDSL.orthogonalPosition(aPoint);
402//-----------------------------------------------------------------------------
403template <typename TCoordinate, typename TInteger, unsigned short adjacency>
405typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Position
406DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::position(const Point& aPoint) const
408 return myDSL.position(aPoint);
411//-----------------------------------------------------------------------------
412template <typename TCoordinate, typename TInteger, unsigned short adjacency>
415DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::before(const Point& aP1, const Point& aP2) const
417 return myDSL.before(aP1, aP2);
420//-----------------------------------------------------------------------------
421template <typename TCoordinate, typename TInteger, unsigned short adjacency>
424DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::beforeOrEqual(const Point& aP1, const Point& aP2) const
426 return myDSL.beforeOrEqual(aP1, aP2);
429//-----------------------------------------------------------------------------
430template <typename TCoordinate, typename TInteger, unsigned short adjacency>
433DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::isInDSL(const Point& aPoint) const
435 return myDSL(aPoint);
438//-----------------------------------------------------------------------------
439template <typename TCoordinate, typename TInteger, unsigned short adjacency>
442DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::isInDSS(const Point& aPoint) const
446 Integer s = position(aPoint);
447 Integer s1 = position(myF);
448 Integer s2 = position(myL);
451 return (aPoint == myF);
453 return ( (s >= s1)&&(s <= s2) );
455 return ( (s >= s2)&&(s <= s1) );
462//----------------------------------------------------------------------------
463template <typename TCoordinate, typename TInteger, unsigned short adjacency>
466DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>::isInDSL(const DSL& aDSL) const
469 // Test whether the DSL and the DSS belong to the same octant or not
470 typename ArithmeticalDSL<TCoordinate,TInteger,adjacency>::Octant::first_type oc;
471 std::vector<Point> LPoints;
472 if(!(myDSL.sameOctant(aDSL, &oc)))
475 // Consider the leaning points of 'this'.
476 LPoints.push_back(myUf);
477 LPoints.push_back(myLf);
478 LPoints.push_back(myUl);
479 LPoints.push_back(myLl);
481 typename std::vector<Point>::const_iterator it = LPoints.begin();
484 while(it != LPoints.end() && inDSL)
497//----------------------------------------------------------------------------
498template <typename TCoordinate, typename TInteger, unsigned short adjacency>
501DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>::isInDSL(const DSL& aDSL, std::vector<Point> &Ulp, std::vector<Point> &Llp, Point &outP) const
503 // Test whether the DSL and the DSS belong to the same octant or not
504 typename ArithmeticalDSL<TCoordinate,TInteger,adjacency>::Octant::first_type oc;
505 std::vector<Point> LPoints;
506 Ulp.clear(); Llp.clear();
508 if(!(myDSL.sameOctant(aDSL, &oc)))
511 // Consider the leaning points of 'this'
512 LPoints.push_back(myLf);
513 LPoints.push_back(myUf);
514 LPoints.push_back(myLl);
515 LPoints.push_back(myUl);
517 typename std::vector<Point>::const_iterator it = LPoints.begin();
520 while(it != LPoints.end() && inDSL)
527 // store the points of the DSS that are leaning points for the DSL
528 if(aDSL.isUpperLeaningPoint(p))
530 if(aDSL.isLowerLeaningPoint(p))
543//-----------------------------------------------------------------------------
544template <typename TCoordinate, typename TInteger, unsigned short adjacency>
546DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
547DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
548::computeUnion(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
551 typename ArithmeticalDSL<TCoordinate,TInteger,adjacency>::Octant::first_type oc;
554 typedef DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency> DSS;
555 typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
557 DSS DSSres(Point(0,0)),DSS1(Point(0,0)),DSS2(Point(0,0));
559 // Test whether the two DSSs belong to the same octant or not
560 if(!(myDSL.sameOctant(aOther.dsl(), &oc)))
561 return ArithmeticalDSS(Point(0,0));
563 if(beforeOrEqual(myL,aOther.front()))
575 if(beforeOrEqual(DSS2.back(),DSS1.front()))
579 Vector step = DSS2.back() - DSS1.front();
580 Coordinate deviation = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(step[0], step[1]);
583 if(oc == 0 || oc == 3 || oc == 4 || oc == 7)
584 dir = Vector(0,1); else dir = Vector(1,0);
586 DGtal::IntegerComputer<Integer> ic;
587 if (deviation <= NumberTraits<Coordinate>::ONE || ic.dotProduct(step,dir) == NumberTraits<TInteger>::ZERO)
588 easyCase = true; else easyCase = false;
592 Point newUf, newUl, newLf, newLl;
593 std::vector<Point> LPoints;
595 // Test whether DSS2 belongs to the DSL supporting DSS1
596 std::vector<Point> UlPoints, LlPoints;
598 bool inDSL = DSS2.isInDSL(DSS1.dsl(),UlPoints,LlPoints,outP);
600 if(inDSL) // DSS2 in included in the DSL supporting DSS1
602 // if DSS1 is totally included into DSS2
603 if(beforeOrEqual(DSS2.back(),DSS1.back()))
607 if(easyCase) // the leaning points are updated in constant time
609 // if a upper leaning point of myDSL was found on aOther
610 if(!UlPoints.empty())
611 newUl = UlPoints.back(); // update the last upper leaning point with the furthest leaning point found
615 // if a lower leaning point of myDSL was found on aOther
616 if(!LlPoints.empty())
617 newLl = LlPoints.back(); // update the last lower leaning point with the furthest leaning point found
620 DSSres = ArithmeticalDSS<TCoordinate,TInteger,adjacency>(DSS1.a(),DSS1.b(),DSS1.back(),DSS2.front(),DSS1.Uf(),newUl,DSS1.Lf(),newLl);
623 else // leaning points may lie between DSS1 and DSS2
624 DSSres = Factory::createDSS(DSS1.a(),DSS1.b(),DSS1.back(),DSS2.front(),DSS1.Uf());
629 // If DSS2 does not belong to the DSL supporting DSS1, the
630 // remainder of the "out" point is used to infer the two candidate
633 // Candidate critical points are stored in Lf, Ll, Uf, Ul
634 if(DSS1.remainder(outP)>=(DSS1.dsl()).myUpperBound+1) // lower exterior point -> the slope decreases
639 else // remainder(outP) < 0: upper exterior point -> the slope increases
646 // Test whether DSS1 belongs to the supporting DSL of DSS2
647 std::vector<Point> UfPoints, LfPoints;
648 inDSL = DSS1.isInDSL(DSS2.dsl(),UfPoints,LfPoints,outP);
650 if(inDSL) // DSS1 in included in the DSL supporting DSS2
652 // if DSS1 is totally included into DSS2
653 if(beforeOrEqual(DSS2.back(),DSS1.back()))
659 // if a upper leaning point of DSS2 supporting DSL was found on DSS1
660 if(!UfPoints.empty())
661 newUf = UfPoints.front(); // update the first upper leaning point with the first leaning point found
665 // if a lower leaning point of DSS2 supporting DSL was found on DSS1
666 if(!LfPoints.empty())
667 newLf = LfPoints.front(); // update the first lower leaning point with the first leaning point found
670 DSSres = ArithmeticalDSS<TCoordinate,TInteger,adjacency>(DSS2.a(),DSS2.b(),DSS1.back(),DSS2.front(),newUf,DSS2.Ul(),newLf,DSS2.Ll());
673 DSSres = Factory::createDSS(DSS2.a(),DSS2.b(),DSS1.back(),DSS2.front(),DSS2.Ul());
678 // If DSS1 does not belong to the DSL supporting DSS2, the
679 // remainder of the "out" point is used to infer the two candidate
683 if(DSS2.remainder(outP)>=(DSS2.dsl()).myUpperBound+1) // lower exterior point -> the slope increases
688 else // remainder(outP) < 0: upper exterior point -> the slope decreases
694 // Four candidate critical points are known in Uf, Lf, Ul, Ll
695 // If only two different candidates, no solution
696 if(aUf==aUl && aLf==aLl)
697 return ArithmeticalDSS(Point(0,0));
699 // Otherwise, compute the exact critical points
700 Integer aB = aLl[0]-aLf[0];
701 Integer aA = aLl[1]-aLf[1];
702 Vector shiftRes = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::shift(aA,aB);
703 Integer aMu = -DSL::remainder(aA,aB,shiftRes) + 1 + DSL::remainder(aA,aB,aLf);
704 DSL DSLLower(aA,aB,aMu); // DSL defined by the two candidate lower
709 aMu = DSL::remainder(aA,aB,aUf);
710 DSL DSLUpper(aA,aB,aMu); // DSL defined by the two candidate upper
713 bool LPointsInDSL = DSLUpper.isInDSL(aLf) && DSLUpper.isInDSL(aLl);
714 bool UPointsInDSL = DSLLower.isInDSL(aUf) && DSLLower.isInDSL(aUl);
716 // DSS1 \cup DSS2 is not part of a DSL
717 if((aUf == aUl && !UPointsInDSL) || (aLf==aLl && !LPointsInDSL) || (!LPointsInDSL && !UPointsInDSL))
718 return ArithmeticalDSS(Point(0,0));
722 if(LPointsInDSL && UPointsInDSL)
723 { // four critical points
724 if(DSLUpper.a() != NumberTraits<TInteger>::ZERO || DSLUpper.b() != NumberTraits<TInteger>::ZERO) // a and b may be null
725 // if Ul and Uf are confounded
737 { // three critical points only
742 // test which one of Lf or Ll is closer from the DSL defined by the two
743 // upper leaning points
745 if(DSLUpper.remainder(aLf) > DSLUpper.remainder(aLl))
746 aLl = aLf; // Ll is not a critical point
748 aLf = aLl; // Lf is not a critical point
750 else // UPointsInDSL == true
754 // test which one of Uf or Ul is closer from the DSL defined by the two
755 // lower leaning points
756 if(DSLLower.remainder(aUf) < DSLLower.remainder(aUl))
763 // if the two DSSs are connected, or DSS1 last point and DSS2 first
764 // point have the same abscissa (or ordinate, depending on the octant) we are done.
766 return ArithmeticalDSS<TCoordinate, TInteger, adjacency>(aRes,bRes,DSS1.back(),DSS2.front(),aUf,aUl,aLf,aLl);
769 // Compute vectors of maximal and minimal slope
770 shiftRes = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::shift(aRes,bRes);
771 Vector Vup = (aLl-shiftRes) - aUf;
772 Vector Vlow = aUl - (aLf-shiftRes);
774 typedef DGtal::SternBrocot<TInteger, TInteger> SB; // type of Stern-Brocot tree
775 typedef typename SB::Fraction Fraction; // type of fraction
777 DGtal::functors::Abs<TInteger> absComputer;
779 Integer p = absComputer(Vlow[0])<absComputer(Vlow[1])?absComputer(Vlow[0]):absComputer(Vlow[1]);
780 Integer q = absComputer(Vlow[0])<absComputer(Vlow[1])?absComputer(Vlow[1]):absComputer(Vlow[0]);
781 Fraction low(p,q); // a fraction p/q is such that p>0, q>0 and p<=q
783 p = absComputer(Vup[0])<absComputer(Vup[1])?absComputer(Vup[0]):absComputer(Vup[1]);
784 q = absComputer(Vup[0])<absComputer(Vup[1])?absComputer(Vup[1]):absComputer(Vup[0]);
787 // Compute the fraction of smallest denominator between low and up
788 Fraction res = low.simplestFractionInBetween(up);
790 // Replace the result in the good octant
791 Integer aFinal,bFinal;
792 if(aRes < NumberTraits<TInteger>::ZERO)
793 aFinal = (absComputer(aRes)<absComputer(bRes))?-res.p():-res.q();
795 aFinal = (absComputer(aRes)<absComputer(bRes))?res.p():res.q();
797 if(bRes < NumberTraits<TInteger>::ZERO)
798 bFinal = (absComputer(aRes)<absComputer(bRes))?-res.q():-res.p();
800 bFinal = (absComputer(aRes)<absComputer(bRes))?res.q():res.p();
805 DSSres = Factory::createDSS(aFinal,bFinal,DSS1.back(),DSS2.front(),aUf);
809 typedef SimpleMatrix<Integer,2,2> Matrix;
811 A.setComponent(0,0,bFinal);A.setComponent(1,0,aFinal);A.setComponent(0,1,(aUf-aUl)[0]);A.setComponent(1,1,(aUf-aUl)[1]);
812 if(DGtal::SimpleMatrixSpecializations<Matrix,2,2>::determinant(A)>NumberTraits<TInteger>::ZERO)
813 DSSres = Factory::createDSS(aFinal,bFinal,DSS1.back(),DSS2.front(),aUf);
815 DSSres = Factory::createDSS(aFinal,bFinal,DSS1.back(),DSS2.front(),aUl);
822//-----------------------------------------------------------------------------
823template <typename TCoordinate, typename TInteger, unsigned short adjacency>
826DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::operator()(const Point& aPoint) const
828 return isInDSS(aPoint);
831//-----------------------------------------------------------------------------
832template <typename TCoordinate, typename TInteger, unsigned short adjacency>
834typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstIterator
835DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::begin() const
837 return ConstIterator(&myDSL, myF);
840//-----------------------------------------------------------------------------
841template <typename TCoordinate, typename TInteger, unsigned short adjacency>
843typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstIterator
844DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::end() const
846 ConstIterator it(&myDSL, myL);
851//-----------------------------------------------------------------------------
852template <typename TCoordinate, typename TInteger, unsigned short adjacency>
854typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstReverseIterator
855DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::rbegin() const
857 return ConstReverseIterator( end() );
860//-----------------------------------------------------------------------------
861template <typename TCoordinate, typename TInteger, unsigned short adjacency>
863typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstReverseIterator
864DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::rend() const
866 return ConstReverseIterator( begin() );
869//-----------------------------------------------------------------------------
870template <typename TCoordinate, typename TInteger, unsigned short adjacency>
873DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::selfDisplay ( std::ostream & out ) const
875 out << "[ArithmeticalDSS] ";
877 out << "from " << myF << " to " << myL << std::endl;
878 out << "upper leaning points: " << myUf << " " << myUl << std::endl;
879 out << "lower leaning points: " << myLf << " " << myLl;
882//-----------------------------------------------------------------------------
883template <typename TCoordinate, typename TInteger, unsigned short adjacency>
886DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
887isExtendableFront( const Point& aNewPoint ) const
889 Vector step = aNewPoint - myL;
890 Coordinate deviation = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(step[0], step[1]);
892 //if the two last points are confounded OK
893 if (deviation == NumberTraits<Coordinate>::ZERO)
896 //if the two last points are not connected KO
897 else if (deviation > NumberTraits<Coordinate>::ONE)
900 //if the first step does not exist yet OK
901 else if ( (myDSL.mySteps.first[0] == NumberTraits<Coordinate>::ZERO)
902 &&(myDSL.mySteps.first[1] == NumberTraits<Coordinate>::ZERO) )
905 //if the first step exists and
906 //if the second step does not exist
907 else if ( (myDSL.mySteps.second[0] == NumberTraits<Coordinate>::ZERO)
908 &&(myDSL.mySteps.second[1] == NumberTraits<Coordinate>::ZERO) )
910 //if the first step exists and is repeated OK
911 if (step == myDSL.mySteps.first)
915 //if the two steps are compatible OK
916 Vector v = step - myDSL.mySteps.first;
917 if ( DGtal::ArithmeticalDSLKernel<TCoordinate,DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::BackgroundAdjacency>
918 ::norm(v[0], v[1]) == NumberTraits<Coordinate>::ONE )
920 Integer r = remainder(aNewPoint);
921 //if weakly exterior on the left
922 if (r == myDSL.myLowerBound-NumberTraits<Integer>::ONE)
924 //if weakly exterior on the right
927 ASSERT(r == myDSL.myUpperBound+NumberTraits<Integer>::ONE);
935 //if the two steps are initialized
938 //if there are only two steps
939 if ( (step == myDSL.mySteps.first) || (step == myDSL.mySteps.second ) )
941 Integer r = remainder(aNewPoint);
943 //if strongly exterior KO
944 if ( (r < myDSL.myLowerBound-NumberTraits<Integer>::ONE)
945 ||(r > myDSL.myUpperBound+NumberTraits<Integer>::ONE) )
950 if (r == myDSL.myLowerBound)
952 else if (r == myDSL.myUpperBound)
954 else if (r == myDSL.myLowerBound-NumberTraits<Integer>::ONE)
956 else if (r == myDSL.myUpperBound+NumberTraits<Integer>::ONE)
962 //if there are more than two steps KO
968//-----------------------------------------------------------------------------
969template <typename TCoordinate, typename TInteger, unsigned short adjacency>
972DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
973isExtendableBack( const Point& aNewPoint ) const
975 return negate().isExtendableFront( aNewPoint );
978//-----------------------------------------------------------------------------
979template <typename TCoordinate, typename TInteger, unsigned short adjacency>
982DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
983extendFront( const Point& aNewPoint )
985 //true if the DSS can be extended to aNewPoint
989 //code that tells how to update the DSS
990 unsigned short int res = isExtendableFront(aNewPoint);
994 case 1: //first step init
995 myDSL.mySteps.first[0] = (aNewPoint[0] - myL[0]);
996 myDSL.mySteps.first[1] = (aNewPoint[1] - myL[1]);
997 myDSL.myA = myDSL.mySteps.first[1];
998 myDSL.myB = myDSL.mySteps.first[0];
999 myDSL.myLowerBound = remainder(myUf); //myUf doesn't change
1000 myDSL.myUpperBound = remainder(myLf); //myLf doesn't change
1004 myDSL.myShift = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::shift(myDSL.myA, myDSL.myB);
1006 case 2: //first step repeated
1011 case 3: //second step init on the left
1012 myDSL.myA = ((myUl[1] - myUf[1]) + (aNewPoint[1] - myL[1]));
1013 myDSL.myB = ((myUl[0] - myUf[0]) + (aNewPoint[0] - myL[0]));
1014 myDSL.myLowerBound = remainder(myUf); //myUf doesn't change
1015 myDSL.myUpperBound = remainder(myLl); //myLl doesn't change
1019 myDSL.mySteps = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::steps(myDSL.myA, myDSL.myB);
1020 myDSL.myShift = myDSL.mySteps.first-myDSL.mySteps.second;
1022 case 4: //second step init on the right
1023 myDSL.myA = ((myLl[1] - myLf[1]) + (aNewPoint[1] - myL[1]));
1024 myDSL.myB = ((myLl[0] - myLf[0]) + (aNewPoint[0] - myL[0]));
1025 myDSL.myLowerBound = remainder(myUl); //myUl doesn't change
1026 myDSL.myUpperBound = remainder(myLf); //myLf doesn't change
1030 myDSL.mySteps = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::steps(myDSL.myA, myDSL.myB);
1031 myDSL.myShift = myDSL.mySteps.first-myDSL.mySteps.second;
1033 case 5: //weakly interior on the left
1037 case 6: //weakly interior on the right
1041 case 7: //weakly exterior on the left
1042 myDSL.myA = (aNewPoint[1] - myUf[1]);
1043 myDSL.myB = (aNewPoint[0] - myUf[0]);
1044 myDSL.myLowerBound = remainder(myUf); //myUf doesn't change
1045 myDSL.myUpperBound = remainder(myLl); //myLl doesn't change
1050 case 8: //weakly exterior on the right
1051 myDSL.myA = (aNewPoint[1] - myLf[1]);
1052 myDSL.myB = (aNewPoint[0] - myLf[0]);
1053 myDSL.myLowerBound = remainder(myUl); //myUl doesn't change
1054 myDSL.myUpperBound = remainder(myLf); //myLf doesn't change
1059 case 9: //strongly interior
1069//-----------------------------------------------------------------------------
1070template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1073DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1074extendBack( const Point& aNewPoint )
1076 //call extendFront to the opposite
1077 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency> opposite = negate();
1078 bool flag = opposite.extendFront(aNewPoint);
1080 if (flag) //update '*this' if required
1081 *this = opposite.negate();
1086//-----------------------------------------------------------------------------
1087template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1090DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1093 //call retractBack to the opposite
1094 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency> opposite = negate();
1095 bool hasBeenRetracted = opposite.retractBack();
1097 if (hasBeenRetracted) //update '*this' if required
1098 *this = opposite.negate();
1100 return hasBeenRetracted;
1103//-----------------------------------------------------------------------------
1104template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1107DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1115 //next point computation
1116 Point next = *++begin();
1117 //if the next point is the last one
1120 *this = ArithmeticalDSS(next);
1126 //points used to update the DSS
1129 //leaning points / parameters update
1132 bezoutPoint = myUf + myDSL.myShift;
1133 if ( retractUpdateLeaningPoints( Vector(myDSL.myB, myDSL.myA),
1138 retractUpdateParameters( myLf - bezoutPoint );
1143 bezoutPoint = myLf - myDSL.myShift;
1144 if ( retractUpdateLeaningPoints( Vector(myDSL.myB, myDSL.myA),
1149 retractUpdateParameters( myUf - bezoutPoint );
1152 //first point update
1160//-----------------------------------------------------------------------------
1161template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1164DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1165retractUpdateLeaningPoints( const Vector& aDirection,
1166 const Point& aFirst,
1168 const Point& aBezout,
1169 const Point& aFirstAtOppositeSide,
1170 Point& aLastAtOppositeSide,
1171 Point& aFirstAtRemovalSide,
1172 const Point& aLastAtRemovalSide)
1174 if (aFirstAtOppositeSide == aLastAtOppositeSide)
1176 Vector newDirection = aFirstAtOppositeSide - aBezout;
1177 Coordinate k; //number of repetition of newDirection
1179 Vector toLastAtRemovalSide = aLastAtRemovalSide - aFirst;
1180 k = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(toLastAtRemovalSide[1], toLastAtRemovalSide[0])
1181 / DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(newDirection[1], newDirection[0]);
1182 aFirstAtRemovalSide = aLastAtRemovalSide - newDirection*k;
1184 Vector toLast = aLast - aFirstAtOppositeSide;
1185 k = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(toLast[1], toLast[0])
1186 / DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(newDirection[1], newDirection[0]);
1187 aLastAtOppositeSide = aFirstAtOppositeSide + newDirection*k;
1192 aFirstAtRemovalSide += aDirection;
1197//-----------------------------------------------------------------------------
1198template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1201DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1202retractUpdateParameters( const Vector& aNewDirection )
1204 //update of the slope
1205 myDSL.myA = aNewDirection[1];
1206 myDSL.myB = aNewDirection[0];
1208 myDSL.myLowerBound = remainder(myUf);
1209 myDSL.myUpperBound = remainder(myLf);
1210 //update of the steps and shift
1213 ASSERT( myUl == myLl );
1214 myDSL.mySteps = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::steps(myDSL.myA, myDSL.myB);
1215 myDSL.myShift = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::shift(myDSL.myA, myDSL.myB);
1219///////////////////////////////////////////////////////////////////////////////
1220// Drawing services //
1221///////////////////////////////////////////////////////////////////////////////
1222//-----------------------------------------------------------------------------
1223template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1225typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::PointD
1226DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1227project( const Point& aM, double aR ) const
1230 double aa = NumberTraits<Coordinate>::castToDouble(myDSL.myA);
1231 double bb = NumberTraits<Coordinate>::castToDouble(myDSL.myB);
1232 double xm = NumberTraits<Integer>::castToDouble(aM[0]);
1233 double ym = NumberTraits<Integer>::castToDouble(aM[1]);
1236 double d2 = ( aa * aa + bb * bb );
1237 double s = bb * xm + aa * ym;
1238 double xp = ( bb * s + aa * aR ) / d2;
1239 double yp = ( aa * s - bb * aR ) / d2;
1241 return PointD( xp, yp );
1244//-----------------------------------------------------------------------------
1245template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1247typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::PointD
1248DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1249project( const Point& aM, const Point& aP ) const
1251 double r = NumberTraits<Integer>::castToDouble(remainder(aP));
1252 return project(aM,r);
1255//-----------------------------------------------------------------------------
1256template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1259DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1262 return "ArithmeticalDSS";
1265///////////////////////////////////////////////////////////////////////////////
1267///////////////////////////////////////////////////////////////////////////////
1269//-----------------------------------------------------------------------------
1270template <typename TCoordinate, typename TInteger>
1272DGtal::StandardDSS4<TCoordinate, TInteger>::
1273StandardDSS4(const Coordinate& aA, const Coordinate& aB,
1274 const Point& aF, const Point& aL,
1275 const Point& aUf, const Point& aUl,
1276 const Point& aLf, const Point& aLl)
1277 : Super(aA, aB, aF, aL, aUf, aUl, aLf, aLl)
1280//-----------------------------------------------------------------------------
1281template <typename TCoordinate, typename TInteger>
1283DGtal::StandardDSS4<TCoordinate, TInteger>::
1284StandardDSS4(const Point& aF, const Point& aL,
1285 const bool& isOnTheUpperLine)
1286 : Super(aF, aL, isOnTheUpperLine)
1289//-----------------------------------------------------------------------------
1290template <typename TCoordinate, typename TInteger>
1292DGtal::StandardDSS4<TCoordinate, TInteger>::
1293StandardDSS4(const DSL& aDSL,
1294 const Point& aF, const Point& aL)
1295 : Super(aDSL, aF, aL)
1298//-----------------------------------------------------------------------------
1299template <typename TCoordinate, typename TInteger>
1301DGtal::StandardDSS4<TCoordinate, TInteger>::
1302StandardDSS4(const StandardDSS4<TCoordinate, TInteger>& aDSS,
1303 const Point& aF, const Point& aL)
1304 : Super(aDSS, aF, aL)
1307//-----------------------------------------------------------------------------
1308template <typename TCoordinate, typename TInteger>
1309template<typename Iterator>
1311DGtal::StandardDSS4<TCoordinate, TInteger>::
1312StandardDSS4(const Iterator& aItb, const Iterator& aIte)
1316//-----------------------------------------------------------------------------
1317template <typename TCoordinate, typename TInteger>
1319DGtal::StandardDSS4<TCoordinate, TInteger>::
1320StandardDSS4 ( const StandardDSS4 & aOther )
1324//-----------------------------------------------------------------------------
1325template <typename TCoordinate, typename TInteger>
1327DGtal::StandardDSS4<TCoordinate, TInteger>&
1328DGtal::StandardDSS4<TCoordinate, TInteger>::
1329operator= ( const StandardDSS4 & aOther )
1331 if (this != & aOther)
1332 Super::operator=( aOther );
1336//-----------------------------------------------------------------------------
1337template <typename TCoordinate, typename TInteger>
1339DGtal::NaiveDSS8<TCoordinate, TInteger>::
1340NaiveDSS8(const Coordinate& aA, const Coordinate& aB,
1341 const Point& aF, const Point& aL,
1342 const Point& aUf, const Point& aUl,
1343 const Point& aLf, const Point& aLl)
1344 : Super(aA, aB, aF, aL, aUf, aUl, aLf, aLl)
1347//-----------------------------------------------------------------------------
1348template <typename TCoordinate, typename TInteger>
1350DGtal::NaiveDSS8<TCoordinate, TInteger>::
1351NaiveDSS8(const Point& aF, const Point& aL,
1352 const bool& isOnTheUpperLine)
1353 : Super(aF, aL, isOnTheUpperLine)
1356//-----------------------------------------------------------------------------
1357template <typename TCoordinate, typename TInteger>
1359DGtal::NaiveDSS8<TCoordinate, TInteger>::
1360NaiveDSS8(const DSL& aDSL,
1361 const Point& aF, const Point& aL)
1362 : Super(aDSL, aF, aL)
1365//-----------------------------------------------------------------------------
1366template <typename TCoordinate, typename TInteger>
1368DGtal::NaiveDSS8<TCoordinate, TInteger>::
1369NaiveDSS8(const NaiveDSS8<TCoordinate, TInteger>& aDSS,
1370 const Point& aF, const Point& aL)
1371 : Super(aDSS, aF, aL)
1374//-----------------------------------------------------------------------------
1375template <typename TCoordinate, typename TInteger>
1376template<typename Iterator>
1378DGtal::NaiveDSS8<TCoordinate, TInteger>::
1379NaiveDSS8(const Iterator& aItb, const Iterator& aIte)
1383//-----------------------------------------------------------------------------
1384template <typename TCoordinate, typename TInteger>
1386DGtal::NaiveDSS8<TCoordinate, TInteger>::
1387NaiveDSS8 ( const NaiveDSS8 & aOther )
1391//-----------------------------------------------------------------------------
1392template <typename TCoordinate, typename TInteger>
1394DGtal::NaiveDSS8<TCoordinate, TInteger>&
1395DGtal::NaiveDSS8<TCoordinate, TInteger>::
1396operator= ( const NaiveDSS8 & aOther )
1398 if (this != & aOther)
1399 Super::operator=( aOther );
1403///////////////////////////////////////////////////////////////////////////////
1404// Implementation of inline functions //
1406template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1409DGtal::operator<< ( std::ostream & out,
1410 const ArithmeticalDSS<TCoordinate, TInteger, adjacency> & object )
1412 object.selfDisplay( out );
1417///////////////////////////////////////////////////////////////////////////////