src/2geom/basic-intersection.cpp

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