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