DGtal  0.9.2
DSLSubsegment.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 DSLSubsegment.ih
19  * @author Isabelle Sivignon (\c isabelle.sivignon@gipsa-lab.grenoble-inp.fr )
20  * gipsa-lab Grenoble Images Parole Signal Automatique (CNRS, UMR 5216), CNRS, France
21  *
22  * @date 2012/07/17
23  *
24  * Implementation of inline methods defined in DSLSubsegment.h
25  *
26  */
27 
28 ///////////////////////////////////////////////////////////////////////////////
29 // IMPLEMENTATION of inline methods.
30 ///////////////////////////////////////////////////////////////////////////////
31 
32 //////////////////////////////////////////////////////////////////////////////
33 #include <cstdlib>
34 //////////////////////////////////////////////////////////////////////////////
35 
36 // #define DEBUG
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // Implementation of inline methods //
40 
41 
42 template <typename TInteger, typename TNumber>
43 inline
44 DGtal::DSLSubsegment<TInteger,TNumber>::RayC::RayC()
45 {}
46 
47 template <typename TInteger, typename TNumber>
48 inline
49 DGtal::DSLSubsegment<TInteger,TNumber>::RayC::RayC(const Integer x0, const Integer y0)
50 {
51  x = x0;
52  y = y0;
53 }
54 
55 template <typename TInteger, typename TNumber>
56 inline
57 DGtal::DSLSubsegment<TInteger,TNumber>::RayC::RayC(const Integer p , const Integer q, const Integer r, const Integer slope)
58 {
59  x = slope;
60  y = (r+p*x)/q;
61 }
62 
63 template <typename TInteger, typename TNumber>
64 inline
65 DGtal::DSLSubsegment<TInteger,TNumber>::RayC::~RayC()
66 {
67 }
68 
69 
70 template <typename TInteger, typename TNumber>
71 inline
72 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::intersectionVertical(Point &P, Vector &v, Integer n)
73 {
74  if(v[0] == 0)
75  return n;
76  else
77  {
78  Integer val = (n-P[0])/v[0];
79  if(((n-P[0])<0 || v[0] <0) && !((n-P[0])<0 && v[0] <0))
80  val--;
81 
82  return val;
83  }
84 }
85 
86 
87 
88 
89 template <typename TInteger, typename TNumber>
90 inline
91 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::intersection(Point &P, Vector &v, Vector &aL, Integer r)
92 {
93  Number num = P[1]*aL[0]-aL[1]*P[0]-r;
94  Number denom = aL[1]*v[0]-v[1]*aL[0];
95  // using parametric definition of lines and solving for t -> not more precise than the other computation...
96 
97  Integer val = num/denom;
98 
99  // integer rounding is not equivalent to the floor operation when num/denom is negative
100  if((num<0 ||denom <0) && !(num<0 && denom <0))
101  val--;
102 
103  return val;
104 }
105 
106 
107 
108 
109 template <typename TInteger, typename TNumber>
110 inline
111 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::intersection(Point &P, Vector &v, Number s)
112 {
113  // vector to be approximated = (1,s)
114 
115  if(v[0]==0)
116  {
117  Number val = (s*P[0] - P[1])/(v[1]*s);
118  if(ceil(val)-val <= myPrecision) // val is very close and slightly under an integer -> round to the integer
119  return (Integer) ceil(val);
120  else
121  return (Integer) floor(val);
122  }
123  else
124  {
125  Number num = -P[1]*v[0] - P[0]*(v[1] - s*v[0]) + v[1]*P[0];
126  Number denom = v[0]*(v[1] - s*v[0]);
127  Number val = num/denom;
128 
129  Number val3;
130  PointF Q;
131  Q[0] = P[0] + ceil(val)*v[0];
132  Q[1] = P[1] + ceil(val)*v[1];
133 
134  val3 = (Number) Q[0]*s - (Number) Q[1];
135 
136  if(std::abs(val3) <= myPrecision)
137  return (Integer) ceil(val);
138  else
139  return (Integer) floor(val);
140 
141  }
142 
143 }
144 
145 
146 template <typename TInteger, typename TNumber>
147 inline
148 void DGtal::DSLSubsegment<TInteger,TNumber>::update(Vector &u, Point& A, Vector &l, Integer r, Vector *v)
149 {
150  Point AA = A;
151  AA +=(*v);
152  Integer alpha = intersection(AA,u,l,r);
153 
154  u *= alpha;
155  (*v) += u;
156 }
157 
158 
159 
160 template <typename TInteger, typename TNumber>
161 inline
162 void DGtal::DSLSubsegment<TInteger,TNumber>::update(Vector &u, Point &A, Number s, Vector *v)
163 {
164  Point AA = A;
165  AA +=(*v);
166  Integer alpha = intersection(AA,u,s);
167 
168  u *= alpha;
169  (*v) += u;
170 }
171 
172 
173 
174 template <typename TInteger, typename TNumber>
175 inline
176 void DGtal::DSLSubsegment<TInteger,TNumber>::convexHullHarPeled(Vector &l, Integer n, Point *inf, Point *sup)
177 {
178  if(n >= l[0])
179  {
180  *inf = l;
181  *sup = l;
182  }
183 
184  if(l[1] == 0)
185  {
186  *inf = l;
187  *sup = Point(n,1);
188  }
189  else
190  {
191 
192  Point A = Point(0,1);
193  Point B = Point(1,0);
194 
195  int i=0;
196  bool ok = true;
197  Integer alpha1, alpha2, alpha;
198  DGtal::IntegerComputer<Integer> ic;
199 
200  while(ok)
201  {
202  if(i%2==0)
203  {
204  alpha1 = intersection(A,B,l,0);
205  alpha2 = intersectionVertical(A,B,n);
206  alpha = ic.min(alpha1,alpha2);
207  A = A + alpha*B;
208  if(A[0]*l[1]-A[1]*l[0]==0 || alpha == alpha2)
209  ok = false;
210  }
211  if(i%2==1)
212  {
213  alpha1 = intersection(B,A,l,0);
214  alpha2 = intersectionVertical(B,A,n);
215  alpha = ic.min(alpha1,alpha2);
216  B = B + alpha*A;
217  if(B[0]*l[1]-B[1]*l[0]==0 || alpha == alpha2)
218  ok = false;
219  }
220  i++;
221  }
222 
223  *inf = B;
224  *sup = A;
225  }
226 }
227 
228 
229 
230 
231 
232 template <typename TInteger, typename TNumber>
233 inline
234 void DGtal::DSLSubsegment<TInteger,TNumber>::convexHullApprox(Vector &l, Integer r, Integer n, Point *inf, Point *sup)
235 {
236 
237  //std::cout << "In Convex hull approx" << std::endl;
238 
239  //std::cout << l << " " << n << std::endl;
240 
241  if(n >= l[0])
242  {
243  *inf = l;
244  *sup = l;
245  }
246  else
247  if(l[1] == 0)
248  {
249  *inf = l;
250  *sup = Point(n,1);
251  }
252  else
253  {
254 
255  Point A = Point(0,1);
256  Point B = Point(1,0);
257 
258  Vector u = Vector(1,-1);
259  Vector v = Vector(1,0);
260 
261  Vector normal = Vector(-l[1],l[0]);
262 
263  Integer alpha;
264 
265  bool ok = true;
266  bool found = false;
267  bool first = true;
268  DGtal::IntegerComputer<TInteger> ic;
269 
270  while(ok)
271  {
272  if(v[0]*normal[0]+v[1]*normal[1] ==0) // should never happen
273  {
274  ok = false;
275  found = true;
276  }
277  else
278  if(v[0]*normal[0]+v[1]*normal[1] < 0)
279  {
280  alpha = ic.min(intersection(A,v,l,r), intersectionVertical(A,v,n));
281 
282  if(alpha >= 1 || first)
283  {
284  Vector vv = v;
285  vv *= alpha;
286  A += vv;
287  u = B - A;
288  update(u,A,l,r,&v);
289 
290  }
291  else
292  ok = false;
293  first = false;
294  }
295  else
296  if(v[0]*normal[0]+v[1]*normal[1] > 0)
297  {
298  alpha = ic.min(intersection(B,v,l,r),intersectionVertical(B,v,n));
299  if(alpha >= 1)
300  {
301  Vector vv = v;
302  vv *= alpha;
303  B += vv;
304  u = B - A;
305  update(u,A,l,r,&v);
306  }
307  else
308  ok = false;
309  }
310  }
311 
312  if(found)
313  {
314  *sup = v;
315  *inf = v;
316  }
317  else
318  {
319  *sup = A;
320  *inf = B;
321  assert(A[0] != B[0] || A[1] != B[1]);
322  }
323  }
324 }
325 
326 template <typename TInteger, typename TNumber>
327 inline
328 void DGtal::DSLSubsegment<TInteger,TNumber>::convexHullApproxTwoPoints(Vector &l, Integer r, Integer n, Point *inf, Point *sup, Point *prevInf, Point *prevSup, bool inv)
329 {
330  Point A,B,Aold,Bold;
331  Vector u,v;
332 
333  if(inv && r==0)
334  {
335  A = Point(0,0);
336  B = Point(1,0);
337  u = Vector(1,0);
338  v = Vector(0,1);
339  Aold = Point(0,0);
340  Bold = B;
341  }
342  else
343  {
344  A = Point(0,1);
345  B = Point(0,0);
346  u = Vector(0,-1);
347  v = Vector(1,0);
348  Aold = A;
349  Bold = B;
350  }
351 
352  Vector vv;
353  Vector normal = Vector(-l[1],l[0]);
354 
355  Integer alpha;
356 
357  bool ok = true;
358  bool first = true;
359  DGtal::IntegerComputer<Integer> ic;
360 
361  while(ok)
362  {
363 
364  if(v[0]*normal[0]+v[1]*normal[1]< 0)
365  {
366  alpha = ic.min(intersection(A,v,l,r), intersectionVertical(A,v,n));
367 
368  if(alpha >= 1 || first)
369  {
370  Aold = A;
371  vv = v;
372  vv *= alpha;
373  A += vv;
374  //test if A is on the slope
375  if(A[0]*l[1]-A[1]*l[0]+r ==0)
376  {
377  ok = false;
378  }
379  else
380  {
381  u = B - A;
382  update(u,A,l,r,&v);
383  }
384  }
385  else
386  ok = false;
387 
388  first = false;
389  }
390  else
391  if(v[0]*normal[0]+v[1]*normal[1]> 0)
392  {
393  alpha = ic.min(intersection(B,v,l,r),intersectionVertical(B,v,n));
394  if(alpha >= 1)
395  {
396  Bold = B;
397  vv = v;
398  vv *= alpha;
399  B += vv;
400  // test if B is on the slope
401  if(B[0]*l[1]-B[1]*l[0] +r==0 )
402  {
403  ok = false;
404  }
405  else
406  {
407  u = B - A;
408  update(u,A,l,r,&v);
409  }
410  }
411  else
412  ok = false;
413  }
414  else
415  ok = false;
416  }
417 
418 
419  if(!inv)
420  {
421  // if A is on the line, A is a vertex for both the lower and the upper convex hulls
422  if(A[0]*l[1]-A[1]*l[0]+r ==0)
423  {
424  *sup = A;
425  *inf = A;
426  *prevSup = Aold;
427  *prevInf = B;
428  }
429  else
430 
431  if(B[0]*l[1]-B[1]*l[0]+r ==0)
432  {
433  *inf = B;
434  *prevInf = Bold;
435  if(B[0]>A[0])
436  {
437  *sup = B;
438  *prevSup = A;
439  }
440  else
441  {
442  *sup = A;
443  *prevSup = Aold;
444  }
445  }
446  else
447  {
448  *sup = A;
449  *inf = B;
450  *prevInf = Bold;
451  *prevSup = Aold;
452  }
453  }
454  else
455  {
456  if(B[0]*l[1]-B[1]*l[0]+r ==0)
457  {
458  *sup = B;
459  *inf = B;
460  *prevSup = A;
461  *prevInf = Bold;
462  }
463  else
464  if(A[0]*l[1]-A[1]*l[0]+r ==0)
465  {
466  *sup = A;
467  *prevSup = Aold;
468  if(A[0]>B[0])
469  {
470  *inf = A;
471  *prevInf = B;
472  }
473  else
474  {
475  *inf = B;
476  *prevInf = Bold;
477  }
478  }
479  else
480  {
481  *sup = A;
482  *inf = B;
483  *prevInf = Bold;
484  *prevSup = Aold;
485  }
486 
487 
488  }
489 
490 }
491 
492 
493 
494 
495 
496 template <typename TInteger, typename TNumber>
497 inline
498 void DGtal::DSLSubsegment<TInteger,TNumber>::convexHullApprox(Number s, Integer n, Point *inf, Point *sup)
499 {
500  // slope to be approximated = (1,s)
501 
502  //std::cout << "In Convex hull approx" << std::endl;
503 
504  Point A = Point(0,1);
505  Point B = Point(1,0);
506 
507  Vector u = Vector(1,-1);
508  Vector v = Vector(1,0);
509 
510  VectorF normal = VectorF(-s,1);
511 
512  Integer alpha;
513 
514  bool ok = true;
515  bool found = false;
516 
517  DGtal::IntegerComputer<Integer> ic;
518 
519 
520  while(ok)
521  {
522  if(std::abs(v[0]*normal[0]+v[1]*normal[1]) <= myPrecision)
523  {
524  ok = false;
525  found = true;
526  }
527  else
528  if(v[0]*normal[0]+v[1]*normal[1] < 0)
529  {
530  alpha = ic.min(intersection(A,v,s), intersectionVertical(A,v,n));
531  if(alpha >= 1)
532  {
533  Vector vv = v;
534  vv *= alpha;
535  A += vv;
536  // test if A is on the slope
537  if(std::abs(A[0]*normal[0]+A[1]*normal[1]) <= myPrecision)
538  {
539  ok = false;
540  v = A;
541  found = true;
542  }
543  else
544  {
545  u = B - A;
546  update(u,A,s,&v);
547  }
548  }
549  else
550  ok = false;
551  }
552  else
553  if(v[0]*normal[0]+v[1]*normal[1] > 0)
554  {
555  alpha = ic.min(intersection(B,v,s),intersectionVertical(B,v,n));
556  if(alpha >= 1)
557  {
558  Vector vv = v;
559  vv *= alpha;
560  B += vv;
561  // test if B is on the slope
562  if(std::abs(B[0]*normal[0]+B[1]*normal[1]) <= myPrecision)
563  {
564  ok = false;
565  v = B;
566  found = true;
567  }
568  else
569  {
570  u = B - A;
571  update(u,A,s,&v);
572  }
573  }
574  else
575  ok = false;
576  }
577  }
578 
579  if(found)
580  {
581  Integer g = ic.gcd(v[0],v[1]);
582  v[0] = v[0]/g;
583  v[1] = v[1]/g;
584  *sup = v;
585  *inf = v;
586  }
587  else
588  {
589  *sup = A;
590  *inf = B;
591  assert(A[0] != B[0] || A[1] != B[1]);
592  }
593 }
594 
595 
596 
597 
598 // Compute the term following f=p/q in the Farey Series of order n. We
599 // have -q'p+p'q = 1, q' max such that q'<=n
600 // Complexity of extendedEuclid
601 template <typename TInteger, typename TNumber>
602 inline
603 typename DGtal::DSLSubsegment<TInteger,TNumber>::Point DGtal::DSLSubsegment<TInteger,TNumber>::nextTermInFareySeriesEuclid(Integer fp, Integer fq, Integer n)
604 {
605  Integer u,v;
606  // u*p+v*q = 1
607 
608  DGtal::IntegerComputer<Integer> ic;
609  Point p;
610  p = ic.extendedEuclid(fp,fq,1);
611 
612  u = p[0];
613  v = p[1];
614 
615  Integer pp = v;
616  Integer qq = -u;
617 
618  pp = pp + floor((n+u)/fq)*fp;
619  qq = qq + floor((n+u)/fq)*fq;
620 
621  return Point(qq,pp);
622 }
623 
624 
625 
626 template <typename TInteger, typename TNumber>
627 inline
628 typename DGtal::DSLSubsegment<TInteger,TNumber>::RayC DGtal::DSLSubsegment<TInteger,TNumber>::rayOfHighestSlope(Integer p, Integer q, Integer r, Integer smallestSlope, Integer n)
629  {
630  return RayC(p,q,r,n-(n-smallestSlope)%q);
631  }
632 
633 
634 
635 template <typename TInteger, typename TNumber>
636 inline
637 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::slope(Integer p, Integer q, Integer r, Number a, Number b, Number mu)
638 {
639  BOOST_CONCEPT_ASSERT((concepts::CInteger<TNumber>));
640  return (Integer) ceil((LongDoubleType) (r*b-mu*q)/(-p*b+a*q));
641 }
642 
643 template <typename TInteger, typename TNumber>
644 inline
645 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::slope(Integer p, Integer q, Integer r, Number alpha, Number beta)
646 {
647 
648  Number val = (r-beta*q)/(-p+alpha*q);
649 
650  // if the value is very close to a slightly smaller integer, we round to the close integer
651  if(val-floor(val) <= myPrecision)
652  return((Integer) floor(val));
653  else
654  return((Integer) ceil(val));
655 }
656 
657 
658 template <typename TInteger, typename TNumber>
659 inline
660 typename DGtal::DSLSubsegment<TInteger,TNumber>::Position DGtal::DSLSubsegment<TInteger,TNumber>::positionWrtRay(RayC &r, Number a, Number b, Number mu)
661 {
662  BOOST_CONCEPT_ASSERT((concepts::CInteger<TNumber>));
663  Integer v = -a*r.x + r.y*b - mu;
664 
665 
666  if(v == 0)
667  return ONTO;
668  else
669  if(v > 0)
670  return BELOW;
671 
672  return ABOVE;
673 
674 }
675 
676 
677 template <typename TInteger, typename TNumber>
678 inline
679 typename DGtal::DSLSubsegment<TInteger,TNumber>::Position DGtal::DSLSubsegment<TInteger,TNumber>::positionWrtRay(RayC &r, Number alpha, Number beta)
680 {
681  Number v = -alpha*r.x + r.y - beta;
682 
683  if(std::abs(v) <= myPrecision )
684  {
685  return ONTO;
686  }
687  else
688  if(v > 0)
689  return BELOW;
690 
691  return ABOVE;
692 }
693 
694 
695 
696 template <typename TInteger, typename TNumber>
697 inline
698 typename DGtal::DSLSubsegment<TInteger,TNumber>::RayC DGtal::DSLSubsegment<TInteger,TNumber>::smartRayOfSmallestSlope(Integer fp, Integer fq, Integer gp, Integer gq, Integer r)
699  {
700  // Version 1: using floor operator
701  Integer rr;
702  if(r == fq)
703  rr = gq;
704  else
705  rr = (r*gq)/fq;
706 
707  // Compute the slope of the line through (f=p/Q,r/q) and
708  // (g=p'/q',rr/q')
709  Integer x = (r*gq - rr*fq);
710 
711  Integer y = r*gp-rr*fp; // after simplification of the above formula
712 
713  return RayC(x,y);
714 
715  }
716 
717 
718 template <typename TInteger, typename TNumber>
719 inline
720 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::smartFirstDichotomy(Integer fp, Integer fq, Integer gp, Integer gq, Number a, Number b, Number mu, Integer n, bool *flagRayFound)
721 {
722  BOOST_CONCEPT_ASSERT((concepts::CInteger<Number>));
723 
724  RayC myRay;
725  Position myPosition;
726  Integer r = 0;
727  Integer lup = fq;
728  Integer ldown = 0;
729 
730  *flagRayFound = false;
731 
732  while((lup>=2 || ldown >=2) && !*flagRayFound)
733  {
734  myRay = smartRayOfSmallestSlope(fp,fq,gp,gq,r);
735  myPosition = positionWrtRay(myRay,a,b,mu);
736 
737 
738 // #ifdef DEBUG
739 // std::cout << "height " << r << " Ray " << myRay.x << " " << myRay.y << std::endl;
740 // std::cout << "myPosition " << myPosition << std::endl;
741 // #endif
742 
743 
744  if(myPosition == ONTO)
745  *flagRayFound = true;
746  else
747  if(myPosition == ABOVE)
748  {
749  r = r + ((lup%2 == 0)?lup/2:(lup-1)/2);
750  ldown = ((lup%2 ==0)?lup/2:(lup-1)/2);
751  lup = ((lup%2==0)?lup/2:(lup+1)/2);
752  }
753  else
754  {
755  r = r - ((ldown%2==0)?ldown/2:(ldown+1)/2);
756  lup = ((ldown%2==0)?ldown/2:(ldown+1)/2);
757  ldown = ((ldown%2==0)?ldown/2:(ldown-1)/2);
758  }
759 
760 
761  }
762 
763  // If the point is not on a ray of smallest slope
764  if(!*flagRayFound)
765  {
766  myRay = smartRayOfSmallestSlope(fp,fq,gp,gq,r);
767  myPosition = positionWrtRay(myRay,a,b,mu);
768 
769  if(myPosition != ONTO)
770  {
771  if(myPosition == ABOVE)
772  r++;
773 
774  RayC SteepestRay = rayOfHighestSlope(fp, fq,r,(smartRayOfSmallestSlope(fp,fq,gp,gq,r)).x,n);
775  Position pos = positionWrtRay(SteepestRay,a,b,mu);
776 
777  if(pos==BELOW)
778  {
779  r--;
780  *flagRayFound = true;
781 
782  }
783  }
784  else
785  *flagRayFound = true;
786 
787  }
788 
789 // #ifdef DEBUG
790 // std::cout << r << std::endl;
791 // #endif
792 
793  return r;
794  }
795 
796 
797 template <typename TInteger, typename TNumber>
798 inline
799 typename DGtal::DSLSubsegment<TInteger,TNumber>::Integer DGtal::DSLSubsegment<TInteger,TNumber>::smartFirstDichotomy(Integer fp, Integer fq, Integer gp, Integer gq, Number alpha, Number beta, Integer n, bool *flagRayFound)
800 {
801  RayC myRay;
802  Position myPosition;
803  Integer r = 0;
804  Integer lup = fq;
805  Integer ldown = 0;
806 
807  *flagRayFound = false;
808 
809 // #ifdef DEBUG
810 // std::cout << fp << " " << fq << " " << gp << " " << gq << std::endl;
811 // #endif
812 
813 
814  while((lup>=2 || ldown >=2) && !*flagRayFound)
815  {
816  myRay = smartRayOfSmallestSlope(fp,fq,gp,gq,r);
817  myPosition = positionWrtRay(myRay,alpha,beta);
818 
819 
820  // #ifdef DEBUG
821  // std::cout << "height " << r << " Ray " << myRay.x << " " << myRay.y << std::endl;
822  // std::cout << "myPosition " << myPosition << std::endl;
823  // #endif
824 
825 
826  if(myPosition == ONTO)
827  *flagRayFound = true;
828  else
829  if(myPosition == ABOVE)
830  {
831  r = r + ((lup%2 == 0)?lup/2:(lup-1)/2);
832  ldown = ((lup%2 ==0)?lup/2:(lup-1)/2);
833  lup = ((lup%2==0)?lup/2:(lup+1)/2);
834  }
835  else
836  {
837  r = r - ((ldown%2==0)?ldown/2:(ldown+1)/2);
838  lup = ((ldown%2==0)?ldown/2:(ldown+1)/2);
839  ldown = ((ldown%2==0)?ldown/2:(ldown-1)/2);
840  }
841 
842  }
843 
844  // If the point is not on a ray of smallest slope
845  if(!*flagRayFound)
846  {
847  myRay = smartRayOfSmallestSlope(fp,fq,gp,gq,r);
848  myPosition = positionWrtRay(myRay,alpha,beta);
849 
850  if(myPosition != ONTO)
851  {
852  if(myPosition == ABOVE)
853  r++;
854 
855  RayC SteepestRay = rayOfHighestSlope(fp, fq,r,(smartRayOfSmallestSlope(fp,fq,gp,gq,r)).x,n);
856  Position pos = positionWrtRay(SteepestRay,alpha,beta);
857 
858  if(pos == BELOW)
859  {
860  r--;
861  *flagRayFound = true;
862  }
863  }
864  else
865  *flagRayFound = true;
866 
867  }
868 
869  return r;
870 }
871 
872 
873 
874 template <typename TInteger, typename TNumber>
875 inline
876 typename DGtal::DSLSubsegment<TInteger,TNumber>::RayC DGtal::DSLSubsegment<TInteger,TNumber>::localizeRay(Integer fp, Integer fq, Integer gp, Integer gq, Integer r, Number a, Number b, Number mu, Integer n)
877 {
878  boost::ignore_unused_variable_warning( n );
879  BOOST_CONCEPT_ASSERT((concepts::CInteger<Number>));
880  Number alpha = slope(fp, fq,r,a,b,mu);
881 
882 
883  Integer smallestSlope = smartRayOfSmallestSlope(fp,fq,gp,gq,r).x;
884 
885 
886  if(alpha%(fq) == smallestSlope)
887  {
888  return RayC(fp,fq,r,alpha);
889  }
890  else
891  if(alpha%(fq) < smallestSlope)
892  {
893  return RayC(fp,fq,r,alpha + smallestSlope - alpha%(fq));
894  }
895  else
896  // alphaInt%(f.q()) > smallestSlope
897  {
898  return RayC(fp,fq,r,alpha - (alpha%(fq) - smallestSlope) + fq);
899  }
900  }
901 
902 // With a dichotomy
903 template <typename TInteger, typename TNumber>
904 inline
905 typename DGtal::DSLSubsegment<TInteger,TNumber>::RayC DGtal::DSLSubsegment<TInteger,TNumber>::localizeRay(Integer fp, Integer fq, Integer gp, Integer gq, Integer r, Number alpha, Number beta, Integer n)
906 {
907 
908  Integer smallestSlope = smartRayOfSmallestSlope(fp,fq,gp,gq,r).x;
909 
910  Integer kmax = floor((n-smallestSlope)/fq) +1;
911 
912  Integer k = 0;
913  Integer lup = 0;
914  Integer ldown = kmax;
915 
916  Position myPosition;
917 
918  RayC aRay;
919  bool found = false;
920 
921  while((lup>=2 || ldown >=2) && !found)
922  {
923  aRay = RayC(fp,fq,r,smallestSlope+k*fq);
924  myPosition = positionWrtRay(aRay,alpha,beta);
925 
926 #ifdef DEBUG
927  std::cout << "k " << k << " Ray " << aRay.x << " " << aRay.y << std::endl;
928  std::cout << "myPosition " << myPosition << std::endl;
929 #endif
930 
931  if(myPosition == ONTO)
932  found = true;
933  else
934  if(myPosition == ABOVE)
935  {
936  k = k - ((lup%2 == 0)?lup/2:(lup-1)/2);
937  ldown = ((lup%2 ==0)?lup/2:(lup-1)/2);
938  lup = ((lup%2==0)?lup/2:(lup+1)/2);
939  }
940  else
941  {
942  k = k + ((ldown%2==0)?ldown/2:(ldown+1)/2);
943  lup = ((ldown%2==0)?ldown/2:(ldown+1)/2);
944  ldown = ((ldown%2==0)?ldown/2:(ldown-1)/2);
945  }
946  }
947 
948  if(!found)
949  {
950  aRay = RayC(fp,fq,r,smallestSlope+k*fq);
951  myPosition = positionWrtRay(aRay,alpha,beta);
952 
953  if(myPosition == ABOVE || myPosition == ONTO)
954  return aRay;
955  else
956  return RayC(fp,fq,r,smallestSlope+(k+1)*fq);
957  }
958  else
959  return aRay;
960 
961 }
962 
963 
964 
965 
966 template <typename TInteger, typename TNumber>
967 inline
968 typename DGtal::DSLSubsegment<TInteger,TNumber>::RayC DGtal::DSLSubsegment<TInteger,TNumber>::raySup(Integer fp, Integer fq, RayC r)
969  {
970  RayC rr;
971  // r is the highest ray
972  if(r.x - fq < 0)
973  return r;
974  else
975  {
976  rr.x = r.x - fq;
977  rr.y = r.y - fp;
978  }
979 
980  return rr;
981  }
982 
983 
984 template <typename TInteger, typename TNumber>
985 inline
986 void DGtal::DSLSubsegment<TInteger,TNumber>::shortFindSolution(Integer fp, Integer fq, Integer gp, Integer gq, RayC r, Integer n, Integer *resAlphaP, Integer *resAlphaQ, Integer *resBetaP) // resBetaQ = resAlphaQ
987 {
988  Point inf, sup , tmpsup, tmpinf;
989  DGtal::IntegerComputer<Integer> ic;
990 
991  //Compute the lower edge of the facet
992  if(fq <= ic.max(r.x,n-r.x))
993  {
994  inf[0] = fq;
995  inf[1] = fp;
996  }
997  else
998  convexHullApprox(Vector(fq,fp),0,ic.max(r.x,n-r.x),&inf,&tmpsup);
999 
1000  if(gq <= ic.max(r.x,n-r.x))
1001  {
1002  sup[0] = gq;
1003  sup[1] = gp;
1004  }
1005  else
1006  convexHullApprox(Vector(gq,gp),0,ic.max(r.x,n-r.x),&tmpinf,&sup);
1007 
1008 
1009  if(r.x-inf[0] < 0) // R is the ray of smallest slope in inf
1010  {
1011  *resAlphaP = inf[1];
1012  *resAlphaQ = inf[0];
1013  *resBetaP = -(*resAlphaP)*r.x+(*resAlphaQ)*r.y;
1014  }
1015  else
1016  if(r.x+sup[0] > n) // R is the ray of highest slope in sup
1017  {
1018  *resAlphaP = sup[1];
1019  *resAlphaQ = sup[0];
1020  *resBetaP = -(*resAlphaP)*r.x+(*resAlphaQ)*r.y;
1021  }
1022  else //the facet is upper triangular,
1023  {
1024  Integer g = ic.gcd(inf[1] + sup[1],inf[0]+sup[0]);
1025 
1026  *resAlphaP = (inf[1]+sup[1])/g;
1027  *resAlphaQ = (inf[0]+sup[0])/g;
1028  *resBetaP = -(*resAlphaP)*r.x+(*resAlphaQ)*r.y;
1029  }
1030 
1031 }
1032 
1033 
1034 
1035 
1036 template <typename TInteger, typename TNumber>
1037 inline
1038 void DGtal::DSLSubsegment<TInteger,TNumber>::findSolutionWithoutFractions(Integer fp, Integer fq, Integer gp, Integer gq, RayC r, Integer n, Integer *resAlphaP, Integer *resAlphaQ, Integer *resBetaP, bool found) // resBetaQ = resAlphaQ
1039 {
1040  Point inf, sup;
1041  DGtal::IntegerComputer<Integer> ic;
1042  RayC rr;
1043 
1044 
1045  if(found == false)
1046  // r is not the highest ray on A
1047  {
1048 
1049  if(gq <= ic.max(r.x,n-r.x))
1050  { // B is a multiple point, r is the lowest ray on B
1051  // (otherwise, there would be an intersection in [AB]
1052  // between the lowest ray on B and the ray above r on A.
1053  // Thus B is the solution
1054  *resAlphaP = gp;
1055  *resAlphaQ = gq;
1056  *resBetaP = -gp*r.x + r.y*gq;
1057  }
1058  else
1059  { // compute the ray juste above r
1060  rr = raySup(fp,fq,r);
1061  // compute the fraction following f in the Farey Series
1062  // given by the slope of rr
1063  if(ic.max(rr.x,n-rr.x)>ic.max(r.x,n-r.x))
1064  {
1065  Point next = nextTermInFareySeriesEuclid(fp,fq,ic.max(rr.x,n-rr.x));
1066  *resAlphaP = next[1];
1067  *resAlphaQ = next[0];
1068  *resBetaP = (Integer) -(*resAlphaP)*r.x+(*resAlphaQ)*r.y;
1069  }
1070  else
1071  {
1072  Point next = nextTermInFareySeriesEuclid(fp,fq,ic.max(r.x,n-r.x));
1073  *resAlphaP = next[1];
1074  *resAlphaQ = next[0];
1075  *resBetaP = (Integer) -(*resAlphaP)*r.x+(*resAlphaQ)*r.y ;
1076  }
1077  }
1078  }
1079  else
1080  {
1081  if(fq <= ic.max(r.x,n-r.x))
1082  { // A is a multiple point
1083  // A is the solution
1084  *resAlphaP = fp;
1085  *resAlphaQ = fq;
1086  *resBetaP = -fp*r.x +fq*r.y;
1087  }
1088  else
1089  if(gq <= ic.max(r.x,n-r.x) && (r.x + gq) >n)
1090  { // B is a multiple point and r is the lowest ray
1091  // B is the solution
1092  *resAlphaP = gp;
1093  *resAlphaQ = gq;
1094  *resBetaP = -gp*r.x +gq*r.y;
1095  }
1096  else
1097  {
1098  // A is not a multiple point. Check if the vertex above A
1099  // is a multiple point
1100  Integer h = -fp*r.x + r.y*fq +1;
1101  rr = smartRayOfSmallestSlope(fp,fq,gp,gq,h);
1102  if(fq<=ic.max(rr.x,n-rr.x))
1103  { //the vertex above A is a
1104  // multiple point -> A is the solution
1105  *resAlphaP = fp;
1106  *resAlphaQ = fq;
1107  *resBetaP = -fp*r.x +fq*r.y;
1108  }
1109  else
1110  {
1111  Vector v(fq,fp);
1112  convexHullHarPeled(v,ic.max(r.x,n-r.x),&inf,&sup);
1113  // Let C be the point on r with abscissa equal to inf.
1114  rr = raySup(inf[1],inf[0],r); // ray above r passing through C
1115  if(rr.x == r.x) // r is the highest ray passing through C
1116  { // C is the solution
1117  *resAlphaP = inf[1];
1118  *resAlphaQ = inf[0];
1119  *resBetaP = -inf[1]*r.x + inf[0]*r.y;
1120  }
1121  else
1122  {
1123  if(ic.max(rr.x,n-rr.x)>ic.max(r.x,n-r.x))
1124  {
1125  // the solution is given by the fraction following inf
1126  // in the Farey Series of order
1127  // max(x-inf.q(),n-(x.inf.q())), on the ray rr
1128 
1129  // we compute the mediant on inf and sup: if
1130  // the denominator is lower or equal to the
1131  // order given by the ray rr, then the
1132  // mediant is the solution. Otherwise, sup is
1133  // the solution.
1134  Integer g = ic.gcd(inf[1] + sup[1],inf[0]+sup[0]);
1135  if((inf[0]+sup[0])/g <= ic.max(rr.x,n-rr.x))
1136  {
1137  *resAlphaP = (inf[1]+sup[1])/g;
1138  *resAlphaQ = (inf[0]+sup[0])/g;
1139  }
1140  else
1141  {
1142  *resAlphaP = sup[1];
1143  *resAlphaQ = sup[0];
1144  }
1145  *resBetaP = -(*resAlphaP)*rr.x+(*resAlphaQ)*rr.y - 1;
1146  }
1147  else
1148  {
1149  *resAlphaP = sup[1];
1150  *resAlphaQ = sup[0];
1151  *resBetaP = -(*resAlphaP)*r.x+(*resAlphaQ)*r.y;
1152  }
1153  }
1154 
1155 
1156  }
1157 
1158  }
1159 
1160  }
1161 }
1162 
1163 
1164 // Constructor in the case of integer input parameters
1165 template <typename TInteger, typename TNumber>
1166 inline
1167 DGtal::DSLSubsegment<TInteger,TNumber>::DSLSubsegment(Number a, Number b, Number mu, Point &A, Point &B, std::string type)
1168 {
1169  try
1170  {
1171  if(type.compare("farey")==0)
1172  this->DSLSubsegmentFareyFan(a,b,mu,A,B);
1173  else
1174  if(type.compare("localCH")==0)
1175  this->DSLSubsegmentLocalCH(a,b,mu,A,B);
1176  else
1177  throw std::invalid_argument("ERROR: Unknow string to specify the algorithm. \"farey\" is used hereafter.");
1178  }
1179  catch(std::invalid_argument & e)
1180  {
1181  std::cerr << e.what() << std::endl;
1182  this->DSLSubsegmentFareyFan(a,b,mu,A,B);
1183  }
1184 }
1185 
1186 
1187 
1188 template <typename TInteger, typename TNumber>
1189 inline
1190 void DGtal::DSLSubsegment<TInteger,TNumber>::DSLSubsegmentFareyFan(Number a, Number b, Number mu, Point &A, Point &B)
1191 {
1192  Integer n = B[0] - A[0];
1193 
1194 
1195  if(n >= 2*b)
1196  {
1197  myA = a;
1198  myB = b;
1199  myMu = mu;
1200  //std::cout << "DSS parameters are DSL parameters." << std::endl;
1201  ////std::cout << " " << a << " " << b << " " << mu;
1202  }
1203  else
1204  {
1205  // A becomes the origin // mu must be between 0 and b
1206  mu += a*A[0] - A[1]*b;
1207  Point inf, sup;
1208 
1209  Point inf2, sup2;
1210 
1211  Integer fp,fq,gp,gq;
1212 
1213  if(b>n)
1214  {
1215  Vector v(b,a);
1216  convexHullHarPeled(v,n,&inf,&sup);
1217  fp = inf[1];
1218  fq = inf[0];
1219  gp = sup[1];
1220  gq = sup[0];
1221  }
1222  else
1223  {
1224  Point next = nextTermInFareySeriesEuclid(a,b,n);
1225  fp = a;
1226  fq = b;
1227  gp = next[1];
1228  gq = next[0];
1229  }
1230 
1231  bool found;
1232 
1233  // Find the height in the ladder
1234  // Returns the height h such that:
1235  // - param is in between the rays passing through the point (inf =
1236  // p/q, h/q)
1237  // ==> found is set to false
1238  // - or param is above the ray of smallest slope passing through
1239  // (inf = p/q, h/q) but below all the rays passing through (p/q,
1240  // h+1/q) ==> found is set to true
1241 
1242 
1243  Integer h = smartFirstDichotomy(fp,fq,gp,gq,a,b,mu,n,&found);
1244 
1245  RayC r;
1246 
1247 
1248  if(found)
1249  {
1250  r = smartRayOfSmallestSlope(fp,fq,gp,gq,h);
1251  }
1252  else
1253  {
1254  r = localizeRay(fp,fq,gp,gq,h,a,b,mu,n);
1255  }
1256 
1257  Integer resAlphaP=0, resAlphaQ=0, resBetaP=0;
1258  findSolutionWithoutFractions(fp,fq, gp, gq, r, n, &resAlphaP, &resAlphaQ, &resBetaP, found);
1259 
1260  myA = resAlphaP;
1261  myB = resAlphaQ;
1262  myMu = resBetaP - myA*A[0] + myB*A[1];
1263 
1264  }
1265 
1266 
1267 }
1268 
1269 
1270 
1271 template <typename TInteger, typename TNumber>
1272 inline
1273 void DGtal::DSLSubsegment<TInteger,TNumber>::lowerConvexHull(Vector &l, Integer mu, Point &A, Point &B,
1274  Point *prevInfL, Point *infL, Point *infR, Point *prevInfR)
1275 {
1276 
1277  Integer rA = l[1]*A[0]-l[0]*A[1] + mu;
1278  Integer rB = l[1]*B[0]-l[0]*B[1] + mu;
1279 
1280  Point prevSupL,supL,supR,prevSupR;
1281 
1282  // from left to right
1283  convexHullApproxTwoPoints(l,rA,B[0]-A[0],infL,&supL,prevInfL,&prevSupL,0); // computation of the left part of the convex hulls
1284 
1285 
1286  *infL += A;
1287  supL += A;
1288 
1289  *prevInfL +=A;
1290  prevSupL +=A;
1291 
1292 
1293  // from right to left
1294  if(rB !=0)
1295  {
1296  convexHullApproxTwoPoints(l,l[0]-rB,B[0]-A[0],infR,&supR,prevInfR,&prevSupR,1);
1297 
1298  Point B1 = B + Point(0,1);
1299 
1300  *infR = B1 - *infR;
1301  supR = B1 - supR;
1302 
1303  // swap inf and sup
1304  Point tmp;
1305  tmp = *infR;
1306  *infR = supR;
1307  supR = tmp;
1308 
1309 
1310  *prevInfR = B1 - *prevInfR;
1311  prevSupR = B1 - prevSupR;
1312 
1313  tmp = *prevInfR;
1314  *prevInfR = prevSupR;
1315  prevSupR = tmp;
1316  }
1317  else
1318  {
1319  convexHullApproxTwoPoints(l,0,B[0]-A[0],infR,&supR,prevInfR,&prevSupR,1);
1320 
1321  *infR = B - *infR;
1322  supR = B - supR;
1323 
1324  Point tmp;
1325  tmp = *infR;
1326  *infR = supR;
1327  supR = tmp;
1328 
1329 
1330  *prevInfR = B - *prevInfR;
1331  prevSupR = B - prevSupR;
1332 
1333  tmp = *prevInfR;
1334  *prevInfR = prevSupR;
1335  prevSupR = tmp;
1336 
1337  }
1338 
1339 
1340 
1341 }
1342 
1343 
1344 
1345 
1346 // Contructor calling the localConvexHull approx instead of the walk in the Farey Fan -> should be a new class.
1347 template <typename TInteger, typename TNumber>
1348 inline
1349 void DGtal::DSLSubsegment<TInteger,TNumber>::DSLSubsegmentLocalCH(Number a, Number b, Number mu, Point &A, Point &B)
1350 {
1351  Integer rB = a*B[0]-b*B[1] + mu;
1352 
1353  Vector l(b,a);
1354 
1355  Point infL,supL,infR,supR,prevInfL,prevSupL,prevInfR,prevSupR;
1356 
1357  lowerConvexHull(l,mu,A,B,&prevInfL,&infL,&infR,&prevInfR);
1358 
1359  Point pz(0,0);
1360  Point BA=B-A;
1361  lowerConvexHull(l,l[0]-rB-1,pz,BA,&prevSupR, &supR, &supL, &prevSupL);
1362 
1363  prevSupR = B+Point(0,1)-prevSupR;
1364  supR = B+Point(0,1)-supR;
1365  prevSupL = B+Point(0,1)-prevSupL;
1366  supL = B+Point(0,1)-supL;
1367 
1368  DGtal::IntegerComputer<Integer> ic;
1369 
1370  myA=0, myB=0, myMu=0;
1371 
1372  // we are left with four possible lines given by the couples (prevInfL, infL=infR) (infL = infR,prevInfR)
1373  // (prevSupL, supL = supR) (supL = supR, prevSupR)
1374  // the solution is given by the slope with bigger b
1375 
1376 
1377  Point inf[4];
1378 
1379  inf[0] = prevInfL;
1380  inf[1] = infL;
1381  inf[2] = infR;
1382  inf[3] = prevInfR;
1383 
1384 
1385 
1386  Point sup[4];
1387 
1388  sup[0] = prevSupL;
1389  sup[1] = supL;
1390  sup[2] = supR;
1391  sup[3] = prevSupR;
1392 
1393 
1394 
1395  Integer atmp, btmp;
1396  for(int i=0;i<3;i++)
1397  {
1398  btmp = inf[i+1][0] - inf[i][0];
1399  atmp = inf[i+1][1] - inf[i][1];
1400  Integer g;
1401  if(atmp==0 || btmp==0)
1402  g = ic.max((Integer) 1,ic.max(atmp,btmp));
1403  else
1404  g = ic.gcd(btmp,atmp);
1405  btmp = btmp/g;
1406  atmp = atmp/g;
1407  if(btmp > myB)
1408  {
1409  myB =btmp;
1410  myA =atmp;
1411  }
1412  }
1413 
1414 
1415  for(int i=0;i<3;i++)
1416  {
1417  btmp = sup[i+1][0] - sup[i][0];
1418  atmp = sup[i+1][1] - sup[i][1];
1419  Integer g;
1420  if(atmp==0 || btmp==0)
1421  g = ic.max((Integer) 1,ic.max(atmp,btmp));
1422  else
1423  g = ic.gcd(btmp,atmp);
1424  btmp = btmp/g;
1425  atmp = atmp/g;
1426  if(btmp > myB)
1427  {
1428  myB =btmp;
1429  myA =atmp;
1430  }
1431  }
1432 
1433 
1434 
1435  myMu = -myA*infL[0] + myB*infL[1];
1436 
1437 }
1438 
1439 
1440 
1441 template <typename TInteger, typename TNumber>
1442 inline
1443 DGtal::DSLSubsegment<TInteger,TNumber>::DSLSubsegment(Number alpha, Number beta,
1444  Point &A, Point &B, Number precision)
1445 {
1446  Integer n = B[0] - A[0];
1447 
1448  myPrecision = precision;
1449 
1450  // A becomes the origin
1451 
1452  beta +=alpha*A[0] - A[1];
1453 
1454  Point inf, sup;
1455  Integer fp,fq,gp,gq;
1456 
1457 
1458 #ifdef DEBUG
1459  std::cout << "\n \nDSLSubsegment with floats, n=" << n << " precision = " << myPrecision << std::endl;
1460  std::cout << std::setprecision(15) << "alpha = " << alpha << " beta = " << beta << std::endl;
1461 #endif
1462 
1463 
1464  convexHullApprox(alpha,n,&inf,&sup);
1465  fp = inf[1];
1466  fq = inf[0];
1467  gp = sup[1];
1468  gq = sup[0];
1469 
1470 
1471 #ifdef DEBUG
1472  std::cout << "f and g = " << fp << " " << fq << " " << gp << " " << gq << std::endl;
1473 #endif
1474 
1475 
1476  if(fp == gp && fq == gq)
1477  {
1478  Point next = nextTermInFareySeriesEuclid(fp,fq,n);
1479  gp = next[1];
1480  gq = next[0];
1481 
1482  }
1483 
1484 
1485  bool found;
1486 
1487  // Find the height in the ladder
1488  // Returns the height h such that:
1489  // - param is in between the rays passing through the point (inf =
1490  // p/q, h/q)
1491  // ==> found is set to false
1492  // - or param is above the ray of smallest slope passing through
1493  // (inf = p/q, h/q) but below all the rays passing through (p/q,
1494  // h+1/q) ==> found is set to true
1495 
1496 
1497  Integer h = smartFirstDichotomy(fp,fq,gp,gq,alpha,beta,n,&found);
1498 
1499  RayC r;
1500 
1501  if(found)
1502  r = smartRayOfSmallestSlope(fp,fq,gp,gq,h);
1503  else
1504  r = localizeRay(fp,fq,gp,gq,h,alpha,beta,n);
1505 
1506 
1507  Integer resAlphaP=0, resAlphaQ=0, resBetaP=0;
1508  findSolutionWithoutFractions(fp,fq, gp, gq, r, n, &resAlphaP, &resAlphaQ, &resBetaP, found);
1509 
1510 
1511  myA = resAlphaP;
1512  myB = resAlphaQ;
1513  myMu = resBetaP - myA*A[0] + myB*A[1];
1514 
1515 }
1516 
1517 
1518 
1519 //------------------------- Accessors --------------------------
1520 
1521 template <typename TInteger, typename TNumber>
1522 inline
1523 TInteger
1524 DGtal::DSLSubsegment<TInteger,TNumber>::getA() const {
1525  return myA;
1526 }
1527 
1528 
1529 template <typename TInteger, typename TNumber>
1530 inline
1531 TInteger
1532 DGtal::DSLSubsegment<TInteger,TNumber>::getB() const {
1533  return myB;
1534 }
1535 
1536 
1537 
1538 template <typename TInteger, typename TNumber>
1539 inline
1540 TInteger
1541 DGtal::DSLSubsegment<TInteger,TNumber>::getMu() const {
1542  return myMu;
1543 }