src/libnr/nr-path.cpp

   1 #define __NR_PATH_C__
   2
   3 /*
   4  * Pixel buffer rendering library
   5  *
   6  * Authors:
   7  *   Lauris Kaplinski <lauris@kaplinski.com>
   8  *
   9  * This code is in public domain
  10  */
  11
  12 #include <glib/gmem.h>
  13 #include "n-art-bpath.h"
  14 #include "nr-rect.h"
  15 #include "nr-path.h"
  16
  17 static void nr_curve_bbox (NR::Coord x000, NR::Coord y000, NR::Coord x001, NR::Coord y001, NR::Coord x011, NR::Coord y011, NR::Coord x111, NR::Coord y111, NRRect *bbox);
  18
  19 static void nr_curve_bbox(NR::Point const p000, NR::Point const p001,
  20                           NR::Point const p011, NR::Point const p111,
  21                           NRRect *bbox);
  22
  23 NRBPath *nr_path_duplicate_transform(NRBPath *d, NRBPath *s, NRMatrix const *transform)
  24 {
  25         int i;
  26
  27         if (!s->path) {
  28                 d->path = NULL;
  29                 return d;
  30         }
  31
  32         i = 0;
  33         while (s->path[i].code != NR_END) i += 1;
  34
  35         d->path = g_new (NArtBpath, i + 1);
  36
  37         i = 0;
  38         while (s->path[i].code != NR_END) {
  39                 d->path[i].code = s->path[i].code;
  40                 if (s->path[i].code == NR_CURVETO) {
  41                         d->path[i].x1 = NR_MATRIX_DF_TRANSFORM_X (transform, s->path[i].x1, s->path[i].y1);
  42                         d->path[i].y1 = NR_MATRIX_DF_TRANSFORM_Y (transform, s->path[i].x1, s->path[i].y1);
  43                         d->path[i].x2 = NR_MATRIX_DF_TRANSFORM_X (transform, s->path[i].x2, s->path[i].y2);
  44                         d->path[i].y2 = NR_MATRIX_DF_TRANSFORM_Y (transform, s->path[i].x2, s->path[i].y2);
  45                 }
  46                 d->path[i].x3 = NR_MATRIX_DF_TRANSFORM_X (transform, s->path[i].x3, s->path[i].y3);
  47                 d->path[i].y3 = NR_MATRIX_DF_TRANSFORM_Y (transform, s->path[i].x3, s->path[i].y3);
  48                 i += 1;
  49         }
  50         d->path[i].code = NR_END;
  51
  52         return d;
  53 }
  54
  55 NRBPath *nr_path_duplicate_transform(NRBPath *d, NRBPath *s, NR::Matrix const transform) {
  56         NRMatrix tr = transform;
  57         return nr_path_duplicate_transform(d, s, &tr);
  58 }
  59
  60 NArtBpath* nr_artpath_affine(NArtBpath *s, NR::Matrix const &aff) {
  61         NRBPath bp, abp;
  62         bp.path = s;
  63         nr_path_duplicate_transform(&abp, &bp, aff);
  64         return abp.path;
  65 }
  66
  67 static void
  68 nr_line_wind_distance (NR::Coord x0, NR::Coord y0, NR::Coord x1, NR::Coord y1, NR::Point &pt, int *wind, NR::Coord *best)
  69 {
  70         NR::Coord Ax, Ay, Bx, By, Dx, Dy, s;
  71         NR::Coord dist2;
  72
  73         /* Find distance */
  74         Ax = x0;
  75         Ay = y0;
  76         Bx = x1;
  77         By = y1;
  78         Dx = x1 - x0;
  79         Dy = y1 - y0;
  80         const NR::Coord Px = pt[NR::X];
  81         const NR::Coord Py = pt[NR::Y];
  82
  83         if (best) {
  84                 s = ((Px - Ax) * Dx + (Py - Ay) * Dy) / (Dx * Dx + Dy * Dy);
  85                 if (s <= 0.0) {
  86                         dist2 = (Px - Ax) * (Px - Ax) + (Py - Ay) * (Py - Ay);
  87                 } else if (s >= 1.0) {
  88                         dist2 = (Px - Bx) * (Px - Bx) + (Py - By) * (Py - By);
  89                 } else {
  90                         NR::Coord Qx, Qy;
  91                         Qx = Ax + s * Dx;
  92                         Qy = Ay + s * Dy;
  93                         dist2 = (Px - Qx) * (Px - Qx) + (Py - Qy) * (Py - Qy);
  94                 }
  95
  96                 if (dist2 < (*best * *best)) *best = sqrt (dist2);
  97         }
  98
  99         if (wind) {
 100                 /* Find wind */
 101                 if ((Ax >= Px) && (Bx >= Px)) return;
 102                 if ((Ay >= Py) && (By >= Py)) return;
 103                 if ((Ay < Py) && (By < Py)) return;
 104                 if (Ay == By) return;
 105                 /* Ctach upper y bound */
 106                 if (Ay == Py) {
 107                         if (Ax < Px) *wind -= 1;
 108                         return;
 109                 } else if (By == Py) {
 110                         if (Bx < Px) *wind += 1;
 111                         return;
 112                 } else {
 113                         NR::Coord Qx;
 114                         /* Have to calculate intersection */
 115                         Qx = Ax + Dx * (Py - Ay) / Dy;
 116                         if (Qx < Px) {
 117                                 *wind += (Dy > 0.0) ? 1 : -1;
 118                         }
 119                 }
 120         }
 121 }
 122
 123 static void
 124 nr_curve_bbox_wind_distance (NR::Coord x000, NR::Coord y000,
 125                              NR::Coord x001, NR::Coord y001,
 126                              NR::Coord x011, NR::Coord y011,
 127                              NR::Coord x111, NR::Coord y111,
 128                              NR::Point &pt,
 129                              NRRect *bbox, int *wind, NR::Coord *best,
 130                              NR::Coord tolerance)
 131 {
 132         NR::Coord x0, y0, x1, y1, len2;
 133         int needdist, needwind, needline;
 134
 135         const NR::Coord Px = pt[NR::X];
 136         const NR::Coord Py = pt[NR::Y];
 137
 138         needdist = 0;
 139         needwind = 0;
 140         needline = 0;
 141
 142         if (bbox) nr_curve_bbox (x000, y000, x001, y001, x011, y011, x111, y111, bbox);
 143
 144         x0 = MIN (x000, x001);
 145         x0 = MIN (x0, x011);
 146         x0 = MIN (x0, x111);
 147         y0 = MIN (y000, y001);
 148         y0 = MIN (y0, y011);
 149         y0 = MIN (y0, y111);
 150         x1 = MAX (x000, x001);
 151         x1 = MAX (x1, x011);
 152         x1 = MAX (x1, x111);
 153         y1 = MAX (y000, y001);
 154         y1 = MAX (y1, y011);
 155         y1 = MAX (y1, y111);
 156
 157         if (best) {
 158                 /* Quicly adjust to endpoints */
 159                 len2 = (x000 - Px) * (x000 - Px) + (y000 - Py) * (y000 - Py);
 160                 if (len2 < (*best * *best)) *best = (NR::Coord) sqrt (len2);
 161                 len2 = (x111 - Px) * (x111 - Px) + (y111 - Py) * (y111 - Py);
 162                 if (len2 < (*best * *best)) *best = (NR::Coord) sqrt (len2);
 163
 164                 if (((x0 - Px) < *best) && ((y0 - Py) < *best) && ((Px - x1) < *best) && ((Py - y1) < *best)) {
 165                         /* Point is inside sloppy bbox */
 166                         /* Now we have to decide, whether subdivide */
 167                         /* fixme: (Lauris) */
 168                         if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
 169                                 needdist = 1;
 170                         } else {
 171                                 needline = 1;
 172                         }
 173                 }
 174         }
 175         if (!needdist && wind) {
 176                 if ((y1 >= Py) && (y0 < Py) && (x0 < Px)) {
 177                         /* Possible intersection at the left */
 178                         /* Now we have to decide, whether subdivide */
 179                         /* fixme: (Lauris) */
 180                         if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
 181                                 needwind = 1;
 182                         } else {
 183                                 needline = 1;
 184                         }
 185                 }
 186         }
 187
 188         if (needdist || needwind) {
 189                 NR::Coord x00t, x0tt, xttt, x1tt, x11t, x01t;
 190                 NR::Coord y00t, y0tt, yttt, y1tt, y11t, y01t;
 191                 NR::Coord s, t;
 192
 193                 t = 0.5;
 194                 s = 1 - t;
 195
 196                 x00t = s * x000 + t * x001;
 197                 x01t = s * x001 + t * x011;
 198                 x11t = s * x011 + t * x111;
 199                 x0tt = s * x00t + t * x01t;
 200                 x1tt = s * x01t + t * x11t;
 201                 xttt = s * x0tt + t * x1tt;
 202
 203                 y00t = s * y000 + t * y001;
 204                 y01t = s * y001 + t * y011;
 205                 y11t = s * y011 + t * y111;
 206                 y0tt = s * y00t + t * y01t;
 207                 y1tt = s * y01t + t * y11t;
 208                 yttt = s * y0tt + t * y1tt;
 209
 210                 nr_curve_bbox_wind_distance (x000, y000, x00t, y00t, x0tt, y0tt, xttt, yttt, pt, NULL, wind, best, tolerance);
 211                 nr_curve_bbox_wind_distance (xttt, yttt, x1tt, y1tt, x11t, y11t, x111, y111, pt, NULL, wind, best, tolerance);
 212         } else if (1 || needline) {
 213                 nr_line_wind_distance (x000, y000, x111, y111, pt, wind, best);
 214         }
 215 }
 216
 217 void
 218 nr_path_matrix_point_bbox_wind_distance (NRBPath *bpath, NR::Matrix const &m, NR::Point &pt,
 219                                              NRRect *bbox, int *wind, NR::Coord *dist,
 220                                                  NR::Coord tolerance, NR::Rect *viewbox)
 221 {
 222         if (!bpath->path) {
 223                 if (wind) *wind = 0;
 224                 if (dist) *dist = NR_HUGE;
 225                 return;
 226         }
 227
 228         NR::Coord x0 = 0;
 229         NR::Coord y0 = 0;
 230         NR::Coord x1 = 0;
 231         NR::Coord y1 = 0;
 232         NR::Coord x2 = 0;
 233         NR::Coord y2 = 0;
 234         NR::Coord x3 = 0;
 235         NR::Coord y3 = 0;
 236
 237         // remembering the start of subpath
 238         NR::Coord x_start = 0, y_start = 0; bool start_set = false;
 239         NR::Rect swept;
 240
 241         for (const NArtBpath *p = bpath->path; p->code != NR_END; p+= 1) {
 242                 switch (p->code) {
 243                 case NR_MOVETO_OPEN:
 244                 case NR_MOVETO:
 245                         if (start_set) { // this is a new subpath
 246                                 if (wind && (x0 != x_start || y0 != y_start)) // for correct fill picking, each subpath must be closed
 247                                         nr_line_wind_distance (x0, y0, x_start, y_start, pt, wind, dist);
 248                         }
 249                         x0 = m[0] * p->x3 + m[2] * p->y3 + m[4];
 250                         y0 = m[1] * p->x3 + m[3] * p->y3 + m[5];
 251                         x_start = x0; y_start = y0; start_set = true;
 252                         if (bbox) {
 253                                 bbox->x0 = (NR::Coord) MIN (bbox->x0, x0);
 254                                 bbox->y0 = (NR::Coord) MIN (bbox->y0, y0);
 255                                 bbox->x1 = (NR::Coord) MAX (bbox->x1, x0);
 256                                 bbox->y1 = (NR::Coord) MAX (bbox->y1, y0);
 257                         }
 258                         break;
 259                 case NR_LINETO:
 260                         x3 = m[0] * p->x3 + m[2] * p->y3 + m[4];
 261                         y3 = m[1] * p->x3 + m[3] * p->y3 + m[5];
 262                         if (bbox) {
 263                                 bbox->x0 = (NR::Coord) MIN (bbox->x0, x3);
 264                                 bbox->y0 = (NR::Coord) MIN (bbox->y0, y3);
 265                                 bbox->x1 = (NR::Coord) MAX (bbox->x1, x3);
 266                                 bbox->y1 = (NR::Coord) MAX (bbox->y1, y3);
 267                         }
 268                         if (dist || wind) {
 269                                 if (wind) { // we need to pick fill, so do what we're told
 270                                         nr_line_wind_distance (x0, y0, x3, y3, pt, wind, dist);
 271                                 } else { // only stroke is being picked; skip this segment if it's totally outside the viewbox
 272                                         swept = NR::Rect(NR::Point(x0, y0), NR::Point(x3, y3));
 273                                         //std::cout << "swept: " << swept;
 274                                         //std::cout << "view: " << *viewbox;
 275                                         //std::cout << "intersects? " << (swept.intersects(*viewbox)? "YES" : "NO") << "\n";
 276                                         if (!viewbox || swept.intersects(*viewbox))
 277                                         nr_line_wind_distance (x0, y0, x3, y3, pt, wind, dist);
 278                                 }
 279                         }
 280                         x0 = x3;
 281                         y0 = y3;
 282                         break;
 283                 case NR_CURVETO:
 284                         {
 285                         x3 = m[0] * p->x3 + m[2] * p->y3 + m[4];
 286                         y3 = m[1] * p->x3 + m[3] * p->y3 + m[5];
 287                         x1 = m[0] * p->x1 + m[2] * p->y1 + m[4];
 288                         y1 = m[1] * p->x1 + m[3] * p->y1 + m[5];
 289                         x2 = m[0] * p->x2 + m[2] * p->y2 + m[4];
 290                         y2 = m[1] * p->x2 + m[3] * p->y2 + m[5];
 291
 292                         swept = NR::Rect(NR::Point(x0, y0), NR::Point(x3, y3));
 293                         swept.expandTo(NR::Point(x1, y1));
 294                         swept.expandTo(NR::Point(x2, y2));
 295
 296                         if (!viewbox || swept.intersects(*viewbox)) { // we see this segment, so do full processing
 297                                 nr_curve_bbox_wind_distance (
 298                        x0, y0,
 299                        x1, y1,
 300                        x2, y2,
 301                                                      x3, y3,
 302                                                      pt,
 303                                                      bbox, wind, dist, tolerance);
 304                         } else {
 305                                 if (wind) { // if we need fill, we can just pretend it's a straight line
 306                                         nr_line_wind_distance (x0, y0, x3, y3, pt, wind, dist);
 307                                 } else { // otherwise, skip it completely
 308                                 }
 309                         }
 310                         x0 = x3;
 311                         y0 = y3;
 312                         }
 313                         break;
 314                 default:
 315                         break;
 316                 }
 317         }
 318
 319         if (start_set) {
 320                 if (wind && (x0 != x_start || y0 != y_start)) // for correct picking, each subpath must be closed
 321                         nr_line_wind_distance (x0, y0, x_start, y_start, pt, wind, dist);
 322         }
 323 }
 324
 325 static void
 326 nr_curve_bbox(NR::Point const p000, NR::Point const p001,
 327               NR::Point const p011, NR::Point const p111,
 328               NRRect *bbox)
 329 {
 330         using NR::X;
 331         using NR::Y;
 332
 333         nr_curve_bbox(p000[X], p000[Y],
 334                       p001[X], p001[Y],
 335                       p011[X], p011[Y],
 336                       p111[X], p111[Y],
 337                       bbox);
 338 }
 339
 340 /* Fast bbox calculation */
 341 /* Thanks to Nathan Hurst for suggesting it */
 342
 343 static void
 344 nr_curve_bbox (NR::Coord x000, NR::Coord y000, NR::Coord x001, NR::Coord y001, NR::Coord x011, NR::Coord y011, NR::Coord x111, NR::Coord y111, NRRect *bbox)
 345 {
 346         NR::Coord a, b, c, D;
 347
 348         bbox->x0 = (NR::Coord) MIN (bbox->x0, x111);
 349         bbox->y0 = (NR::Coord) MIN (bbox->y0, y111);
 350         bbox->x1 = (NR::Coord) MAX (bbox->x1, x111);
 351         bbox->y1 = (NR::Coord) MAX (bbox->y1, y111);
 352
 353         /*
 354          * xttt = s * (s * (s * x000 + t * x001) + t * (s * x001 + t * x011)) + t * (s * (s * x001 + t * x011) + t * (s * x011 + t * x111))
 355          * xttt = s * (s2 * x000 + s * t * x001 + t * s * x001 + t2 * x011) + t * (s2 * x001 + s * t * x011 + t * s * x011 + t2 * x111)
 356          * xttt = s * (s2 * x000 + 2 * st * x001 + t2 * x011) + t * (s2 * x001 + 2 * st * x011 + t2 * x111)
 357          * xttt = s3 * x000 + 2 * s2t * x001 + st2 * x011 + s2t * x001 + 2st2 * x011 + t3 * x111
 358          * xttt = s3 * x000 + 3s2t * x001 + 3st2 * x011 + t3 * x111
 359          * xttt = s3 * x000 + (1 - s) 3s2 * x001 + (1 - s) * (1 - s) * 3s * x011 + (1 - s) * (1 - s) * (1 - s) * x111
 360          * xttt = s3 * x000 + (3s2 - 3s3) * x001 + (3s - 6s2 + 3s3) * x011 + (1 - 2s + s2 - s + 2s2 - s3) * x111
 361          * xttt = (x000 - 3 * x001 + 3 * x011 -     x111) * s3 +
 362          *        (       3 * x001 - 6 * x011 + 3 * x111) * s2 +
 363          *        (                  3 * x011 - 3 * x111) * s  +
 364          *        (                                 x111)
 365          * xttt' = (3 * x000 - 9 * x001 +  9 * x011 - 3 * x111) * s2 +
 366          *         (           6 * x001 - 12 * x011 + 6 * x111) * s  +
 367          *         (                       3 * x011 - 3 * x111)
 368          */
 369
 370         a = 3 * x000 - 9 * x001 + 9 * x011 - 3 * x111;
 371         b = 6 * x001 - 12 * x011 + 6 * x111;
 372         c = 3 * x011 - 3 * x111;
 373
 374         /*
 375          * s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a;
 376          */
 377         if (fabs (a) < NR_EPSILON) {
 378                 /* s = -c / b */
 379                 if (fabs (b) > NR_EPSILON) {
 380                         double s, t, xttt;
 381                         s = -c / b;
 382                         if ((s > 0.0) && (s < 1.0)) {
 383                                 t = 1.0 - s;
 384                                 xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
 385                                 bbox->x0 = (float) MIN (bbox->x0, xttt);
 386                                 bbox->x1 = (float) MAX (bbox->x1, xttt);
 387                         }
 388                 }
 389         } else {
 390                 /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
 391                 D = b * b - 4 * a * c;
 392                 if (D >= 0.0) {
 393                         NR::Coord d, s, t, xttt;
 394                         /* Have solution */
 395                         d = sqrt (D);
 396                         s = (-b + d) / (2 * a);
 397                         if ((s > 0.0) && (s < 1.0)) {
 398                                 t = 1.0 - s;
 399                                 xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
 400                                 bbox->x0 = (NR::Coord) MIN (bbox->x0, xttt);
 401                                 bbox->x1 = (NR::Coord) MAX (bbox->x1, xttt);
 402                         }
 403                         s = (-b - d) / (2 * a);
 404                         if ((s > 0.0) && (s < 1.0)) {
 405                                 t = 1.0 - s;
 406                                 xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
 407                                 bbox->x0 = (NR::Coord) MIN (bbox->x0, xttt);
 408                                 bbox->x1 = (NR::Coord) MAX (bbox->x1, xttt);
 409                         }
 410                 }
 411         }
 412
 413         a = 3 * y000 - 9 * y001 + 9 * y011 - 3 * y111;
 414         b = 6 * y001 - 12 * y011 + 6 * y111;
 415         c = 3 * y011 - 3 * y111;
 416
 417         if (fabs (a) < NR_EPSILON) {
 418                 /* s = -c / b */
 419                 if (fabs (b) > NR_EPSILON) {
 420                         double s, t, yttt;
 421                         s = -c / b;
 422                         if ((s > 0.0) && (s < 1.0)) {
 423                                 t = 1.0 - s;
 424                                 yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
 425                                 bbox->y0 = (float) MIN (bbox->y0, yttt);
 426                                 bbox->y1 = (float) MAX (bbox->y1, yttt);
 427                         }
 428                 }
 429         } else {
 430                 /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
 431                 D = b * b - 4 * a * c;
 432                 if (D >= 0.0) {
 433                         NR::Coord d, s, t, yttt;
 434                         /* Have solution */
 435                         d = sqrt (D);
 436                         s = (-b + d) / (2 * a);
 437                         if ((s > 0.0) && (s < 1.0)) {
 438                                 t = 1.0 - s;
 439                                 yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
 440                                 bbox->y0 = (NR::Coord) MIN (bbox->y0, yttt);
 441                                 bbox->y1 = (NR::Coord) MAX (bbox->y1, yttt);
 442                         }
 443                         s = (-b - d) / (2 * a);
 444                         if ((s > 0.0) && (s < 1.0)) {
 445                                 t = 1.0 - s;
 446                                 yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
 447                                 bbox->y0 = (NR::Coord) MIN (bbox->y0, yttt);
 448                                 bbox->y1 = (NR::Coord) MAX (bbox->y1, yttt);
 449                         }
 450                 }
 451         }
 452 }
 453
 454 void
 455 nr_path_matrix_bbox_union(NRBPath const *bpath, NR::Matrix const &m,
 456                           NRRect *bbox)
 457 {
 458     using NR::X;
 459     using NR::Y;
 460
 461     if (!bpath->path) {
 462         return;
 463     }
 464
 465     NR::Point prev0(0, 0);
 466     for (NArtBpath const *p = bpath->path; p->code != NR_END; ++p) {
 467         switch (p->code) {
 468             case NR_MOVETO_OPEN:
 469             case NR_MOVETO: {
 470                 NR::Point const c3(p->x3,
 471                                    p->y3);
 472                 prev0 = c3 * m;
 473                 nr_rect_union_pt(bbox, prev0);
 474                 break;
 475             }
 476
 477             case NR_LINETO: {
 478                 NR::Point const c3(p->x3,
 479                                    p->y3);
 480                 NR::Point const endPt( c3 * m );
 481                 nr_rect_union_pt(bbox, endPt);
 482                 prev0 = endPt;
 483                 break;
 484             }
 485
 486             case NR_CURVETO: {
 487                 NR::Point const c1(p->x1, p->y1);
 488                 NR::Point const c2(p->x2, p->y2);
 489                 NR::Point const c3(p->x3, p->y3);
 490                 NR::Point const endPt( c3 * m );
 491                 nr_curve_bbox(prev0,
 492                               c1 * m,
 493                               c2 * m,
 494                               endPt,
 495                               bbox);
 496                 prev0 = endPt;
 497                 break;
 498             }
 499
 500             default:
 501                 break;
 502         }
 503     }
 504 }
 505
 506 NArtBpath *nr_path_from_rect(NRRect const &r)
 507 {
 508     NArtBpath *path = g_new (NArtBpath, 6);
 509
 510     path[0].code = NR_MOVETO;
 511     path[0].setC(3, NR::Point(r.x0, r.y0));
 512     path[1].code = NR_LINETO;
 513     path[1].setC(3, NR::Point(r.x1, r.y0));
 514     path[2].code = NR_LINETO;
 515     path[2].setC(3, NR::Point(r.x1, r.y1));
 516     path[3].code = NR_LINETO;
 517     path[3].setC(3, NR::Point(r.x0, r.y1));
 518     path[4].code = NR_LINETO;
 519     path[4].setC(3, NR::Point(r.x0, r.y0));
 520     path[5].code = NR_END;
 521
 522     return path;
 523 }
 524