DGtal 1.3.0
Loading...
Searching...
No Matches
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
42template <typename TInteger, typename TNumber>
43inline
44DGtal::DSLSubsegment<TInteger,TNumber>::RayC::RayC()
45{}
46
47template <typename TInteger, typename TNumber>
48inline
49DGtal::DSLSubsegment<TInteger,TNumber>::RayC::RayC(const Integer x0, const Integer y0)
50{
51 x = x0;
52 y = y0;
53}
54
55template <typename TInteger, typename TNumber>
56inline
57DGtal::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
63template <typename TInteger, typename TNumber>
64inline
65DGtal::DSLSubsegment<TInteger,TNumber>::RayC::~RayC()
66{
67}
68
69
70template <typename TInteger, typename TNumber>
71inline
72typename 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
89template <typename TInteger, typename TNumber>
90inline
91typename 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
109template <typename TInteger, typename TNumber>
110inline
111typename 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
146template <typename TInteger, typename TNumber>
147inline
148void 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
160template <typename TInteger, typename TNumber>
161inline
162void 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
174template <typename TInteger, typename TNumber>
175inline
176void 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
232template <typename TInteger, typename TNumber>
233inline
234void 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
326template <typename TInteger, typename TNumber>
327inline
328void 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
496template <typename TInteger, typename TNumber>
497inline
498void 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
601template <typename TInteger, typename TNumber>
602inline
603typename 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
626template <typename TInteger, typename TNumber>
627inline
628typename 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
635template <typename TInteger, typename TNumber>
636inline
637typename 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
643template <typename TInteger, typename TNumber>
644inline
645typename 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
658template <typename TInteger, typename TNumber>
659inline
660typename 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
677template <typename TInteger, typename TNumber>
678inline
679typename 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
696template <typename TInteger, typename TNumber>
697inline
698typename 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
718template <typename TInteger, typename TNumber>
719inline
720typename 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
797template <typename TInteger, typename TNumber>
798inline
799typename 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
874template <typename TInteger, typename TNumber>
875inline
876typename 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
903template <typename TInteger, typename TNumber>
904inline
905typename 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
966template <typename TInteger, typename TNumber>
967inline
968typename 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
984template <typename TInteger, typename TNumber>
985inline
986void 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
1036template <typename TInteger, typename TNumber>
1037inline
1038void 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
1165template <typename TInteger, typename TNumber>
1166inline
1167DGtal::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
1188template <typename TInteger, typename TNumber>
1189inline
1190void 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
1271template <typename TInteger, typename TNumber>
1272inline
1273void 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.
1347template <typename TInteger, typename TNumber>
1348inline
1349void 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
1441template <typename TInteger, typename TNumber>
1442inline
1443DGtal::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
1521template <typename TInteger, typename TNumber>
1522inline
1523TInteger
1524DGtal::DSLSubsegment<TInteger,TNumber>::getA() const {
1525 return myA;
1526}
1527
1528
1529template <typename TInteger, typename TNumber>
1530inline
1531TInteger
1532DGtal::DSLSubsegment<TInteger,TNumber>::getB() const {
1533 return myB;
1534}
1535
1536
1537
1538template <typename TInteger, typename TNumber>
1539inline
1540TInteger
1541DGtal::DSLSubsegment<TInteger,TNumber>::getMu() const {
1542 return myMu;
1543}
Generated on Sun Nov 27 2022 15:15:35 for DGtal by doxygen 1.9.5