src/2geom/basic-intersection.cpp

   1 #include "basic-intersection.h"
   2
   3 unsigned intersect_steps = 0;
   4
   5 using std::vector;
   6 namespace Geom {
   7
   8 class OldBezier {
   9 public:
  10     std::vector<Geom::Point> p;
  11     OldBezier() {
  12     }
  13     void split(double t, OldBezier &a, OldBezier &b) const;
  14
  15     ~OldBezier() {}
  16
  17     void bounds(double &minax, double &maxax,
  18                 double &minay, double &maxay) {
  19         // Compute bounding box for a
  20         minax = p[0][X];         // These are the most likely to be extremal
  21         maxax = p.back()[X];
  22         if( minax > maxax )
  23             std::swap(minax, maxax);
  24         for(unsigned i = 1; i < p.size()-1; i++) {
  25             if( p[i][X] < minax )
  26                 minax = p[i][X];
  27             else if( p[i][X] > maxax )
  28                 maxax = p[i][X];
  29         }
  30
  31         minay = p[0][Y];         // These are the most likely to be extremal
  32         maxay = p.back()[Y];
  33         if( minay > maxay )
  34             std::swap(minay, maxay);
  35         for(unsigned i = 1; i < p.size()-1; i++) {
  36             if( p[i][Y] < minay )
  37                 minay = p[i][Y];
  38             else if( p[i][Y] > maxay )
  39                 maxay = p[i][Y];
  40         }
  41
  42     }
  43
  44 };
  45
  46 static std::vector<std::pair<double, double> >
  47 find_intersections( OldBezier a, OldBezier b);
  48
  49 static std::vector<std::pair<double, double> >
  50 find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A);
  51
  52 std::vector<std::pair<double, double> >
  53 find_intersections( vector<Geom::Point> const & A,
  54                     vector<Geom::Point> const & B) {
  55     OldBezier a, b;
  56     a.p = A;
  57     b.p = B;
  58     return find_intersections(a,b);
  59 }
  60
  61 std::vector<std::pair<double, double> >
  62 find_self_intersections(OldBezier const &Sb) {
  63     throwNotImplemented();
  64 }
  65
  66 std::vector<std::pair<double, double> >
  67 find_self_intersections(D2<SBasis> const & A) {
  68     OldBezier Sb;
  69     Sb.p = sbasis_to_bezier(A);
  70     return find_self_intersections(Sb, A);
  71 }
  72
  73
  74 static std::vector<std::pair<double, double> >
  75 find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
  76
  77
  78     vector<double> dr = roots(derivative(A[X]));
  79     {
  80         vector<double> dyr = roots(derivative(A[Y]));
  81         dr.insert(dr.begin(), dyr.begin(), dyr.end());
  82     }
  83     dr.push_back(0);
  84     dr.push_back(1);
  85     // We want to be sure that we have no empty segments
  86     sort(dr.begin(), dr.end());
  87     unique(dr.begin(), dr.end());
  88
  89     std::vector<std::pair<double, double> > all_si;
  90
  91     vector<OldBezier> pieces;
  92     {
  93         OldBezier in = Sb, l, r;
  94         for(unsigned i = 0; i < dr.size()-1; i++) {
  95             in.split((dr[i+1]-dr[i]) / (1 - dr[i]), l, r);
  96             pieces.push_back(l);
  97             in = r;
  98         }
  99     }
 100     for(unsigned i = 0; i < dr.size()-1; i++) {
 101         for(unsigned j = i+1; j < dr.size()-1; j++) {
 102             std::vector<std::pair<double, double> > section =
 103                 find_intersections( pieces[i], pieces[j]);
 104             for(unsigned k = 0; k < section.size(); k++) {
 105                 double l = section[k].first;
 106                 double r = section[k].second;
 107 // XXX: This condition will prune out false positives, but it might create some false negatives.  Todo: Confirm it is correct.
 108                 if(j == i+1)
 109                     if((l == 1) && (r == 0))
 110                         continue;
 111                 all_si.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1],
 112                                                 (1-r)*dr[j] + r*dr[j+1]));
 113             }
 114         }
 115     }
 116
 117     // Because i is in order, all_si should be roughly already in order?
 118     //sort(all_si.begin(), all_si.end());
 119     //unique(all_si.begin(), all_si.end());
 120
 121     return all_si;
 122 }
 123
 124 /* The value of 1.0 / (1L<<14) is enough for most applications */
 125 const double INV_EPS = (1L<<14);
 126
 127 /*
 128  * split the curve at the midpoint, returning an array with the two parts
 129  * Temporary storage is minimized by using part of the storage for the result
 130  * to hold an intermediate value until it is no longer needed.
 131  */
 132 void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
 133     const unsigned sz = p.size();
 134     Geom::Point Vtemp[sz][sz];
 135
 136     /* Copy control points      */
 137     std::copy(p.begin(), p.end(), Vtemp[0]);
 138
 139     /* Triangle computation     */
 140     for (unsigned i = 1; i < sz; i++) {
 141         for (unsigned j = 0; j < sz - i; j++) {
 142             Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
 143         }
 144     }
 145
 146     left.p.resize(sz);
 147     right.p.resize(sz);
 148     for (unsigned j = 0; j < sz; j++)
 149         left.p[j]  = Vtemp[j][0];
 150     for (unsigned j = 0; j < sz; j++)
 151         right.p[j] = Vtemp[sz-1-j][j];
 152 }
 153
 154
 155 /*
 156  * Test the bounding boxes of two OldBezier curves for interference.
 157  * Several observations:
 158  *      First, it is cheaper to compute the bounding box of the second curve
 159  *      and test its bounding box for interference than to use a more direct
 160  *      approach of comparing all control points of the second curve with
 161  *      the various edges of the bounding box of the first curve to test
 162  *      for interference.
 163  *      Second, after a few subdivisions it is highly probable that two corners
 164  *      of the bounding box of a given Bezier curve are the first and last
 165  *      control point.  Once this happens once, it happens for all subsequent
 166  *      subcurves.  It might be worth putting in a test and then short-circuit
 167  *      code for further subdivision levels.
 168  *      Third, in the final comparison (the interference test) the comparisons
 169  *      should both permit equality.  We want to find intersections even if they
 170  *      occur at the ends of segments.
 171  *      Finally, there are tighter bounding boxes that can be derived. It isn't
 172  *      clear whether the higher probability of rejection (and hence fewer
 173  *      subdivisions and tests) is worth the extra work.
 174  */
 175
 176 bool intersect_BB( OldBezier a, OldBezier b ) {
 177     double minax, maxax, minay, maxay;
 178     a.bounds(minax, maxax, minay, maxay);
 179     double minbx, maxbx, minby, maxby;
 180     b.bounds(minbx, maxbx, minby, maxby);
 181     // Test bounding box of b against bounding box of a
 182     // Not >= : need boundary case
 183     return not( ( minax > maxbx ) || ( minay > maxby )
 184                 || ( minbx > maxax ) || ( minby > maxay ) );
 185 }
 186
 187 /*
 188  * Recursively intersect two curves keeping track of their real parameters
 189  * and depths of intersection.
 190  * The results are returned in a 2-D array of doubles indicating the parameters
 191  * for which intersections are found.  The parameters are in the order the
 192  * intersections were found, which is probably not in sorted order.
 193  * When an intersection is found, the parameter value for each of the two
 194  * is stored in the index elements array, and the index is incremented.
 195  *
 196  * If either of the curves has subdivisions left before it is straight
 197  *      (depth > 0)
 198  * that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
 199  * the depth(s) is (are) decremented, and the parameter value(s) corresponding
 200  * to the midpoints(s) is (are) computed.
 201  * Then each of the subcurves of one curve is intersected with each of the
 202  * subcurves of the other curve, first by testing the bounding boxes for
 203  * interference.  If there is any bounding box interference, the corresponding
 204  * subcurves are recursively intersected.
 205  *
 206  * If neither curve has subdivisions left, the line segments from the first
 207  * to last control point of each segment are intersected.  (Actually the
 208  * only the parameter value corresponding to the intersection point is found).
 209  *
 210  * The apriori flatness test is probably more efficient than testing at each
 211  * level of recursion, although a test after three or four levels would
 212  * probably be worthwhile, since many curves become flat faster than their
 213  * asymptotic rate for the first few levels of recursion.
 214  *
 215  * The bounding box test fails much more frequently than it succeeds, providing
 216  * substantial pruning of the search space.
 217  *
 218  * Each (sub)curve is subdivided only once, hence it is not possible that for
 219  * one final line intersection test the subdivision was at one level, while
 220  * for another final line intersection test the subdivision (of the same curve)
 221  * was at another.  Since the line segments share endpoints, the intersection
 222  * is robust: a near-tangential intersection will yield zero or two
 223  * intersections.
 224  */
 225 void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
 226                            OldBezier b, double u0, double u1, int depthb,
 227                            std::vector<std::pair<double, double> > &parameters)
 228 {
 229     intersect_steps ++;
 230     if( deptha > 0 )
 231     {
 232         OldBezier A[2];
 233         a.split(0.5, A[0], A[1]);
 234         double tmid = (t0+t1)*0.5;
 235         deptha--;
 236         if( depthb > 0 )
 237         {
 238             OldBezier B[2];
 239             b.split(0.5, B[0], B[1]);
 240             double umid = (u0+u1)*0.5;
 241             depthb--;
 242             if( intersect_BB( A[0], B[0] ) )
 243                 recursively_intersect( A[0], t0, tmid, deptha,
 244                                       B[0], u0, umid, depthb,
 245                                       parameters );
 246             if( intersect_BB( A[1], B[0] ) )
 247                 recursively_intersect( A[1], tmid, t1, deptha,
 248                                       B[0], u0, umid, depthb,
 249                                       parameters );
 250             if( intersect_BB( A[0], B[1] ) )
 251                 recursively_intersect( A[0], t0, tmid, deptha,
 252                                       B[1], umid, u1, depthb,
 253                                       parameters );
 254             if( intersect_BB( A[1], B[1] ) )
 255                 recursively_intersect( A[1], tmid, t1, deptha,
 256                                       B[1], umid, u1, depthb,
 257                                       parameters );
 258         }
 259         else
 260         {
 261             if( intersect_BB( A[0], b ) )
 262                 recursively_intersect( A[0], t0, tmid, deptha,
 263                                       b, u0, u1, depthb,
 264                                       parameters );
 265             if( intersect_BB( A[1], b ) )
 266                 recursively_intersect( A[1], tmid, t1, deptha,
 267                                       b, u0, u1, depthb,
 268                                       parameters );
 269         }
 270     }
 271     else
 272         if( depthb > 0 )
 273         {
 274             OldBezier B[2];
 275             b.split(0.5, B[0], B[1]);
 276             double umid = (u0 + u1)*0.5;
 277             depthb--;
 278             if( intersect_BB( a, B[0] ) )
 279                 recursively_intersect( a, t0, t1, deptha,
 280                                       B[0], u0, umid, depthb,
 281                                       parameters );
 282             if( intersect_BB( a, B[1] ) )
 283                 recursively_intersect( a, t0, t1, deptha,
 284                                       B[0], umid, u1, depthb,
 285                                       parameters );
 286         }
 287         else // Both segments are fully subdivided; now do line segments
 288         {
 289             double xlk = a.p.back()[X] - a.p[0][X];
 290             double ylk = a.p.back()[Y] - a.p[0][Y];
 291             double xnm = b.p.back()[X] - b.p[0][X];
 292             double ynm = b.p.back()[Y] - b.p[0][Y];
 293             double xmk = b.p[0][X] - a.p[0][X];
 294             double ymk = b.p[0][Y] - a.p[0][Y];
 295             double det = xnm * ylk - ynm * xlk;
 296             if( 1.0 + det == 1.0 )
 297                 return;
 298             else
 299             {
 300                 double detinv = 1.0 / det;
 301                 double s = ( xnm * ymk - ynm *xmk ) * detinv;
 302                 double t = ( xlk * ymk - ylk * xmk ) * detinv;
 303                 if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
 304                     return;
 305                 parameters.push_back(std::pair<double, double>(t0 + s * ( t1 - t0 ),
 306                                                          u0 + t * ( u1 - u0 )));
 307             }
 308         }
 309 }
 310
 311 inline double log4( double x ) { return log(x)/log(4.); }
 312
 313 /*
 314  * Wang's theorem is used to estimate the level of subdivision required,
 315  * but only if the bounding boxes interfere at the top level.
 316  * Assuming there is a possible intersection, recursively_intersect is
 317  * used to find all the parameters corresponding to intersection points.
 318  * these are then sorted and returned in an array.
 319  */
 320
 321 double Lmax(Point p) {
 322     return std::max(fabs(p[X]), fabs(p[Y]));
 323 }
 324
 325 unsigned wangs_theorem(OldBezier a) {
 326     return 12; // seems a good approximation!
 327     double la1 = Lmax( ( a.p[2] - a.p[1] ) - (a.p[1] - a.p[0]) );
 328     double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
 329     double l0 = std::max(la1, la2);
 330     unsigned ra;
 331     if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
 332         ra = 0;
 333     else
 334         ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
 335     std::cout << ra << std::endl;
 336     return ra;
 337 }
 338
 339 std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezier b)
 340 {
 341     std::vector<std::pair<double, double> > parameters;
 342     if( intersect_BB( a, b ) )
 343     {
 344         recursively_intersect( a, 0., 1., wangs_theorem(a),
 345                                b, 0., 1., wangs_theorem(b),
 346                                parameters);
 347     }
 348     std::sort(parameters.begin(), parameters.end());
 349     return parameters;
 350 }
 351 };
 352
 353 /*
 354   Local Variables:
 355   mode:c++
 356   c-file-style:"stroustrup"
 357   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 358   indent-tabs-mode:nil
 359   fill-column:99
 360   End:
 361 */
 362 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :