src/helper/geom.cpp

   1 #define INKSCAPE_HELPER_GEOM_CPP
   2
   3 /**
   4  * Specific geometry functions for Inkscape, not provided my lib2geom.
   5  *
   6  * Author:
   7  *   Johan Engelen <goejendaagh@zonnet.nl>
   8  *
   9  * Copyright (C) 2008 Johan Engelen
  10  *
  11  * Released under GNU GPL
  12  */
  13
  14 #include "helper/geom.h"
  15 #include <typeinfo>
  16 #include <2geom/pathvector.h>
  17 #include <2geom/path.h>
  18 #include <2geom/bezier-curve.h>
  19 #include <2geom/hvlinesegment.h>
  20 #include <2geom/transforms.h>
  21 #include <2geom/rect.h>
  22 #include <2geom/coord.h>
  23 #include <2geom/sbasis-to-bezier.h>
  24 #include <libnr/nr-convert2geom.h>
  25 #include <glibmm.h>
  26
  27 using Geom::X;
  28 using Geom::Y;
  29
  30 #define NR_HUGE   1e18
  31
  32 //#################################################################################
  33 // BOUNDING BOX CALCULATIONS
  34
  35 /* Fast bbox calculation */
  36 /* Thanks to Nathan Hurst for suggesting it */
  37 static void
  38 cubic_bbox (Geom::Coord x000, Geom::Coord y000, Geom::Coord x001, Geom::Coord y001, Geom::Coord x011, Geom::Coord y011, Geom::Coord x111, Geom::Coord y111, Geom::Rect &bbox)
  39 {
  40     Geom::Coord a, b, c, D;
  41
  42     bbox[0].extendTo(x111);
  43     bbox[1].extendTo(y111);
  44
  45     // It already contains (x000,y000) and (x111,y111)
  46     // All points of the Bezier lie in the convex hull of (x000,y000), (x001,y001), (x011,y011) and (x111,y111)
  47     // So, if it also contains (x001,y001) and (x011,y011) we don't have to compute anything else!
  48     // Note that we compute it for the X and Y range separately to make it easier to use them below
  49     bool containsXrange = bbox[0].contains(x001) && bbox[0].contains(x011);
  50     bool containsYrange = bbox[1].contains(y001) && bbox[1].contains(y011);
  51
  52     /*
  53      * xttt = s * (s * (s * x000 + t * x001) + t * (s * x001 + t * x011)) + t * (s * (s * x001 + t * x011) + t * (s * x011 + t * x111))
  54      * xttt = s * (s2 * x000 + s * t * x001 + t * s * x001 + t2 * x011) + t * (s2 * x001 + s * t * x011 + t * s * x011 + t2 * x111)
  55      * xttt = s * (s2 * x000 + 2 * st * x001 + t2 * x011) + t * (s2 * x001 + 2 * st * x011 + t2 * x111)
  56      * xttt = s3 * x000 + 2 * s2t * x001 + st2 * x011 + s2t * x001 + 2st2 * x011 + t3 * x111
  57      * xttt = s3 * x000 + 3s2t * x001 + 3st2 * x011 + t3 * x111
  58      * xttt = s3 * x000 + (1 - s) 3s2 * x001 + (1 - s) * (1 - s) * 3s * x011 + (1 - s) * (1 - s) * (1 - s) * x111
  59      * xttt = s3 * x000 + (3s2 - 3s3) * x001 + (3s - 6s2 + 3s3) * x011 + (1 - 2s + s2 - s + 2s2 - s3) * x111
  60      * xttt = (x000 - 3 * x001 + 3 * x011 -     x111) * s3 +
  61      *        (       3 * x001 - 6 * x011 + 3 * x111) * s2 +
  62      *        (                  3 * x011 - 3 * x111) * s  +
  63      *        (                                 x111)
  64      * xttt' = (3 * x000 - 9 * x001 +  9 * x011 - 3 * x111) * s2 +
  65      *         (           6 * x001 - 12 * x011 + 6 * x111) * s  +
  66      *         (                       3 * x011 - 3 * x111)
  67      */
  68
  69     if (!containsXrange) {
  70         a = 3 * x000 - 9 * x001 + 9 * x011 - 3 * x111;
  71         b = 6 * x001 - 12 * x011 + 6 * x111;
  72         c = 3 * x011 - 3 * x111;
  73
  74         /*
  75         * s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a;
  76         */
  77         if (fabs (a) < Geom::EPSILON) {
  78             /* s = -c / b */
  79             if (fabs (b) > Geom::EPSILON) {
  80                 double s, t, xttt;
  81                 s = -c / b;
  82                 if ((s > 0.0) && (s < 1.0)) {
  83                     t = 1.0 - s;
  84                     xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
  85                     bbox[0].extendTo(xttt);
  86                 }
  87             }
  88         } else {
  89             /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
  90             D = b * b - 4 * a * c;
  91             if (D >= 0.0) {
  92                 Geom::Coord d, s, t, xttt;
  93                 /* Have solution */
  94                 d = sqrt (D);
  95                 s = (-b + d) / (2 * a);
  96                 if ((s > 0.0) && (s < 1.0)) {
  97                     t = 1.0 - s;
  98                     xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
  99                     bbox[0].extendTo(xttt);
 100                 }
 101                 s = (-b - d) / (2 * a);
 102                 if ((s > 0.0) && (s < 1.0)) {
 103                     t = 1.0 - s;
 104                     xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
 105                     bbox[0].extendTo(xttt);
 106                 }
 107             }
 108         }
 109     }
 110
 111     if (!containsYrange) {
 112         a = 3 * y000 - 9 * y001 + 9 * y011 - 3 * y111;
 113         b = 6 * y001 - 12 * y011 + 6 * y111;
 114         c = 3 * y011 - 3 * y111;
 115
 116         if (fabs (a) < Geom::EPSILON) {
 117             /* s = -c / b */
 118             if (fabs (b) > Geom::EPSILON) {
 119                 double s, t, yttt;
 120                 s = -c / b;
 121                 if ((s > 0.0) && (s < 1.0)) {
 122                     t = 1.0 - s;
 123                     yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
 124                     bbox[1].extendTo(yttt);
 125                 }
 126             }
 127         } else {
 128             /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
 129             D = b * b - 4 * a * c;
 130             if (D >= 0.0) {
 131                 Geom::Coord d, s, t, yttt;
 132                 /* Have solution */
 133                 d = sqrt (D);
 134                 s = (-b + d) / (2 * a);
 135                 if ((s > 0.0) && (s < 1.0)) {
 136                     t = 1.0 - s;
 137                     yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
 138                     bbox[1].extendTo(yttt);
 139                 }
 140                 s = (-b - d) / (2 * a);
 141                 if ((s > 0.0) && (s < 1.0)) {
 142                     t = 1.0 - s;
 143                     yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
 144                     bbox[1].extendTo(yttt);
 145                 }
 146             }
 147         }
 148     }
 149 }
 150
 151 Geom::Rect
 152 bounds_fast_transformed(Geom::PathVector const & pv, Geom::Matrix const & t)
 153 {
 154     return bounds_exact_transformed(pv, t); //use this as it is faster for now! :)
 155 //    return Geom::bounds_fast(pv * t);
 156 }
 157
 158 Geom::Rect
 159 bounds_exact_transformed(Geom::PathVector const & pv, Geom::Matrix const & t)
 160 {
 161     Geom::Rect bbox;
 162
 163     if (pv.empty())
 164         return bbox;
 165
 166     Geom::Point initial = pv.front().initialPoint() * t;
 167     bbox = Geom::Rect(initial, initial);        // obtain well defined bbox as starting point to unionWith
 168
 169     for (Geom::PathVector::const_iterator it = pv.begin(); it != pv.end(); ++it) {
 170         bbox.expandTo(it->initialPoint() * t);
 171
 172         // don't loop including closing segment, since that segment can never increase the bbox
 173         for (Geom::Path::const_iterator cit = it->begin(); cit != it->end_open(); ++cit) {
 174             Geom::Curve const &c = *cit;
 175
 176             if( dynamic_cast<Geom::LineSegment const*>(&c) ||
 177                 dynamic_cast<Geom::HLineSegment const*>(&c) ||
 178                 dynamic_cast<Geom::VLineSegment const*>(&c)    )
 179             {
 180                 bbox.expandTo( c.finalPoint() * t );
 181             }
 182             else if(Geom::CubicBezier const *cubic_bezier = dynamic_cast<Geom::CubicBezier const  *>(&c))
 183             {
 184                 Geom::Point c0 = (*cubic_bezier)[0] * t;
 185                 Geom::Point c1 = (*cubic_bezier)[1] * t;
 186                 Geom::Point c2 = (*cubic_bezier)[2] * t;
 187                 Geom::Point c3 = (*cubic_bezier)[3] * t;
 188                 cubic_bbox( c0[0], c0[1],
 189                             c1[0], c1[1],
 190                             c2[0], c2[1],
 191                             c3[0], c3[1],
 192                             bbox );
 193             }
 194             else
 195             {
 196                 // should handle all not-so-easy curves:
 197                 Geom::Curve *ctemp = cit->transformed(t);
 198                 bbox.unionWith( ctemp->boundsExact());
 199                 delete ctemp;
 200             }
 201         }
 202     }
 203     //return Geom::bounds_exact(pv * t);
 204     return bbox;
 205 }
 206
 207
 208
 209 static void
 210 geom_line_wind_distance (Geom::Coord x0, Geom::Coord y0, Geom::Coord x1, Geom::Coord y1, Geom::Point const &pt, int *wind, Geom::Coord *best)
 211 {
 212     Geom::Coord Ax, Ay, Bx, By, Dx, Dy, s;
 213     Geom::Coord dist2;
 214
 215     /* Find distance */
 216     Ax = x0;
 217     Ay = y0;
 218     Bx = x1;
 219     By = y1;
 220     Dx = x1 - x0;
 221     Dy = y1 - y0;
 222     const Geom::Coord Px = pt[X];
 223     const Geom::Coord Py = pt[Y];
 224
 225     if (best) {
 226         s = ((Px - Ax) * Dx + (Py - Ay) * Dy) / (Dx * Dx + Dy * Dy);
 227         if (s <= 0.0) {
 228             dist2 = (Px - Ax) * (Px - Ax) + (Py - Ay) * (Py - Ay);
 229         } else if (s >= 1.0) {
 230             dist2 = (Px - Bx) * (Px - Bx) + (Py - By) * (Py - By);
 231         } else {
 232             Geom::Coord Qx, Qy;
 233             Qx = Ax + s * Dx;
 234             Qy = Ay + s * Dy;
 235             dist2 = (Px - Qx) * (Px - Qx) + (Py - Qy) * (Py - Qy);
 236         }
 237
 238         if (dist2 < (*best * *best)) *best = sqrt (dist2);
 239     }
 240
 241     if (wind) {
 242         /* Find wind */
 243         if ((Ax >= Px) && (Bx >= Px)) return;
 244         if ((Ay >= Py) && (By >= Py)) return;
 245         if ((Ay < Py) && (By < Py)) return;
 246         if (Ay == By) return;
 247         /* Ctach upper y bound */
 248         if (Ay == Py) {
 249             if (Ax < Px) *wind -= 1;
 250             return;
 251         } else if (By == Py) {
 252             if (Bx < Px) *wind += 1;
 253             return;
 254         } else {
 255             Geom::Coord Qx;
 256             /* Have to calculate intersection */
 257             Qx = Ax + Dx * (Py - Ay) / Dy;
 258             if (Qx < Px) {
 259                 *wind += (Dy > 0.0) ? 1 : -1;
 260             }
 261         }
 262     }
 263 }
 264
 265 static void
 266 geom_cubic_bbox_wind_distance (Geom::Coord x000, Geom::Coord y000,
 267                  Geom::Coord x001, Geom::Coord y001,
 268                  Geom::Coord x011, Geom::Coord y011,
 269                  Geom::Coord x111, Geom::Coord y111,
 270                  Geom::Point const &pt,
 271                  Geom::Rect *bbox, int *wind, Geom::Coord *best,
 272                  Geom::Coord tolerance)
 273 {
 274     Geom::Coord x0, y0, x1, y1, len2;
 275     int needdist, needwind, needline;
 276
 277     const Geom::Coord Px = pt[X];
 278     const Geom::Coord Py = pt[Y];
 279
 280     needdist = 0;
 281     needwind = 0;
 282     needline = 0;
 283
 284     if (bbox) cubic_bbox (x000, y000, x001, y001, x011, y011, x111, y111, *bbox);
 285
 286     x0 = MIN (x000, x001);
 287     x0 = MIN (x0, x011);
 288     x0 = MIN (x0, x111);
 289     y0 = MIN (y000, y001);
 290     y0 = MIN (y0, y011);
 291     y0 = MIN (y0, y111);
 292     x1 = MAX (x000, x001);
 293     x1 = MAX (x1, x011);
 294     x1 = MAX (x1, x111);
 295     y1 = MAX (y000, y001);
 296     y1 = MAX (y1, y011);
 297     y1 = MAX (y1, y111);
 298
 299     if (best) {
 300         /* Quickly adjust to endpoints */
 301         len2 = (x000 - Px) * (x000 - Px) + (y000 - Py) * (y000 - Py);
 302         if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);
 303         len2 = (x111 - Px) * (x111 - Px) + (y111 - Py) * (y111 - Py);
 304         if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);
 305
 306         if (((x0 - Px) < *best) && ((y0 - Py) < *best) && ((Px - x1) < *best) && ((Py - y1) < *best)) {
 307             /* Point is inside sloppy bbox */
 308             /* Now we have to decide, whether subdivide */
 309             /* fixme: (Lauris) */
 310             if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
 311                 needdist = 1;
 312             } else {
 313                 needline = 1;
 314             }
 315         }
 316     }
 317     if (!needdist && wind) {
 318         if ((y1 >= Py) && (y0 < Py) && (x0 < Px)) {
 319             /* Possible intersection at the left */
 320             /* Now we have to decide, whether subdivide */
 321             /* fixme: (Lauris) */
 322             if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
 323                 needwind = 1;
 324             } else {
 325                 needline = 1;
 326             }
 327         }
 328     }
 329
 330     if (needdist || needwind) {
 331         Geom::Coord x00t, x0tt, xttt, x1tt, x11t, x01t;
 332         Geom::Coord y00t, y0tt, yttt, y1tt, y11t, y01t;
 333         Geom::Coord s, t;
 334
 335         t = 0.5;
 336         s = 1 - t;
 337
 338         x00t = s * x000 + t * x001;
 339         x01t = s * x001 + t * x011;
 340         x11t = s * x011 + t * x111;
 341         x0tt = s * x00t + t * x01t;
 342         x1tt = s * x01t + t * x11t;
 343         xttt = s * x0tt + t * x1tt;
 344
 345         y00t = s * y000 + t * y001;
 346         y01t = s * y001 + t * y011;
 347         y11t = s * y011 + t * y111;
 348         y0tt = s * y00t + t * y01t;
 349         y1tt = s * y01t + t * y11t;
 350         yttt = s * y0tt + t * y1tt;
 351
 352         geom_cubic_bbox_wind_distance (x000, y000, x00t, y00t, x0tt, y0tt, xttt, yttt, pt, NULL, wind, best, tolerance);
 353         geom_cubic_bbox_wind_distance (xttt, yttt, x1tt, y1tt, x11t, y11t, x111, y111, pt, NULL, wind, best, tolerance);
 354     } else if (1 || needline) {
 355         geom_line_wind_distance (x000, y000, x111, y111, pt, wind, best);
 356     }
 357 }
 358
 359 static void
 360 geom_curve_bbox_wind_distance(Geom::Curve const & c, Geom::Matrix const &m,
 361                  Geom::Point const &pt,
 362                  Geom::Rect *bbox, int *wind, Geom::Coord *dist,
 363                  Geom::Coord tolerance, Geom::Rect const *viewbox,
 364                  Geom::Point &p0) // pass p0 through as it represents the last endpoint added (the finalPoint of last curve)
 365 {
 366     if( dynamic_cast<Geom::LineSegment const*>(&c) ||
 367         dynamic_cast<Geom::HLineSegment const*>(&c) ||
 368         dynamic_cast<Geom::VLineSegment const*>(&c) )
 369     {
 370         Geom::Point pe = c.finalPoint() * m;
 371         if (bbox) {
 372             bbox->expandTo(pe);
 373         }
 374         if (dist || wind) {
 375             if (wind) { // we need to pick fill, so do what we're told
 376                 geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);
 377             } else { // only stroke is being picked; skip this segment if it's totally outside the viewbox
 378                 Geom::Rect swept(p0, pe);
 379                 if (!viewbox || swept.intersects(*viewbox))
 380                     geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);
 381             }
 382         }
 383         p0 = pe;
 384     }
 385     else if(Geom::CubicBezier const *cubic_bezier = dynamic_cast<Geom::CubicBezier const  *>(&c)) {
 386         Geom::Point p1 = (*cubic_bezier)[1] * m;
 387         Geom::Point p2 = (*cubic_bezier)[2] * m;
 388         Geom::Point p3 = (*cubic_bezier)[3] * m;
 389
 390         // get approximate bbox from handles (convex hull property of beziers):
 391         Geom::Rect swept(p0, p3);
 392         swept.expandTo(p1);
 393         swept.expandTo(p2);
 394
 395         if (!viewbox || swept.intersects(*viewbox)) { // we see this segment, so do full processing
 396             geom_cubic_bbox_wind_distance ( p0[X], p0[Y],
 397                                             p1[X], p1[Y],
 398                                             p2[X], p2[Y],
 399                                             p3[X], p3[Y],
 400                                             pt,
 401                                             bbox, wind, dist, tolerance);
 402         } else {
 403             if (wind) { // if we need fill, we can just pretend it's a straight line
 404                 geom_line_wind_distance (p0[X], p0[Y], p3[X], p3[Y], pt, wind, dist);
 405             } else { // otherwise, skip it completely
 406             }
 407         }
 408         p0 = p3;
 409     } else {
 410         //this case handles sbasis as well as all other curve types
 411         Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(c.toSBasis(), 0.1);
 412
 413         //recurse to convert the new path resulting from the sbasis to svgd
 414         for(Geom::Path::iterator iter = sbasis_path.begin(); iter != sbasis_path.end(); ++iter) {
 415             geom_curve_bbox_wind_distance(*iter, m, pt, bbox, wind, dist, tolerance, viewbox, p0);
 416         }
 417     }
 418 }
 419
 420 /* Calculates...
 421    and returns ... in *wind and the distance to ... in *dist.
 422    Returns bounding box in *bbox if bbox!=NULL.
 423  */
 424 void
 425 pathv_matrix_point_bbox_wind_distance (Geom::PathVector const & pathv, Geom::Matrix const &m, Geom::Point const &pt,
 426                          Geom::Rect *bbox, int *wind, Geom::Coord *dist,
 427                          Geom::Coord tolerance, Geom::Rect const *viewbox)
 428 {
 429     if (pathv.empty()) {
 430         if (wind) *wind = 0;
 431         if (dist) *dist = NR_HUGE;
 432         return;
 433     }
 434
 435     // remember last point of last curve
 436     Geom::Point p0(0,0);
 437
 438     // remembering the start of subpath
 439     Geom::Point p_start(0,0);
 440     bool start_set = false;
 441
 442     for (Geom::PathVector::const_iterator it = pathv.begin(); it != pathv.end(); ++it) {
 443
 444         if (start_set) { // this is a new subpath
 445             if (wind && (p0 != p_start)) // for correct fill picking, each subpath must be closed
 446                 geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);
 447         }
 448         p0 = it->initialPoint() * m;
 449         p_start = p0;
 450         start_set = true;
 451         if (bbox) {
 452             bbox->expandTo(p0);
 453         }
 454
 455         // loop including closing segment if path is closed
 456         for (Geom::Path::const_iterator cit = it->begin(); cit != it->end_default(); ++cit) {
 457             geom_curve_bbox_wind_distance(*cit, m, pt, bbox, wind, dist, tolerance, viewbox, p0);
 458         }
 459     }
 460
 461     if (start_set) {
 462         if (wind && (p0 != p_start)) // for correct picking, each subpath must be closed
 463             geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);
 464     }
 465 }
 466
 467 // temporary wrapper
 468 void
 469 pathv_matrix_point_bbox_wind_distance (Geom::PathVector const & pathv, NR::Matrix const &m, NR::Point const &pt,
 470                          NR::Rect *bbox, int *wind, NR::Coord *dist,
 471                          NR::Coord tolerance, NR::Rect const *viewbox)
 472 {
 473     Geom::Rect _bbox;
 474     if (bbox)
 475         _bbox = to_2geom(*bbox);
 476     Geom::Coord _dist;
 477     if (dist)
 478         _dist = *dist;
 479     Geom::Rect _viewbox;
 480     if (viewbox)
 481         _viewbox = to_2geom(*viewbox);
 482
 483     pathv_matrix_point_bbox_wind_distance( pathv, to_2geom(m), to_2geom(pt),
 484                                            bbox ? &_bbox : NULL,
 485                                            wind,
 486                                            dist ? &_dist : NULL,
 487                                            tolerance,
 488                                            viewbox ? &_viewbox : NULL );
 489
 490     if (bbox)
 491         *bbox = from_2geom(_bbox);
 492     if (dist)
 493         *dist = _dist;
 494 }
 495 //#################################################################################
 496
 497 /*
 498  * Converts all segments in all paths to Geom::LineSegment or Geom::HLineSegment or
 499  * Geom::VLineSegment or Geom::CubicBezier.
 500  */
 501 Geom::PathVector
 502 pathv_to_linear_and_cubic_beziers( Geom::PathVector const &pathv )
 503 {
 504     Geom::PathVector output;
 505
 506     for (Geom::PathVector::const_iterator pit = pathv.begin(); pit != pathv.end(); ++pit) {
 507         output.push_back( Geom::Path() );
 508         output.back().start( pit->initialPoint() );
 509         output.back().close( pit->closed() );
 510
 511         for (Geom::Path::const_iterator cit = pit->begin(); cit != pit->end_open(); ++cit) {
 512             if( dynamic_cast<Geom::CubicBezier const*>(&*cit) ||
 513                 dynamic_cast<Geom::LineSegment const*>(&*cit) ||
 514                 dynamic_cast<Geom::HLineSegment const*>(&*cit) ||
 515                 dynamic_cast<Geom::VLineSegment const*>(&*cit) )
 516             {
 517                 output.back().append(*cit);
 518             }
 519             else {
 520                 // convert all other curve types to cubicbeziers
 521                 Geom::Path cubicbezier_path = Geom::cubicbezierpath_from_sbasis(cit->toSBasis(), 0.1);
 522                 output.back().append(cubicbezier_path);
 523             }
 524         }
 525     }
 526
 527     return output;
 528 }
 529
 530
 531
 532
 533 /*
 534   Local Variables:
 535   mode:c++
 536   c-file-style:"stroustrup"
 537   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 538   indent-tabs-mode:nil
 539   fill-column:99
 540   End:
 541 */
 542 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :