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