DGtal  0.9.2
ArithmeticalDSS.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 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
21  *
22  * @date 2013/06/28
23  *
24  * Implementation of inline methods defined in ArithmeticalDSS.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 ///////////////////////////////////////////////////////////////////////////////
35 // IMPLEMENTATION of inline methods.
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // ----------------------- Standard services ------------------------------
40 
41 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
42 inline
43 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
44 ::ArithmeticalDSS(const Point& aPoint)
45 :
46 myF(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) )
56 {
57 }
58 
59 //-----------------------------------------------------------------------------
60 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
61 inline
62 DGtal::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)
69  :
70  myF(aF), myL(aL),
71  myUf(aUf), myUl(aUl), myLf(aLf), myLl(aLl),
72  myDSL( aA, aB, aLowerBound, aUpperBound, aSteps, aShift )
73 {
74 }
75 
76 //-----------------------------------------------------------------------------
77 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
78 inline
79 DGtal::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)
84  :
85  myF(aF), myL(aL),
86  myUf(aUf), myUl(aUl), myLf(aLf), myLl(aLl),
87  myDSL( aDSL )
88 {
89 }
90 
91 //-----------------------------------------------------------------------------
92 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
93 inline
94 DGtal::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)
99  :
100  myF(aF), myL(aL),
101  myUf(aUf), myUl(aUl), myLf(aLf), myLl(aLl),
102  myDSL( aA, aB, DSL::remainder( aA, aB, aUf ) )
103 {
104 }
105 
106 
107 //-----------------------------------------------------------------------------
108 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
109 inline
110 DGtal::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 )
114 {
115  typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
116 
117  if (areOnTheUpperLine)
118  *this = Factory::createPattern(aF, aL);
119  else
120  *this = Factory::createReversedPattern(aF, aL);
121 }
122 
123 //-----------------------------------------------------------------------------
124 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
125 inline
126 DGtal::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 )
129 {
130  typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
131  *this = Factory::createSubsegment(aDSL, aF, aL);
132 }
133 
134 //-----------------------------------------------------------------------------
135 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
136 inline
137 DGtal::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 )
141 {
142  typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
143  *this = Factory::createSubsegment(aDSS, aF, aL);
144 }
145 
146 //-----------------------------------------------------------------------------
147 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
148 template <typename Iterator>
149 inline
150 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
151 ::ArithmeticalDSS(const Iterator& aItb, const Iterator& aIte) throw( InputException )
152  : myDSL( NumberTraits<Coordinate>::ZERO, NumberTraits<Coordinate>::ZERO, NumberTraits<Integer>::ZERO )
153 {
154  if (aItb != aIte)
155  {
156  Iterator it = aItb;
157  //construction from the begin iterator
158  Point p = *it;
159  *this = ArithmeticalDSS(p);
160  //extension
161  for (++it; it != aIte; ++it)
162  extendFront(*it);
163  }
164  else
165  {
166  throw InputException();
167  }
168 }
169 
170 //-----------------------------------------------------------------------------
171 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
172 inline
173 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
174 ::ArithmeticalDSS(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther)
175  :
176  myF(aOther.myF), myL(aOther.myL),
177  myUf(aOther.myUf), myUl(aOther.myUl), myLf(aOther.myLf), myLl(aOther.myLl),
178  myDSL(aOther.myDSL)
179 {}
180 
181 //-----------------------------------------------------------------------------
182 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
183 inline
184 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>&
185 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
186 ::operator=(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther)
187 {
188  if ( this != &aOther )
189  {
190  myF = aOther.myF;
191  myL = aOther.myL;
192  myUf = aOther.myUf;
193  myLf = aOther.myLf;
194  myUl = aOther.myUl;
195  myLl = aOther.myLl;
196  myDSL = aOther.myDSL;
197  }
198  return *this;
199 }
200 
201 //-----------------------------------------------------------------------------
202 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
203 inline
204 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
205 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
206 ::negate() const
207 {
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),
212  -myDSL.myShift);
213 }
214 
215 //-----------------------------------------------------------------------------
216 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
217 inline
218 bool
219 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
220 ::equalsTo(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
221 {
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)) );
229 }
230 
231 //-----------------------------------------------------------------------------
232 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
233 inline
234 bool
235 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
236 ::operator==(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
237 {
238  return ( equalsTo(aOther) || equalsTo(aOther.negate()) );
239 }
240 
241 //-----------------------------------------------------------------------------
242 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
243 inline
244 bool
245 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
246 ::operator!=(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
247 {
248  return !( operator==(aOther) );
249 }
250 
251 //-----------------------------------------------------------------------------
252 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
253 inline
254 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::~ArithmeticalDSS()
255 {
256 }
257 
258 //-----------------------------------------------------------------------------
259 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
260 inline
261 bool
262 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::isValid() const
263 {
264  return functions::checkAll(*this);
265 }
266 
267 //-----------------------------------------------------------------------------
268 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
269 inline
270 const typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::DSL&
271 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::dsl() const
272 {
273  return myDSL;
274 }
275 
276 //-----------------------------------------------------------------------------
277 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
278 inline
279 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Coordinate
280 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::a() const
281 {
282  return myDSL.a();
283 }
284 
285 //-----------------------------------------------------------------------------
286 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
287 inline
288 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Coordinate
289 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::b() const
290 {
291  return myDSL.b();
292 }
293 
294 //-----------------------------------------------------------------------------
295 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
296 inline
297 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
298 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::mu() const
299 {
300  return myDSL.mu();
301 }
302 
303 //-----------------------------------------------------------------------------
304 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
305 inline
306 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
307 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::omega() const
308 {
309  return myDSL.omega();
310 }
311 
312 //-----------------------------------------------------------------------------
313 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
314 inline
315 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Vector
316 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::shift() const
317 {
318  return myDSL.shift();
319 }
320 
321 //-----------------------------------------------------------------------------
322 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
323 inline
324 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Steps
325 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::steps() const
326 {
327  return myDSL.steps();
328 }
329 
330 //-----------------------------------------------------------------------------
331 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
332 inline
333 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
334 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::back() const
335 {
336  return myF;
337 }
338 
339 //-----------------------------------------------------------------------------
340 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
341 inline
342 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
343 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::front() const
344 {
345  return myL;
346 }
347 
348 //-----------------------------------------------------------------------------
349 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
350 inline
351 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
352 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Uf() const
353 {
354  return myUf;
355 }
356 
357 //-----------------------------------------------------------------------------
358 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
359 inline
360 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
361 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Ul() const
362 {
363  return myUl;
364 }
365 
366 //-----------------------------------------------------------------------------
367 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
368 inline
369 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
370 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Lf() const
371 {
372  return myLf;
373 }
374 
375 //-----------------------------------------------------------------------------
376 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
377 inline
378 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Point
379 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Ll() const
380 {
381  return myLl;
382 }
383 
384 //-----------------------------------------------------------------------------
385 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
386 inline
387 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
388 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::remainder(const Point& aPoint) const
389 {
390  return myDSL.remainder(aPoint);
391 }
392 
393 //-----------------------------------------------------------------------------
394 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
395 inline
396 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Integer
397 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::orthogonalPosition(const Point& aPoint) const
398 {
399  return myDSL.orthogonalPosition(aPoint);
400 }
401 
402 //-----------------------------------------------------------------------------
403 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
404 inline
405 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::Position
406 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::position(const Point& aPoint) const
407 {
408  return myDSL.position(aPoint);
409 }
410 
411 //-----------------------------------------------------------------------------
412 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
413 inline
414 bool
415 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::before(const Point& aP1, const Point& aP2) const
416 {
417  return myDSL.before(aP1, aP2);
418 }
419 
420 //-----------------------------------------------------------------------------
421 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
422 inline
423 bool
424 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::beforeOrEqual(const Point& aP1, const Point& aP2) const
425 {
426  return myDSL.beforeOrEqual(aP1, aP2);
427 }
428 
429 //-----------------------------------------------------------------------------
430 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
431 inline
432 bool
433 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::isInDSL(const Point& aPoint) const
434 {
435  return myDSL(aPoint);
436 }
437 
438 //-----------------------------------------------------------------------------
439 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
440 inline
441 bool
442 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::isInDSS(const Point& aPoint) const
443 {
444  if (isInDSL(aPoint))
445  {
446  Integer s = position(aPoint);
447  Integer s1 = position(myF);
448  Integer s2 = position(myL);
449 
450  if (s1 == s2)
451  return (aPoint == myF);
452  else if (s1 < s2)
453  return ( (s >= s1)&&(s <= s2) );
454  else
455  return ( (s >= s2)&&(s <= s1) );
456  }
457  else
458  return false;
459 }
460 
461 
462 //----------------------------------------------------------------------------
463 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
464 inline
465 bool
466 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>::isInDSL(const DSL& aDSL) const
467 {
468 
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)))
473  return false;
474 
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);
480 
481  typename std::vector<Point>::const_iterator it = LPoints.begin();
482  Point p;
483  bool inDSL = true;
484  while(it != LPoints.end() && inDSL)
485  {
486  p = *it;
487  if(!aDSL.isInDSL(p))
488  inDSL = false;
489  ++it;
490  }
491 
492  return inDSL;
493 }
494 
495 
496 
497 //----------------------------------------------------------------------------
498 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
499 inline
500 bool
501 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>::isInDSL(const DSL& aDSL, std::vector<Point> &Ulp, std::vector<Point> &Llp, Point &outP) const
502 {
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();
507  outP = Point(0,0);
508  if(!(myDSL.sameOctant(aDSL, &oc)))
509  return false;
510 
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);
516 
517  typename std::vector<Point>::const_iterator it = LPoints.begin();
518  Point p;
519  bool inDSL = true;
520  while(it != LPoints.end() && inDSL)
521  {
522  p = *it;
523  if(!aDSL.isInDSL(p))
524  inDSL = false;
525  else
526  {
527  // store the points of the DSS that are leaning points for the DSL
528  if(aDSL.isUpperLeaningPoint(p))
529  Ulp.push_back(p);
530  if(aDSL.isLowerLeaningPoint(p))
531  Llp.push_back(p);
532  ++it;
533  }
534  }
535 
536  if(!inDSL)
537  outP = *it;
538  return inDSL;
539 }
540 
541 
542 
543 //-----------------------------------------------------------------------------
544 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
545 inline
546 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
547 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>
548 ::computeUnion(const ArithmeticalDSS<TCoordinate, TInteger, adjacency>& aOther) const
549 {
550 
551  typename ArithmeticalDSL<TCoordinate,TInteger,adjacency>::Octant::first_type oc;
552  bool easyCase;
553 
554  typedef DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency> DSS;
555  typedef DGtal::ArithmeticalDSSFactory<TCoordinate, TInteger, adjacency> Factory;
556 
557  DSS DSSres(Point(0,0)),DSS1(Point(0,0)),DSS2(Point(0,0));
558 
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));
562 
563  if(beforeOrEqual(myL,aOther.front()))
564  {
565  DSS1 = DSS(*this);
566  DSS2 = DSS(aOther);
567  }
568  else
569  {
570  DSS1 = DSS(aOther);
571  DSS2 = DSS(*this);
572  }
573 
574  // Test connectivity
575  if(beforeOrEqual(DSS2.back(),DSS1.front()))
576  easyCase = true;
577  else
578  {
579  Vector step = DSS2.back() - DSS1.front();
580  Coordinate deviation = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(step[0], step[1]);
581 
582  Vector dir;
583  if(oc == 0 || oc == 3 || oc == 4 || oc == 7)
584  dir = Vector(0,1); else dir = Vector(1,0);
585 
586  DGtal::IntegerComputer<Integer> ic;
587  if (deviation <= NumberTraits<Coordinate>::ONE || ic.dotProduct(step,dir) == NumberTraits<TInteger>::ZERO)
588  easyCase = true; else easyCase = false;
589  }
590 
591 
592  Point newUf, newUl, newLf, newLl;
593  std::vector<Point> LPoints;
594 
595  // Test whether DSS2 belongs to the DSL supporting DSS1
596  std::vector<Point> UlPoints, LlPoints;
597  Point outP;
598  bool inDSL = DSS2.isInDSL(DSS1.dsl(),UlPoints,LlPoints,outP);
599 
600  if(inDSL) // DSS2 in included in the DSL supporting DSS1
601  {
602  // if DSS1 is totally included into DSS2
603  if(beforeOrEqual(DSS2.back(),DSS1.back()))
604  DSSres = DSS2;
605  else
606  {
607  if(easyCase) // the leaning points are updated in constant time
608  {
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
612  else
613  newUl = DSS1.Ul();
614 
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
618  else
619  newLl = DSS1.Ll();
620  DSSres = ArithmeticalDSS<TCoordinate,TInteger,adjacency>(DSS1.a(),DSS1.b(),DSS1.back(),DSS2.front(),DSS1.Uf(),newUl,DSS1.Lf(),newLl);
621 
622  }
623  else // leaning points may lie between DSS1 and DSS2
624  DSSres = Factory::createDSS(DSS1.a(),DSS1.b(),DSS1.back(),DSS2.front(),DSS1.Uf());
625  }
626  return DSSres;
627  }
628 
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
631  // critical points.
632  Point aUf, aLf;
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
635  {
636  aUf = DSS1.Ul();
637  aLf = DSS1.Lf();
638  }
639  else // remainder(outP) < 0: upper exterior point -> the slope increases
640  {
641  aUf = DSS1.Uf();
642  aLf = DSS1.Ll();
643  }
644 
645 
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);
649 
650  if(inDSL) // DSS1 in included in the DSL supporting DSS2
651  {
652  // if DSS1 is totally included into DSS2
653  if(beforeOrEqual(DSS2.back(),DSS1.back()))
654  DSSres = DSS2;
655  else
656  {
657  if(easyCase)
658  {
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
662  else
663  newUf = DSS2.Uf();
664 
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
668  else
669  newLf = DSS2.Lf();
670  DSSres = ArithmeticalDSS<TCoordinate,TInteger,adjacency>(DSS2.a(),DSS2.b(),DSS1.back(),DSS2.front(),newUf,DSS2.Ul(),newLf,DSS2.Ll());
671  }
672  else
673  DSSres = Factory::createDSS(DSS2.a(),DSS2.b(),DSS1.back(),DSS2.front(),DSS2.Ul());
674  }
675  return DSSres;
676  }
677 
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
680  // critical points.
681  Point aUl, aLl;
682 
683  if(DSS2.remainder(outP)>=(DSS2.dsl()).myUpperBound+1) // lower exterior point -> the slope increases
684  {
685  aUl = DSS2.Uf();
686  aLl = DSS2.Ll();
687  }
688  else // remainder(outP) < 0: upper exterior point -> the slope decreases
689  {
690  aUl = DSS2.Ul();
691  aLl = DSS2.Lf();
692  }
693 
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));
698 
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
705  // leaning points
706 
707  aB = aUl[0]-aUf[0];
708  aA = aUl[1]-aUf[1];
709  aMu = DSL::remainder(aA,aB,aUf);
710  DSL DSLUpper(aA,aB,aMu); // DSL defined by the two candidate upper
711  // leaning points
712 
713  bool LPointsInDSL = DSLUpper.isInDSL(aLf) && DSLUpper.isInDSL(aLl);
714  bool UPointsInDSL = DSLLower.isInDSL(aUf) && DSLLower.isInDSL(aUl);
715 
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));
719 
720  Integer aRes, bRes;
721 
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
726  {
727  aRes = DSLUpper.a();
728  bRes = DSLUpper.b();
729  }
730  else
731  {
732  aRes = DSLLower.a();
733  bRes = DSLLower.b();
734  }
735  }
736  else
737  { // three critical points only
738  if(LPointsInDSL)
739  {
740  aRes = DSLUpper.a();
741  bRes = DSLUpper.b();
742  // test which one of Lf or Ll is closer from the DSL defined by the two
743  // upper leaning points
744 
745  if(DSLUpper.remainder(aLf) > DSLUpper.remainder(aLl))
746  aLl = aLf; // Ll is not a critical point
747  else
748  aLf = aLl; // Lf is not a critical point
749  }
750  else // UPointsInDSL == true
751  {
752  aRes = DSLLower.a();
753  bRes = DSLLower.b();
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))
757  aUl = aUf;
758  else
759  aUf = aUl;
760  }
761  }
762 
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.
765  if(easyCase)
766  return ArithmeticalDSS<TCoordinate, TInteger, adjacency>(aRes,bRes,DSS1.back(),DSS2.front(),aUf,aUl,aLf,aLl);
767 
768 
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);
773 
774  typedef DGtal::SternBrocot<TInteger, TInteger> SB; // type of Stern-Brocot tree
775  typedef typename SB::Fraction Fraction; // type of fraction
776 
777  DGtal::functors::Abs<TInteger> absComputer;
778 
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
782 
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]);
785  Fraction up(p,q);
786 
787  // Compute the fraction of smallest denominator between low and up
788  Fraction res = low.simplestFractionInBetween(up);
789 
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();
794  else
795  aFinal = (absComputer(aRes)<absComputer(bRes))?res.p():res.q();
796 
797  if(bRes < NumberTraits<TInteger>::ZERO)
798  bFinal = (absComputer(aRes)<absComputer(bRes))?-res.q():-res.p();
799  else
800  bFinal = (absComputer(aRes)<absComputer(bRes))?res.q():res.p();
801 
802 
803  if(aUf == aUl)
804  {
805  DSSres = Factory::createDSS(aFinal,bFinal,DSS1.back(),DSS2.front(),aUf);
806  }
807  else
808  {
809  typedef SimpleMatrix<Integer,2,2> Matrix;
810  Matrix A; A.clear();
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);
814  else
815  DSSres = Factory::createDSS(aFinal,bFinal,DSS1.back(),DSS2.front(),aUl);
816  }
817 
818  return DSSres;
819 }
820 
821 
822 //-----------------------------------------------------------------------------
823 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
824 inline
825 bool
826 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::operator()(const Point& aPoint) const
827 {
828  return isInDSS(aPoint);
829 }
830 
831 //-----------------------------------------------------------------------------
832 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
833 inline
834 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstIterator
835 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::begin() const
836 {
837  return ConstIterator(&myDSL, myF);
838 }
839 
840 //-----------------------------------------------------------------------------
841 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
842 inline
843 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstIterator
844 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::end() const
845 {
846  ConstIterator it(&myDSL, myL);
847  it++;
848  return it;
849 }
850 
851 //-----------------------------------------------------------------------------
852 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
853 inline
854 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstReverseIterator
855 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::rbegin() const
856 {
857  return ConstReverseIterator( end() );
858 }
859 
860 //-----------------------------------------------------------------------------
861 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
862 inline
863 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::ConstReverseIterator
864 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::rend() const
865 {
866  return ConstReverseIterator( begin() );
867 }
868 
869 //-----------------------------------------------------------------------------
870 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
871 inline
872 void
873 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::selfDisplay ( std::ostream & out ) const
874 {
875  out << "[ArithmeticalDSS] ";
876  out << myDSL;
877  out << "from " << myF << " to " << myL << std::endl;
878  out << "upper leaning points: " << myUf << " " << myUl << std::endl;
879  out << "lower leaning points: " << myLf << " " << myLl;
880 }
881 
882 //-----------------------------------------------------------------------------
883 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
884 inline
885 unsigned short int
886 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
887 isExtendableFront( const Point& aNewPoint ) const
888 {
889  Vector step = aNewPoint - myL;
890  Coordinate deviation = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::norm(step[0], step[1]);
891 
892  //if the two last points are confounded OK
893  if (deviation == NumberTraits<Coordinate>::ZERO)
894  return 9;
895 
896  //if the two last points are not connected KO
897  else if (deviation > NumberTraits<Coordinate>::ONE)
898  return 0;
899 
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) )
903  return 1;
904 
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) )
909  {
910  //if the first step exists and is repeated OK
911  if (step == myDSL.mySteps.first)
912  return 2;
913  else
914  {
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 )
919  {
920  Integer r = remainder(aNewPoint);
921  //if weakly exterior on the left
922  if (r == myDSL.myLowerBound-NumberTraits<Integer>::ONE)
923  return 3;
924  //if weakly exterior on the right
925  else
926  {
927  ASSERT(r == myDSL.myUpperBound+NumberTraits<Integer>::ONE);
928  return 4;
929  }
930  }
931  else
932  return 0;
933  }
934  }
935  //if the two steps are initialized
936  else
937  {
938  //if there are only two steps
939  if ( (step == myDSL.mySteps.first) || (step == myDSL.mySteps.second ) )
940  {
941  Integer r = remainder(aNewPoint);
942 
943  //if strongly exterior KO
944  if ( (r < myDSL.myLowerBound-NumberTraits<Integer>::ONE)
945  ||(r > myDSL.myUpperBound+NumberTraits<Integer>::ONE) )
946  return 0;
947  //otherwise OK
948  else
949  {
950  if (r == myDSL.myLowerBound)
951  return 5;
952  else if (r == myDSL.myUpperBound)
953  return 6;
954  else if (r == myDSL.myLowerBound-NumberTraits<Integer>::ONE)
955  return 7;
956  else if (r == myDSL.myUpperBound+NumberTraits<Integer>::ONE)
957  return 8;
958  else
959  return 9;
960  }
961  }
962  //if there are more than two steps KO
963  else
964  return 0;
965  }
966 }
967 
968 //-----------------------------------------------------------------------------
969 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
970 inline
971 unsigned short int
972 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
973 isExtendableBack( const Point& aNewPoint ) const
974 {
975  return negate().isExtendableFront( aNewPoint );
976 }
977 
978 //-----------------------------------------------------------------------------
979 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
980 inline
981 bool
982 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
983 extendFront( const Point& aNewPoint )
984 {
985  //true if the DSS can be extended to aNewPoint
986  //false otherwise
987  bool flag = true;
988 
989  //code that tells how to update the DSS
990  unsigned short int res = isExtendableFront(aNewPoint);
991 
992  switch (res)
993  {
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
1001  myL = aNewPoint;
1002  myUl = aNewPoint;
1003  myLl = aNewPoint;
1004  myDSL.myShift = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::shift(myDSL.myA, myDSL.myB);
1005  break;
1006  case 2: //first step repeated
1007  myL = aNewPoint;
1008  myUl = aNewPoint;
1009  myLl = aNewPoint;
1010  break;
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
1016  myL = aNewPoint;
1017  myUl = aNewPoint;
1018  myLf = myLl;
1019  myDSL.mySteps = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::steps(myDSL.myA, myDSL.myB);
1020  myDSL.myShift = myDSL.mySteps.first-myDSL.mySteps.second;
1021  break;
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
1027  myL = aNewPoint;
1028  myLl = aNewPoint;
1029  myUf = myUl;
1030  myDSL.mySteps = DGtal::ArithmeticalDSLKernel<TCoordinate,adjacency>::steps(myDSL.myA, myDSL.myB);
1031  myDSL.myShift = myDSL.mySteps.first-myDSL.mySteps.second;
1032  break;
1033  case 5: //weakly interior on the left
1034  myL = aNewPoint;
1035  myUl = aNewPoint;
1036  break;
1037  case 6: //weakly interior on the right
1038  myL = aNewPoint;
1039  myLl = aNewPoint;
1040  break;
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
1046  myL = aNewPoint;
1047  myUl = aNewPoint;
1048  myLf = myLl;
1049  break;
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
1055  myL = aNewPoint;
1056  myLl = aNewPoint;
1057  myUf = myUl;
1058  break;
1059  case 9: //strongly interior
1060  myL = aNewPoint;
1061  break;
1062  default:
1063  flag = false;
1064  }
1065 
1066  return flag;
1067 }
1068 
1069 //-----------------------------------------------------------------------------
1070 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1071 inline
1072 bool
1073 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1074 extendBack( const Point& aNewPoint )
1075 {
1076  //call extendFront to the opposite
1077  DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency> opposite = negate();
1078  bool flag = opposite.extendFront(aNewPoint);
1079 
1080  if (flag) //update '*this' if required
1081  *this = opposite.negate();
1082 
1083  return flag;
1084 }
1085 
1086 //-----------------------------------------------------------------------------
1087 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1088 inline
1089 bool
1090 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1091 retractFront()
1092 {
1093  //call retractBack to the opposite
1094  DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency> opposite = negate();
1095  bool hasBeenRetracted = opposite.retractBack();
1096 
1097  if (hasBeenRetracted) //update '*this' if required
1098  *this = opposite.negate();
1099 
1100  return hasBeenRetracted;
1101 }
1102 
1103 //-----------------------------------------------------------------------------
1104 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1105 inline
1106 bool
1107 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1108 retractBack()
1109 {
1110  if (myF == myL)
1111  return false;
1112  else
1113  {
1114 
1115  //next point computation
1116  Point next = *++begin();
1117  //if the next point is the last one
1118  if (next == myL)
1119  {
1120  *this = ArithmeticalDSS(next);
1121  return true;
1122  }
1123  //otherwise
1124  else
1125  {
1126  //points used to update the DSS
1127  Point bezoutPoint;
1128 
1129  //leaning points / parameters update
1130  if (myF == myUf)
1131  {
1132  bezoutPoint = myUf + myDSL.myShift;
1133  if ( retractUpdateLeaningPoints( Vector(myDSL.myB, myDSL.myA),
1134  next, myL,
1135  bezoutPoint,
1136  myLf, myLl,
1137  myUf, myUl ) )
1138  retractUpdateParameters( myLf - bezoutPoint );
1139 
1140  }
1141  if (myF == myLf)
1142  {
1143  bezoutPoint = myLf - myDSL.myShift;
1144  if ( retractUpdateLeaningPoints( Vector(myDSL.myB, myDSL.myA),
1145  next, myL,
1146  bezoutPoint,
1147  myUf, myUl,
1148  myLf, myLl ) )
1149  retractUpdateParameters( myUf - bezoutPoint );
1150  }
1151 
1152  //first point update
1153  myF = next;
1154 
1155  return true;
1156  }
1157  }
1158 }
1159 
1160 //-----------------------------------------------------------------------------
1161 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1162 inline
1163 bool
1164 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1165 retractUpdateLeaningPoints( const Vector& aDirection,
1166  const Point& aFirst,
1167  const Point& aLast,
1168  const Point& aBezout,
1169  const Point& aFirstAtOppositeSide,
1170  Point& aLastAtOppositeSide,
1171  Point& aFirstAtRemovalSide,
1172  const Point& aLastAtRemovalSide)
1173 {
1174  if (aFirstAtOppositeSide == aLastAtOppositeSide)
1175  {
1176  Vector newDirection = aFirstAtOppositeSide - aBezout;
1177  Coordinate k; //number of repetition of newDirection
1178 
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;
1183 
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;
1188  return true;
1189  }
1190  else
1191  {
1192  aFirstAtRemovalSide += aDirection;
1193  return false;
1194  }
1195 }
1196 
1197 //-----------------------------------------------------------------------------
1198 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1199 inline
1200 void
1201 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1202 retractUpdateParameters( const Vector& aNewDirection )
1203 {
1204  //update of the slope
1205  myDSL.myA = aNewDirection[1];
1206  myDSL.myB = aNewDirection[0];
1207  //intercepts
1208  myDSL.myLowerBound = remainder(myUf);
1209  myDSL.myUpperBound = remainder(myLf);
1210  //update of the steps and shift
1211  if (myUf == myLf)
1212  {
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);
1216  }
1217 }
1218 
1219 ///////////////////////////////////////////////////////////////////////////////
1220 // Drawing services //
1221 ///////////////////////////////////////////////////////////////////////////////
1222 //-----------------------------------------------------------------------------
1223 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1224 inline
1225 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::PointD
1226 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1227 project( const Point& aM, double aR ) const
1228 {
1229  //casts
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]);
1234 
1235  //computation
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;
1240 
1241  return PointD( xp, yp );
1242 }
1243 
1244 //-----------------------------------------------------------------------------
1245 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1246 inline
1247 typename DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::PointD
1248 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1249 project( const Point& aM, const Point& aP ) const
1250 {
1251  double r = NumberTraits<Integer>::castToDouble(remainder(aP));
1252  return project(aM,r);
1253 }
1254 
1255 //-----------------------------------------------------------------------------
1256 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1257 inline
1258 std::string
1259 DGtal::ArithmeticalDSS<TCoordinate, TInteger, adjacency>::
1260 className() const
1261 {
1262  return "ArithmeticalDSS";
1263 }
1264 
1265 ///////////////////////////////////////////////////////////////////////////////
1266 // Aliases //
1267 ///////////////////////////////////////////////////////////////////////////////
1268 
1269 //-----------------------------------------------------------------------------
1270 template <typename TCoordinate, typename TInteger>
1271 inline
1272 DGtal::StandardDSS4<TCoordinate, TInteger>::
1273 StandardDSS4(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)
1278 {}
1279 
1280 //-----------------------------------------------------------------------------
1281 template <typename TCoordinate, typename TInteger>
1282 inline
1283 DGtal::StandardDSS4<TCoordinate, TInteger>::
1284 StandardDSS4(const Point& aF, const Point& aL,
1285  const bool& isOnTheUpperLine)
1286  : Super(aF, aL, isOnTheUpperLine)
1287 {}
1288 
1289 //-----------------------------------------------------------------------------
1290 template <typename TCoordinate, typename TInteger>
1291 inline
1292 DGtal::StandardDSS4<TCoordinate, TInteger>::
1293 StandardDSS4(const DSL& aDSL,
1294  const Point& aF, const Point& aL)
1295  : Super(aDSL, aF, aL)
1296 {}
1297 
1298 //-----------------------------------------------------------------------------
1299 template <typename TCoordinate, typename TInteger>
1300 inline
1301 DGtal::StandardDSS4<TCoordinate, TInteger>::
1302 StandardDSS4(const StandardDSS4<TCoordinate, TInteger>& aDSS,
1303  const Point& aF, const Point& aL)
1304  : Super(aDSS, aF, aL)
1305 {}
1306 
1307 //-----------------------------------------------------------------------------
1308 template <typename TCoordinate, typename TInteger>
1309 template<typename Iterator>
1310 inline
1311 DGtal::StandardDSS4<TCoordinate, TInteger>::
1312 StandardDSS4(const Iterator& aItb, const Iterator& aIte) throw( InputException )
1313  : Super(aItb, aIte)
1314 {}
1315 
1316 //-----------------------------------------------------------------------------
1317 template <typename TCoordinate, typename TInteger>
1318 inline
1319 DGtal::StandardDSS4<TCoordinate, TInteger>::
1320 StandardDSS4 ( const StandardDSS4 & aOther )
1321  : Super( aOther )
1322 {}
1323 
1324 //-----------------------------------------------------------------------------
1325 template <typename TCoordinate, typename TInteger>
1326 inline
1327 DGtal::StandardDSS4<TCoordinate, TInteger>&
1328 DGtal::StandardDSS4<TCoordinate, TInteger>::
1329 operator= ( const StandardDSS4 & aOther )
1330 {
1331  if (this != & aOther)
1332  Super::operator=( aOther );
1333  return *this;
1334 }
1335 
1336 //-----------------------------------------------------------------------------
1337 template <typename TCoordinate, typename TInteger>
1338 inline
1339 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1340 NaiveDSS8(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)
1345 {}
1346 
1347 //-----------------------------------------------------------------------------
1348 template <typename TCoordinate, typename TInteger>
1349 inline
1350 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1351 NaiveDSS8(const Point& aF, const Point& aL,
1352  const bool& isOnTheUpperLine)
1353  : Super(aF, aL, isOnTheUpperLine)
1354 {}
1355 
1356 //-----------------------------------------------------------------------------
1357 template <typename TCoordinate, typename TInteger>
1358 inline
1359 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1360 NaiveDSS8(const DSL& aDSL,
1361  const Point& aF, const Point& aL)
1362  : Super(aDSL, aF, aL)
1363 {}
1364 
1365 //-----------------------------------------------------------------------------
1366 template <typename TCoordinate, typename TInteger>
1367 inline
1368 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1369 NaiveDSS8(const NaiveDSS8<TCoordinate, TInteger>& aDSS,
1370  const Point& aF, const Point& aL)
1371  : Super(aDSS, aF, aL)
1372 {}
1373 
1374 //-----------------------------------------------------------------------------
1375 template <typename TCoordinate, typename TInteger>
1376 template<typename Iterator>
1377 inline
1378 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1379 NaiveDSS8(const Iterator& aItb, const Iterator& aIte) throw( InputException )
1380  : Super(aItb, aIte)
1381 {}
1382 
1383 //-----------------------------------------------------------------------------
1384 template <typename TCoordinate, typename TInteger>
1385 inline
1386 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1387 NaiveDSS8 ( const NaiveDSS8 & aOther )
1388  : Super( aOther )
1389 {}
1390 
1391 //-----------------------------------------------------------------------------
1392 template <typename TCoordinate, typename TInteger>
1393 inline
1394 DGtal::NaiveDSS8<TCoordinate, TInteger>&
1395 DGtal::NaiveDSS8<TCoordinate, TInteger>::
1396 operator= ( const NaiveDSS8 & aOther )
1397 {
1398  if (this != & aOther)
1399  Super::operator=( aOther );
1400  return *this;
1401 }
1402 
1403 ///////////////////////////////////////////////////////////////////////////////
1404 // Implementation of inline functions //
1405 
1406 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
1407 inline
1408 std::ostream&
1409 DGtal::operator<< ( std::ostream & out,
1410  const ArithmeticalDSS<TCoordinate, TInteger, adjacency> & object )
1411 {
1412  object.selfDisplay( out );
1413  return out;
1414 }
1415 
1416 // //
1417 ///////////////////////////////////////////////////////////////////////////////
1418 
1419