src/2geom/bezier-clipping.cpp

   1 /*
   2  * Implement the Bezier clipping algorithm for finding
   3  * Bezier curve intersection points and collinear normals
   4  *
   5  * Authors:
   6  *      Marco Cecchetti <mrcekets at gmail.com>
   7  *
   8  * Copyright 2008  authors
   9  *
  10  * This library is free software; you can redistribute it and/or
  11  * modify it either under the terms of the GNU Lesser General Public
  12  * License version 2.1 as published by the Free Software Foundation
  13  * (the "LGPL") or, at your option, under the terms of the Mozilla
  14  * Public License Version 1.1 (the "MPL"). If you do not alter this
  15  * notice, a recipient may use your version of this file under either
  16  * the MPL or the LGPL.
  17  *
  18  * You should have received a copy of the LGPL along with this library
  19  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  21  * You should have received a copy of the MPL along with this library
  22  * in the file COPYING-MPL-1.1
  23  *
  24  * The contents of this file are subject to the Mozilla Public License
  25  * Version 1.1 (the "License"); you may not use this file except in
  26  * compliance with the License. You may obtain a copy of the License at
  27  * http://www.mozilla.org/MPL/
  28  *
  29  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  30  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  31  * the specific language governing rights and limitations.
  32  */
  33
  34
  35
  36
  37 #include <2geom/basic-intersection.h>
  38 #include <2geom/choose.h>
  39 #include <2geom/point.h>
  40 #include <2geom/interval.h>
  41 #include <2geom/bezier.h>
  42 //#include <2geom/convex-cover.h>
  43 #include <2geom/numeric/matrix.h>
  44
  45 #include <cassert>
  46 #include <vector>
  47 #include <algorithm>
  48 #include <utility>
  49 //#include <iomanip>
  50
  51
  52
  53
  54 #define VERBOSE 0
  55 #define CHECK 0
  56
  57
  58 namespace Geom {
  59
  60 namespace detail { namespace bezier_clipping {
  61
  62 ////////////////////////////////////////////////////////////////////////////////
  63 // for debugging
  64 //
  65
  66 inline
  67 void print(std::vector<Point> const& cp, const char* msg = "")
  68 {
  69     std::cerr << msg << std::endl;
  70     for (size_t i = 0; i < cp.size(); ++i)
  71         std::cerr << i << " : " << cp[i] << std::endl;
  72 }
  73
  74 template< class charT >
  75 inline
  76 std::basic_ostream<charT> &
  77 operator<< (std::basic_ostream<charT> & os, const Interval & I)
  78 {
  79     os << "[" << I.min() << ", " << I.max() << "]";
  80     return os;
  81 }
  82
  83 inline
  84 double angle (std::vector<Point> const& A)
  85 {
  86     size_t n = A.size() -1;
  87     double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);
  88     return (180 * a / M_PI);
  89 }
  90
  91 inline
  92 size_t get_precision(Interval const& I)
  93 {
  94     double d = I.extent();
  95     double e = 0.1, p = 10;
  96     int n = 0;
  97     while (n < 16 && d < e)
  98     {
  99         p *= 10;
 100         e = 1/p;
 101         ++n;
 102     }
 103     return n;
 104 }
 105
 106 inline
 107 void range_assertion(int k, int m, int n, const char* msg)
 108 {
 109     if ( k < m || k > n)
 110     {
 111         std::cerr << "range assertion failed: \n"
 112                   << msg << std::endl
 113                   << "value: " << k
 114                   << "  range: " << m << ", " << n << std::endl;
 115         assert (k >= m && k <= n);
 116     }
 117 }
 118
 119
 120 ////////////////////////////////////////////////////////////////////////////////
 121 //  convex hull
 122
 123 /*
 124  * return true in case the oriented polyline p0, p1, p2 is a right turn
 125  */
 126 inline
 127 bool is_a_right_turn (Point const& p0, Point const& p1, Point const& p2)
 128 {
 129     if (p1 == p2) return false;
 130     Point q1 = p1 - p0;
 131     Point q2 = p2 - p0;
 132     if (q1 == -q2) return false;
 133     return (cross (q1, q2) < 0);
 134 }
 135
 136 /*
 137  * return true if p < q wrt the lexicographyc order induced by the coordinates
 138  */
 139 struct lex_less
 140 {
 141     bool operator() (Point const& p, Point const& q)
 142     {
 143       return ((p[X] < q[X]) || (p[X] == q[X] && p[Y] < q[Y]));
 144     }
 145 };
 146
 147 /*
 148  * return true if p > q wrt the lexicographyc order induced by the coordinates
 149  */
 150 struct lex_greater
 151 {
 152     bool operator() (Point const& p, Point const& q)
 153     {
 154         return ((p[X] > q[X]) || (p[X] == q[X] && p[Y] > q[Y]));
 155     }
 156 };
 157
 158 /*
 159  * Compute the convex hull of a set of points.
 160  * The implementation is based on the Andrew's scan algorithm
 161  * note: in the Bezier clipping for collinear normals it seems
 162  * to be more stable wrt the Graham's scan algorithm and in general
 163  * a bit quikier
 164  */
 165 void convex_hull (std::vector<Point> & P)
 166 {
 167     size_t n = P.size();
 168     if (n < 2)  return;
 169     std::sort(P.begin(), P.end(), lex_less());
 170     if (n < 4) return;
 171     // upper hull
 172     size_t u = 2;
 173     for (size_t i = 2; i < n; ++i)
 174     {
 175         while (u > 1 && !is_a_right_turn(P[u-2], P[u-1], P[i]))
 176         {
 177             --u;
 178         }
 179         std::swap(P[u], P[i]);
 180         ++u;
 181     }
 182     std::sort(P.begin() + u, P.end(), lex_greater());
 183     std::rotate(P.begin(), P.begin() + 1, P.end());
 184     // lower hull
 185     size_t l = u;
 186     size_t k = u - 1;
 187     for (size_t i = l; i < n; ++i)
 188     {
 189         while (l > k && !is_a_right_turn(P[l-2], P[l-1], P[i]))
 190         {
 191             --l;
 192         }
 193         std::swap(P[l], P[i]);
 194         ++l;
 195     }
 196     P.resize(l);
 197 }
 198
 199
 200 ////////////////////////////////////////////////////////////////////////////////
 201 //  numerical routines
 202
 203 /*
 204  * Compute the binomial coefficient (n, k)
 205  */
 206 inline
 207 double binomial(unsigned int n, unsigned int k)
 208 {
 209     return choose<double>(n, k);
 210 }
 211
 212 /*
 213  * Compute the determinant of the 2x2 matrix with column the point P1, P2
 214  */
 215 inline
 216 double det(Point const& P1, Point const& P2)
 217 {
 218     return P1[X]*P2[Y] - P1[Y]*P2[X];
 219 }
 220
 221 /*
 222  * Solve the linear system [P1,P2] * P = Q
 223  * in case there isn't exactly one solution the routine returns false
 224  */
 225 inline
 226 bool solve(Point & P, Point const& P1, Point const& P2, Point const& Q)
 227 {
 228     double d = det(P1, P2);
 229     if (d == 0)  return false;
 230     d = 1 / d;
 231     P[X] = det(Q, P2) * d;
 232     P[Y] = det(P1, Q) * d;
 233     return true;
 234 }
 235
 236 ////////////////////////////////////////////////////////////////////////////////
 237 // interval routines
 238
 239 /*
 240  * Map the sub-interval I in [0,1] into the interval J and assign it to J
 241  */
 242 inline
 243 void map_to(Interval & J, Interval const& I)
 244 {
 245     double length = J.extent();
 246     J[1] = I.max() * length + J[0];
 247     J[0] = I.min() * length + J[0];
 248 }
 249
 250 /*
 251  * The interval [1,0] is used to represent the empty interval, this routine
 252  * is just an helper function for creating such an interval
 253  */
 254 inline
 255 Interval make_empty_interval()
 256 {
 257     Interval I(0);
 258     I[0] = 1;
 259     return I;
 260 }
 261
 262
 263 ////////////////////////////////////////////////////////////////////////////////
 264 // bezier curve routines
 265
 266 /*
 267  * Return true if all the Bezier curve control points are near,
 268  * false otherwise
 269  */
 270 inline
 271 bool is_constant(std::vector<Point> const& A, double precision = EPSILON)
 272 {
 273     for (unsigned int i = 1; i < A.size(); ++i)
 274     {
 275         if(!are_near(A[i], A[0], precision))
 276             return false;
 277     }
 278     return true;
 279 }
 280
 281 /*
 282  * Compute the hodograph of the bezier curve B and return it in D
 283  */
 284 inline
 285 void derivative(std::vector<Point> & D, std::vector<Point> const& B)
 286 {
 287     D.clear();
 288     size_t sz = B.size();
 289     if (sz == 0) return;
 290     if (sz == 1)
 291     {
 292         D.resize(1, Point(0,0));
 293         return;
 294     }
 295     size_t n = sz-1;
 296     D.reserve(n);
 297     for (size_t i = 0; i < n; ++i)
 298     {
 299         D.push_back(n*(B[i+1] - B[i]));
 300     }
 301 }
 302
 303 /*
 304  * Compute the hodograph of the Bezier curve B rotated of 90 degree
 305  * and return it in D; we have N(t) orthogonal to B(t) for any t
 306  */
 307 inline
 308 void normal(std::vector<Point> & N, std::vector<Point> const& B)
 309 {
 310     derivative(N,B);
 311     for (size_t i = 0; i < N.size(); ++i)
 312     {
 313         N[i] = rot90(N[i]);
 314     }
 315 }
 316
 317 /*
 318  *  Compute the portion of the Bezier curve "B" wrt the interval [0,t]
 319  */
 320 inline
 321 void left_portion(Coord t, std::vector<Point> & B)
 322 {
 323     size_t n = B.size();
 324     for (size_t i = 1; i < n; ++i)
 325     {
 326         for (size_t j = n-1; j > i-1 ; --j)
 327         {
 328             B[j] = lerp(t, B[j-1], B[j]);
 329         }
 330     }
 331 }
 332
 333 /*
 334  *  Compute the portion of the Bezier curve "B" wrt the interval [t,1]
 335  */
 336 inline
 337 void right_portion(Coord t, std::vector<Point> & B)
 338 {
 339     size_t n = B.size();
 340     for (size_t i = 1; i < n; ++i)
 341     {
 342         for (size_t j = 0; j < n-i; ++j)
 343         {
 344             B[j] = lerp(t, B[j], B[j+1]);
 345         }
 346     }
 347 }
 348
 349 /*
 350  *  Compute the portion of the Bezier curve "B" wrt the interval "I"
 351  */
 352 inline
 353 void portion (std::vector<Point> & B , Interval const& I)
 354 {
 355     if (I.min() == 0)
 356     {
 357         if (I.max() == 1)  return;
 358         left_portion(I.max(), B);
 359         return;
 360     }
 361     right_portion(I.min(), B);
 362     if (I.max() == 1)  return;
 363     double t = I.extent() / (1 - I.min());
 364     left_portion(t, B);
 365 }
 366
 367
 368 ////////////////////////////////////////////////////////////////////////////////
 369 // tags
 370
 371 struct intersection_point_tag;
 372 struct collinear_normal_tag;
 373 template <typename Tag>
 374 void clip(Interval & dom,
 375           std::vector<Point> const& A,
 376           std::vector<Point> const& B);
 377 template <typename Tag>
 378 void iterate(std::vector<Interval>& domsA,
 379              std::vector<Interval>& domsB,
 380              std::vector<Point> const& A,
 381              std::vector<Point> const& B,
 382              Interval const& domA,
 383              Interval const& domB,
 384              double precision );
 385
 386
 387 ////////////////////////////////////////////////////////////////////////////////
 388 // intersection
 389
 390 /*
 391  *  Make up an orientation line using the control points c[i] and c[j]
 392  *  the line is returned in the output parameter "l" in the form of a 3 element
 393  *  vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
 394  */
 395 inline
 396 void orientation_line (std::vector<double> & l,
 397                        std::vector<Point> const& c,
 398                        size_t i, size_t j)
 399 {
 400     l[0] = c[j][Y] - c[i][Y];
 401     l[1] = c[i][X] - c[j][X];
 402     l[2] = cross(c[i], c[j]);
 403     double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);
 404     assert (length != 0);
 405     l[0] /= length;
 406     l[1] /= length;
 407     l[2] /= length;
 408 }
 409
 410 /*
 411  * Pick up an orientation line for the Bezier curve "c" and return it in
 412  * the output parameter "l"
 413  */
 414 inline
 415 void pick_orientation_line (std::vector<double> & l,
 416                             std::vector<Point> const& c)
 417 {
 418     size_t i = c.size();
 419     while (--i > 0 && are_near(c[0], c[i]))
 420     {}
 421     if (i == 0)
 422     {
 423         // this should never happen because when a new curve portion is created
 424         // we check that it is not constant;
 425         // however this requires that the precision used in the is_constant
 426         // routine has to be the same used here in the are_near test
 427         assert(i != 0);
 428     }
 429     orientation_line(l, c, 0, i);
 430     //std::cerr << "i = " << i << std::endl;
 431 }
 432
 433 /*
 434  *  Make up an orientation line for constant bezier curve;
 435  *  the orientation line is made up orthogonal to the other curve base line;
 436  *  the line is returned in the output parameter "l" in the form of a 3 element
 437  *  vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
 438  */
 439 inline
 440 void orthogonal_orientation_line (std::vector<double> & l,
 441                                   std::vector<Point> const& c,
 442                                   Point const& p)
 443 {
 444     if (is_constant(c))
 445     {
 446         // this should never happen
 447         assert(!is_constant(c));
 448     }
 449     std::vector<Point> ol(2);
 450     ol[0] = p;
 451     ol[1] = (c.back() - c.front()).cw() + p;
 452     orientation_line(l, ol, 0, 1);
 453 }
 454
 455 /*
 456  *  Compute the signed distance of the point "P" from the normalized line l
 457  */
 458 inline
 459 double distance (Point const& P, std::vector<double> const& l)
 460 {
 461     return l[X] * P[X] + l[Y] * P[Y] + l[2];
 462 }
 463
 464 /*
 465  * Compute the min and max distance of the control points of the Bezier
 466  * curve "c" from the normalized orientation line "l".
 467  * This bounds are returned through the output Interval parameter"bound".
 468  */
 469 inline
 470 void fat_line_bounds (Interval& bound,
 471                       std::vector<Point> const& c,
 472                       std::vector<double> const& l)
 473 {
 474     bound[0] = 0;
 475     bound[1] = 0;
 476     double d;
 477     for (size_t i = 0; i < c.size(); ++i)
 478     {
 479         d = distance(c[i], l);
 480         if (bound[0] > d)  bound[0] = d;
 481         if (bound[1] < d)  bound[1] = d;
 482     }
 483 }
 484
 485 /*
 486  * return the x component of the intersection point between the line
 487  * passing through points p1, p2 and the line Y = "y"
 488  */
 489 inline
 490 double intersect (Point const& p1, Point const& p2, double y)
 491 {
 492     // we are sure that p2[Y] != p1[Y] because this routine is called
 493     // only when the lower or the upper bound is crossed
 494     double dy = (p2[Y] - p1[Y]);
 495     double s = (y - p1[Y]) / dy;
 496     return (p2[X]-p1[X])*s + p1[X];
 497 }
 498
 499 /*
 500  * Clip the Bezier curve "B" wrt the fat line defined by the orientation
 501  * line "l" and the interval range "bound", the new parameter interval for
 502  * the clipped curve is returned through the output parameter "dom"
 503  */
 504 void clip_interval (Interval& dom,
 505                     std::vector<Point> const& B,
 506                     std::vector<double> const& l,
 507                     Interval const& bound)
 508 {
 509     double n = B.size() - 1;  // number of sub-intervals
 510     std::vector<Point> D;     // distance curve control points
 511     D.reserve (B.size());
 512     double d;
 513     for (size_t i = 0; i < B.size(); ++i)
 514     {
 515         d = distance (B[i], l);
 516         D.push_back (Point(i/n, d));
 517     }
 518     //print(D);
 519
 520     convex_hull(D);
 521     std::vector<Point> & p = D;
 522     //print(p);
 523
 524     bool plower, phigher;
 525     bool clower, chigher;
 526     double t, tmin = 1, tmax = 0;
 527 //    std::cerr << "bound : " << bound << std::endl;
 528
 529     plower = (p[0][Y] < bound.min());
 530     phigher = (p[0][Y] > bound.max());
 531     if (!(plower || phigher))  // inside the fat line
 532     {
 533         if (tmin > p[0][X])  tmin = p[0][X];
 534         if (tmax < p[0][X])  tmax = p[0][X];
 535 //        std::cerr << "0 : inside " << p[0]
 536 //                  << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
 537     }
 538
 539     for (size_t i = 1; i < p.size(); ++i)
 540     {
 541         clower = (p[i][Y] < bound.min());
 542         chigher = (p[i][Y] > bound.max());
 543         if (!(clower || chigher))  // inside the fat line
 544         {
 545             if (tmin > p[i][X])  tmin = p[i][X];
 546             if (tmax < p[i][X])  tmax = p[i][X];
 547 //            std::cerr << i << " : inside " << p[i]
 548 //                      << " : tmin = " << tmin << ", tmax = " << tmax
 549 //                      << std::endl;
 550         }
 551         if (clower != plower)  // cross the lower bound
 552         {
 553             t = intersect(p[i-1], p[i], bound.min());
 554             if (tmin > t)  tmin = t;
 555             if (tmax < t)  tmax = t;
 556             plower = clower;
 557 //            std::cerr << i << " : lower " << p[i]
 558 //                      << " : tmin = " << tmin << ", tmax = " << tmax
 559 //                      << std::endl;
 560         }
 561         if (chigher != phigher)  // cross the upper bound
 562         {
 563             t = intersect(p[i-1], p[i], bound.max());
 564             if (tmin > t)  tmin = t;
 565             if (tmax < t)  tmax = t;
 566             phigher = chigher;
 567 //            std::cerr << i << " : higher " << p[i]
 568 //                      << " : tmin = " << tmin << ", tmax = " << tmax
 569 //                      << std::endl;
 570         }
 571     }
 572
 573     // we have to test the closing segment for intersection
 574     size_t last = p.size() - 1;
 575     clower = (p[0][Y] < bound.min());
 576     chigher = (p[0][Y] > bound.max());
 577     if (clower != plower)  // cross the lower bound
 578     {
 579         t = intersect(p[last], p[0], bound.min());
 580         if (tmin > t)  tmin = t;
 581         if (tmax < t)  tmax = t;
 582 //        std::cerr << "0 : lower " << p[0]
 583 //                  << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
 584     }
 585     if (chigher != phigher)  // cross the upper bound
 586     {
 587         t = intersect(p[last], p[0], bound.max());
 588         if (tmin > t)  tmin = t;
 589         if (tmax < t)  tmax = t;
 590 //        std::cerr << "0 : higher " << p[0]
 591 //                  << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
 592     }
 593
 594     dom[0] = tmin;
 595     dom[1] = tmax;
 596 }
 597
 598 /*
 599  *  Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
 600  *  intersection points the new parameter interval for the clipped curve
 601  *  is returned through the output parameter "dom"
 602  */
 603 template <>
 604 inline
 605 void clip<intersection_point_tag> (Interval & dom,
 606                                    std::vector<Point> const& A,
 607                                    std::vector<Point> const& B)
 608 {
 609     std::vector<double> bl(3);
 610     Interval bound;
 611     if (is_constant(A))
 612     {
 613         Point M = middle_point(A.front(), A.back());
 614         orthogonal_orientation_line(bl, B, M);
 615     }
 616     else
 617     {
 618         pick_orientation_line(bl, A);
 619     }
 620     fat_line_bounds(bound, A, bl);
 621     clip_interval(dom, B, bl, bound);
 622 }
 623
 624
 625 ///////////////////////////////////////////////////////////////////////////////
 626 // collinear normal
 627
 628 /*
 629  * Compute a closed focus for the Bezier curve B and return it in F
 630  * A focus is any curve through which all lines perpendicular to B(t) pass.
 631  */
 632 inline
 633 void make_focus (std::vector<Point> & F, std::vector<Point> const& B)
 634 {
 635     assert (B.size() > 2);
 636     size_t n = B.size() - 1;
 637     normal(F, B);
 638     Point c(1, 1);
 639 #if VERBOSE
 640     if (!solve(c, F[0], -F[n-1], B[n]-B[0]))
 641     {
 642         std::cerr << "make_focus: unable to make up a closed focus" << std::endl;
 643     }
 644 #else
 645     solve(c, F[0], -F[n-1], B[n]-B[0]);
 646 #endif
 647 //    std::cerr << "c = " << c << std::endl;
 648
 649
 650     // B(t) + c(t) * N(t)
 651     double n_inv = 1 / (double)(n);
 652     Point c0ni;
 653     F.push_back(c[1] * F[n-1]);
 654     F[n] += B[n];
 655     for (size_t i = n-1; i > 0; --i)
 656     {
 657         F[i] *= -c[0];
 658         c0ni = F[i];
 659         F[i] += (c[1] * F[i-1]);
 660         F[i] *= (i * n_inv);
 661         F[i] -= c0ni;
 662         F[i] += B[i];
 663     }
 664     F[0] *= c[0];
 665     F[0] += B[0];
 666 }
 667
 668 /*
 669  * Compute the projection on the plane (t, d) of the control points
 670  * (t, u, D(t,u)) where D(t,u) = <(B(t) - F(u)), B'(t)> with 0 <= t, u <= 1
 671  * B is a Bezier curve and F is a focus of another Bezier curve.
 672  * See Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping.
 673  */
 674 void distance_control_points (std::vector<Point> & D,
 675                               std::vector<Point> const& B,
 676                               std::vector<Point> const& F)
 677 {
 678     assert (B.size() > 1);
 679     assert (F.size() > 0);
 680     const size_t n = B.size() - 1;
 681     const size_t m = F.size() - 1;
 682     const size_t r = 2 * n - 1;
 683     const double r_inv = 1 / (double)(r);
 684     D.clear();
 685     D.reserve (B.size() * F.size());
 686
 687     std::vector<Point> dB;
 688     dB.reserve(n);
 689     for (size_t k = 0; k < n; ++k)
 690     {
 691         dB.push_back (B[k+1] - B[k]);
 692     }
 693     NL::Matrix dBB(n,B.size());
 694     for (size_t i = 0; i < n; ++i)
 695         for (size_t j = 0; j < B.size(); ++j)
 696             dBB(i,j) = dot (dB[i], B[j]);
 697     NL::Matrix dBF(n, F.size());
 698     for (size_t i = 0; i < n; ++i)
 699         for (size_t j = 0; j < F.size(); ++j)
 700             dBF(i,j) = dot (dB[i], F[j]);
 701
 702     size_t k0, kn, l;
 703     double bc, bri;
 704     Point dij;
 705     std::vector<double> d(F.size());
 706     for (size_t i = 0; i <= r; ++i)
 707     {
 708         for (size_t j = 0; j <= m; ++j)
 709         {
 710             d[j] = 0;
 711         }
 712         k0 = std::max(i, n) - n;
 713         kn = std::min(i, n-1);
 714         bri = n / binomial(r,i);
 715         for (size_t k = k0; k <= kn; ++k)
 716         {
 717             //if (k > i || (i-k) > n) continue;
 718             l = i - k;
 719 #if CHECK
 720             assert (l <= n);
 721 #endif
 722             bc = bri * binomial(n,l) * binomial(n-1, k);
 723             for (size_t j = 0; j <= m; ++j)
 724             {
 725                 //d[j] += bc * dot(dB[k], B[l] - F[j]);
 726                 d[j] += bc * (dBB(k,l) - dBF(k,j));
 727             }
 728         }
 729         double dmin, dmax;
 730         dmin = dmax = d[m];
 731         for (size_t j = 0; j < m; ++j)
 732         {
 733             if (dmin > d[j])  dmin = d[j];
 734             if (dmax < d[j])  dmax = d[j];
 735         }
 736         dij[0] = i * r_inv;
 737         dij[1] = dmin;
 738         D.push_back (dij);
 739         dij[1] = dmax;
 740         D.push_back (dij);
 741     }
 742 }
 743
 744 /*
 745  * Clip the Bezier curve "B" wrt the focus "F"; the new parameter interval for
 746  * the clipped curve is returned through the output parameter "dom"
 747  */
 748 void clip_interval (Interval& dom,
 749                     std::vector<Point> const& B,
 750                     std::vector<Point> const& F)
 751 {
 752     std::vector<Point> D;     // distance curve control points
 753     distance_control_points(D, B, F);
 754     //print(D, "D");
 755 //    ConvexHull chD(D);
 756 //    std::vector<Point>& p = chD.boundary; // convex hull vertices
 757
 758     convex_hull(D);
 759     std::vector<Point> & p = D;
 760     //print(p, "CH(D)");
 761
 762     bool plower, clower;
 763     double t, tmin = 1, tmax = 0;
 764
 765     plower = (p[0][Y] < 0);
 766     if (p[0][Y] == 0)  // on the x axis
 767     {
 768         if (tmin > p[0][X])  tmin = p[0][X];
 769         if (tmax < p[0][X])  tmax = p[0][X];
 770 //        std::cerr << "0 : on x axis " << p[0]
 771 //                  << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
 772     }
 773
 774     for (size_t i = 1; i < p.size(); ++i)
 775     {
 776         clower = (p[i][Y] < 0);
 777         if (p[i][Y] == 0)  // on x axis
 778         {
 779             if (tmin > p[i][X])  tmin = p[i][X];
 780             if (tmax < p[i][X])  tmax = p[i][X];
 781 //            std::cerr << i << " : on x axis " << p[i]
 782 //                      << " : tmin = " << tmin << ", tmax = " << tmax
 783 //                      << std::endl;
 784         }
 785         else if (clower != plower)  // cross the x axis
 786         {
 787             t = intersect(p[i-1], p[i], 0);
 788             if (tmin > t)  tmin = t;
 789             if (tmax < t)  tmax = t;
 790             plower = clower;
 791 //            std::cerr << i << " : lower " << p[i]
 792 //                      << " : tmin = " << tmin << ", tmax = " << tmax
 793 //                      << std::endl;
 794         }
 795     }
 796
 797     // we have to test the closing segment for intersection
 798     size_t last = p.size() - 1;
 799     clower = (p[0][Y] < 0);
 800     if (clower != plower)  // cross the x axis
 801     {
 802         t = intersect(p[last], p[0], 0);
 803         if (tmin > t)  tmin = t;
 804         if (tmax < t)  tmax = t;
 805 //        std::cerr << "0 : lower " << p[0]
 806 //                  << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
 807     }
 808     dom[0] = tmin;
 809     dom[1] = tmax;
 810 }
 811
 812 /*
 813  *  Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
 814  *  points which have collinear normals; the new parameter interval
 815  *  for the clipped curve is returned through the output parameter "dom"
 816  */
 817 template <>
 818 inline
 819 void clip<collinear_normal_tag> (Interval & dom,
 820                                  std::vector<Point> const& A,
 821                                  std::vector<Point> const& B)
 822 {
 823     std::vector<Point> F;
 824     make_focus(F, A);
 825     clip_interval(dom, B, F);
 826 }
 827
 828
 829
 830 const double MAX_PRECISION = 1e-8;
 831 const double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;
 832 const Interval UNIT_INTERVAL(0,1);
 833 const Interval EMPTY_INTERVAL = make_empty_interval();
 834 const Interval H1_INTERVAL(0, 0.5);
 835 const Interval H2_INTERVAL(0.5 + MAX_PRECISION, 1.0);
 836
 837 /*
 838  * iterate
 839  *
 840  * input:
 841  * A, B: control point sets of two bezier curves
 842  * domA, domB: real parameter intervals of the two curves
 843  * precision: required computational precision of the returned parameter ranges
 844  * output:
 845  * domsA, domsB: sets of parameter intervals
 846  *
 847  * The parameter intervals are computed by using a Bezier clipping algorithm,
 848  * in case the clipping doesn't shrink the initial interval more than 20%,
 849  * a subdivision step is performed.
 850  * If during the computation both curves collapse to a single point
 851  * the routine exits indipendently by the precision reached in the computation
 852  * of the curve intervals.
 853  */
 854 template <>
 855 void iterate<intersection_point_tag> (std::vector<Interval>& domsA,
 856                                       std::vector<Interval>& domsB,
 857                                       std::vector<Point> const& A,
 858                                       std::vector<Point> const& B,
 859                                       Interval const& domA,
 860                                       Interval const& domB,
 861                                       double precision )
 862 {
 863     // in order to limit recursion
 864     static size_t counter = 0;
 865     if (domA.extent() == 1 && domB.extent() == 1) counter  = 0;
 866     if (++counter > 100) return;
 867 #if VERBOSE
 868     std::cerr << std::fixed << std::setprecision(16);
 869     std::cerr << ">> curve subdision performed <<" << std::endl;
 870     std::cerr << "dom(A) : " << domA << std::endl;
 871     std::cerr << "dom(B) : " << domB << std::endl;
 872 //    std::cerr << "angle(A) : " << angle(A) << std::endl;
 873 //    std::cerr << "angle(B) : " << angle(B) << std::endl;
 874 #endif
 875
 876     if (precision < MAX_PRECISION)
 877         precision = MAX_PRECISION;
 878
 879     std::vector<Point> pA = A;
 880     std::vector<Point> pB = B;
 881     std::vector<Point>* C1 = &pA;
 882     std::vector<Point>* C2 = &pB;
 883
 884     Interval dompA = domA;
 885     Interval dompB = domB;
 886     Interval* dom1 = &dompA;
 887     Interval* dom2 = &dompB;
 888
 889     Interval dom;
 890
 891     if ( is_constant(A) && is_constant(B) ){
 892         Point M1 = middle_point(C1->front(), C1->back());
 893         Point M2 = middle_point(C2->front(), C2->back());
 894         if (are_near(M1,M2)){
 895             domsA.push_back(domA);
 896             domsB.push_back(domB);
 897         }
 898         return;
 899     }
 900
 901     size_t iter = 0;
 902     while (++iter < 100
 903             && (dompA.extent() >= precision || dompB.extent() >= precision))
 904     {
 905 #if VERBOSE
 906         std::cerr << "iter: " << iter << std::endl;
 907 #endif
 908         clip<intersection_point_tag>(dom, *C1, *C2);
 909
 910         // [1,0] is utilized to represent an empty interval
 911         if (dom == EMPTY_INTERVAL)
 912         {
 913 #if VERBOSE
 914             std::cerr << "dom: empty" << std::endl;
 915 #endif
 916             return;
 917         }
 918 #if VERBOSE
 919         std::cerr << "dom : " << dom << std::endl;
 920 #endif
 921         // all other cases where dom[0] > dom[1] are invalid
 922         if (dom.min() >  dom.max())
 923         {
 924             assert(dom.min() <  dom.max());
 925         }
 926
 927         map_to(*dom2, dom);
 928
 929         portion(*C2, dom);
 930         if (is_constant(*C2) && is_constant(*C1))
 931         {
 932             Point M1 = middle_point(C1->front(), C1->back());
 933             Point M2 = middle_point(C2->front(), C2->back());
 934 #if VERBOSE
 935             std::cerr << "both curves are constant: \n"
 936                       << "M1: " << M1 << "\n"
 937                       << "M2: " << M2 << std::endl;
 938             print(*C2, "C2");
 939             print(*C1, "C1");
 940 #endif
 941             if (are_near(M1,M2))
 942                 break;  // append the new interval
 943             else
 944                 return; // exit without appending any new interval
 945         }
 946
 947
 948         // if we have clipped less than 20% than we need to subdive the curve
 949         // with the largest domain into two sub-curves
 950         if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
 951         {
 952 #if VERBOSE
 953             std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
 954             std::cerr << "angle(pA) : " << angle(pA) << std::endl;
 955             std::cerr << "angle(pB) : " << angle(pB) << std::endl;
 956 #endif
 957             std::vector<Point> pC1, pC2;
 958             Interval dompC1, dompC2;
 959             if (dompA.extent() > dompB.extent())
 960             {
 961                 pC1 = pC2 = pA;
 962                 portion(pC1, H1_INTERVAL);
 963                 portion(pC2, H2_INTERVAL);
 964                 dompC1 = dompC2 = dompA;
 965                 map_to(dompC1, H1_INTERVAL);
 966                 map_to(dompC2, H2_INTERVAL);
 967                 iterate<intersection_point_tag>(domsA, domsB, pC1, pB,
 968                                                 dompC1, dompB, precision);
 969                 iterate<intersection_point_tag>(domsA, domsB, pC2, pB,
 970                                                 dompC2, dompB, precision);
 971             }
 972             else
 973             {
 974                 pC1 = pC2 = pB;
 975                 portion(pC1, H1_INTERVAL);
 976                 portion(pC2, H2_INTERVAL);
 977                 dompC1 = dompC2 = dompB;
 978                 map_to(dompC1, H1_INTERVAL);
 979                 map_to(dompC2, H2_INTERVAL);
 980                 iterate<intersection_point_tag>(domsB, domsA, pC1, pA,
 981                                                 dompC1, dompA, precision);
 982                 iterate<intersection_point_tag>(domsB, domsA, pC2, pA,
 983                                                 dompC2, dompA, precision);
 984             }
 985             return;
 986         }
 987
 988         std::swap(C1, C2);
 989         std::swap(dom1, dom2);
 990 #if VERBOSE
 991         std::cerr << "dom(pA) : " << dompA << std::endl;
 992         std::cerr << "dom(pB) : " << dompB << std::endl;
 993 #endif
 994     }
 995     domsA.push_back(dompA);
 996     domsB.push_back(dompB);
 997 }
 998
 999
1000 /*
1001  * iterate
1002  *
1003  * input:
1004  * A, B: control point sets of two bezier curves
1005  * domA, domB: real parameter intervals of the two curves
1006  * precision: required computational precision of the returned parameter ranges
1007  * output:
1008  * domsA, domsB: sets of parameter intervals
1009  *
1010  * The parameter intervals are computed by using a Bezier clipping algorithm,
1011  * in case the clipping doesn't shrink the initial interval more than 20%,
1012  * a subdivision step is performed.
1013  * If during the computation one of the two curve interval length becomes less
1014  * than MAX_PRECISION the routine exits indipendently by the precision reached
1015  * in the computation of the other curve interval.
1016  */
1017 template <>
1018 void iterate<collinear_normal_tag> (std::vector<Interval>& domsA,
1019                                     std::vector<Interval>& domsB,
1020                                     std::vector<Point> const& A,
1021                                     std::vector<Point> const& B,
1022                                     Interval const& domA,
1023                                     Interval const& domB,
1024                                     double precision)
1025 {
1026     // in order to limit recursion
1027     static size_t counter = 0;
1028     if (domA.extent() == 1 && domB.extent() == 1) counter  = 0;
1029     if (++counter > 100) return;
1030 #if VERBOSE
1031     std::cerr << std::fixed << std::setprecision(16);
1032     std::cerr << ">> curve subdision performed <<" << std::endl;
1033     std::cerr << "dom(A) : " << domA << std::endl;
1034     std::cerr << "dom(B) : " << domB << std::endl;
1035 //    std::cerr << "angle(A) : " << angle(A) << std::endl;
1036 //    std::cerr << "angle(B) : " << angle(B) << std::endl;
1037 #endif
1038
1039     if (precision < MAX_PRECISION)
1040         precision = MAX_PRECISION;
1041
1042     std::vector<Point> pA = A;
1043     std::vector<Point> pB = B;
1044     std::vector<Point>* C1 = &pA;
1045     std::vector<Point>* C2 = &pB;
1046
1047     Interval dompA = domA;
1048     Interval dompB = domB;
1049     Interval* dom1 = &dompA;
1050     Interval* dom2 = &dompB;
1051
1052     Interval dom;
1053
1054     size_t iter = 0;
1055     while (++iter < 100
1056             && (dompA.extent() >= precision || dompB.extent() >= precision))
1057     {
1058 #if VERBOSE
1059         std::cerr << "iter: " << iter << std::endl;
1060 #endif
1061         clip<collinear_normal_tag>(dom, *C1, *C2);
1062
1063         // [1,0] is utilized to represent an empty interval
1064         if (dom == EMPTY_INTERVAL)
1065         {
1066 #if VERBOSE
1067             std::cerr << "dom: empty" << std::endl;
1068 #endif
1069             return;
1070         }
1071 #if VERBOSE
1072         std::cerr << "dom : " << dom << std::endl;
1073 #endif
1074         // all other cases where dom[0] > dom[1] are invalid
1075         if (dom.min() >  dom.max())
1076         {
1077             assert(dom.min() <  dom.max());
1078         }
1079
1080         map_to(*dom2, dom);
1081
1082         // it's better to stop before losing computational precision
1083         if (iter > 1 && (dom2->extent() <= MAX_PRECISION))
1084         {
1085 #if VERBOSE
1086             std::cerr << "beyond max precision limit" << std::endl;
1087 #endif
1088             break;
1089         }
1090
1091         portion(*C2, dom);
1092         if (iter > 1 && is_constant(*C2))
1093         {
1094 #if VERBOSE
1095             std::cerr << "new curve portion pC1 is constant" << std::endl;
1096 #endif
1097             break;
1098         }
1099
1100
1101         // if we have clipped less than 20% than we need to subdive the curve
1102         // with the largest domain into two sub-curves
1103         if ( dom.extent() > MIN_CLIPPED_SIZE_THRESHOLD)
1104         {
1105 #if VERBOSE
1106             std::cerr << "clipped less than 20% : " << dom.extent() << std::endl;
1107             std::cerr << "angle(pA) : " << angle(pA) << std::endl;
1108             std::cerr << "angle(pB) : " << angle(pB) << std::endl;
1109 #endif
1110             std::vector<Point> pC1, pC2;
1111             Interval dompC1, dompC2;
1112             if (dompA.extent() > dompB.extent())
1113             {
1114                 if ((dompA.extent() / 2) < MAX_PRECISION)
1115                 {
1116                     break;
1117                 }
1118                 pC1 = pC2 = pA;
1119                 portion(pC1, H1_INTERVAL);
1120                 if (false && is_constant(pC1))
1121                 {
1122 #if VERBOSE
1123                     std::cerr << "new curve portion pC1 is constant" << std::endl;
1124 #endif
1125                     break;
1126                 }
1127                 portion(pC2, H2_INTERVAL);
1128                 if (is_constant(pC2))
1129                 {
1130 #if VERBOSE
1131                     std::cerr << "new curve portion pC2 is constant" << std::endl;
1132 #endif
1133                     break;
1134                 }
1135                 dompC1 = dompC2 = dompA;
1136                 map_to(dompC1, H1_INTERVAL);
1137                 map_to(dompC2, H2_INTERVAL);
1138                 iterate<collinear_normal_tag>(domsA, domsB, pC1, pB,
1139                                               dompC1, dompB, precision);
1140                 iterate<collinear_normal_tag>(domsA, domsB, pC2, pB,
1141                                               dompC2, dompB, precision);
1142             }
1143             else
1144             {
1145                 if ((dompB.extent() / 2) < MAX_PRECISION)
1146                 {
1147                     break;
1148                 }
1149                 pC1 = pC2 = pB;
1150                 portion(pC1, H1_INTERVAL);
1151                 if (is_constant(pC1))
1152                 {
1153 #if VERBOSE
1154                     std::cerr << "new curve portion pC1 is constant" << std::endl;
1155 #endif
1156                     break;
1157                 }
1158                 portion(pC2, H2_INTERVAL);
1159                 if (is_constant(pC2))
1160                 {
1161 #if VERBOSE
1162                     std::cerr << "new curve portion pC2 is constant" << std::endl;
1163 #endif
1164                     break;
1165                 }
1166                 dompC1 = dompC2 = dompB;
1167                 map_to(dompC1, H1_INTERVAL);
1168                 map_to(dompC2, H2_INTERVAL);
1169                 iterate<collinear_normal_tag>(domsB, domsA, pC1, pA,
1170                                               dompC1, dompA, precision);
1171                 iterate<collinear_normal_tag>(domsB, domsA, pC2, pA,
1172                                               dompC2, dompA, precision);
1173             }
1174             return;
1175         }
1176
1177         std::swap(C1, C2);
1178         std::swap(dom1, dom2);
1179 #if VERBOSE
1180         std::cerr << "dom(pA) : " << dompA << std::endl;
1181         std::cerr << "dom(pB) : " << dompB << std::endl;
1182 #endif
1183     }
1184     domsA.push_back(dompA);
1185     domsB.push_back(dompB);
1186 }
1187
1188
1189 /*
1190  * get_solutions
1191  *
1192  *  input: A, B       - set of control points of two Bezier curve
1193  *  input: precision  - required precision of computation
1194  *  input: clip       - the routine used for clipping
1195  *  output: xs        - set of pairs of parameter values
1196  *                      at which the clipping algorithm converges
1197  *
1198  *  This routine is based on the Bezier Clipping Algorithm,
1199  *  see: Sederberg - Computer Aided Geometric Design
1200  */
1201 template <typename Tag>
1202 void get_solutions (std::vector< std::pair<double, double> >& xs,
1203                     std::vector<Point> const& A,
1204                     std::vector<Point> const& B,
1205                     double precision)
1206 {
1207     std::pair<double, double> ci;
1208     std::vector<Interval> domsA, domsB;
1209     iterate<Tag> (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);
1210     if (domsA.size() != domsB.size())
1211     {
1212         assert (domsA.size() == domsB.size());
1213     }
1214     xs.clear();
1215     xs.reserve(domsA.size());
1216     for (size_t i = 0; i < domsA.size(); ++i)
1217     {
1218 #if VERBOSE
1219         std::cerr << i << " : domA : " << domsA[i] << std::endl;
1220         std::cerr << "extent A: " << domsA[i].extent() << "  ";
1221         std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;
1222         std::cerr << i << " : domB : " << domsB[i] << std::endl;
1223         std::cerr << "extent B: " << domsB[i].extent() << "  ";
1224         std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;
1225 #endif
1226         ci.first = domsA[i].middle();
1227         ci.second = domsB[i].middle();
1228         xs.push_back(ci);
1229     }
1230 }
1231
1232 } /* end namespace bezier_clipping */ } /* end namespace detail */
1233
1234
1235 /*
1236  * find_collinear_normal
1237  *
1238  *  input: A, B       - set of control points of two Bezier curve
1239  *  input: precision  - required precision of computation
1240  *  output: xs        - set of pairs of parameter values
1241  *                      at which there are collinear normals
1242  *
1243  *  This routine is based on the Bezier Clipping Algorithm,
1244  *  see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
1245  */
1246 void find_collinear_normal (std::vector< std::pair<double, double> >& xs,
1247                             std::vector<Point> const& A,
1248                             std::vector<Point> const& B,
1249                             double precision)
1250 {
1251     using detail::bezier_clipping::get_solutions;
1252     using detail::bezier_clipping::collinear_normal_tag;
1253     get_solutions<collinear_normal_tag>(xs, A, B, precision);
1254 }
1255
1256
1257 /*
1258  * find_intersections_bezier_clipping
1259  *
1260  *  input: A, B       - set of control points of two Bezier curve
1261  *  input: precision  - required precision of computation
1262  *  output: xs        - set of pairs of parameter values
1263  *                      at which crossing happens
1264  *
1265  *  This routine is based on the Bezier Clipping Algorithm,
1266  *  see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
1267  */
1268 void find_intersections_bezier_clipping (std::vector< std::pair<double, double> >& xs,
1269                          std::vector<Point> const& A,
1270                          std::vector<Point> const& B,
1271                          double precision)
1272 {
1273     using detail::bezier_clipping::get_solutions;
1274     using detail::bezier_clipping::intersection_point_tag;
1275     get_solutions<intersection_point_tag>(xs, A, B, precision);
1276 }
1277
1278 }  // end namespace Geom
1279
1280
1281
1282
1283 /*
1284   Local Variables:
1285   mode:c++
1286   c-file-style:"stroustrup"
1287   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1288   indent-tabs-mode:nil
1289   fill-column:99
1290   End:
1291 */
1292 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :