1/*********************************************************************** 2 * * 3 * $Id: hpgsbezier.c 270 2006-01-29 21:12:23Z softadm $ 4 * * 5 * hpgs - HPGl Script, a hpgl/2 interpreter, which uses a Postscript * 6 * API for rendering a scene and thus renders to a variety of * 7 * devices and fileformats. * 8 * * 9 * (C) 2004-2006 ev-i Informationstechnologie GmbH http://www.ev-i.at * 10 * * 11 * Author: Wolfgang Glas * 12 * * 13 * hpgs is free software; you can redistribute it and/or * 14 * modify it under the terms of the GNU Lesser General Public * 15 * License as published by the Free Software Foundation; either * 16 * version 2.1 of the License, or (at your option) any later version. * 17 * * 18 * hpgs is distributed in the hope that it will be useful, * 19 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * 21 * Lesser General Public License for more details. * 22 * * 23 * You should have received a copy of the GNU Lesser General Public * 24 * License along with this library; if not, write to the * 25 * Free Software Foundation, Inc., 59 Temple Place, Suite 330, * 26 * Boston, MA 02111-1307 USA * 27 * * 28 *********************************************************************** 29 * * 30 * Bezier spline length and scanline cut implementation. * 31 * * 32 ***********************************************************************/ 33 34#include <hpgspaint.h> 35#include <math.h> 36#include <assert.h> 37 38// #define HPGS_BEZIER_DEBUG 39 40/*! Calculates the derivation of the bezier curve starting 41 at point \c i of \c path at curve parameter \c t. 42*/ 43void hpgs_bezier_path_delta(const hpgs_paint_path *path, int i, 44 double t, hpgs_point *p) 45{ 46 double tp,tm; 47 48 tp = t + 0.5; 49 tm = 0.5 - t; 50 51 p->x = 3.0* ((path->points[i+1].p.x-path->points[i].p.x)*tm*tm+ 52 2.0 * (path->points[i+2].p.x-path->points[i+1].p.x) * tm*tp+ 53 (path->points[i+3].p.x-path->points[i+2].p.x) * tp*tp); 54 55 p->y = 3.0 * ((path->points[i+1].p.y-path->points[i].p.y)*tm*tm+ 56 2.0 * (path->points[i+2].p.y-path->points[i+1].p.y) * tm*tp+ 57 (path->points[i+3].p.y-path->points[i+2].p.y) * tp*tp); 58} 59 60/*! Calculates the derivation of the bezier curve starting 61 at point \c i of \c path at curve parameter \c t. 62*/ 63double hpgs_bezier_path_delta_x(const hpgs_paint_path *path, int i, double t) 64{ 65 double tp,tm; 66 67 tp = t + 0.5; 68 tm = 0.5 - t; 69 70 return 3.0* ((path->points[i+1].p.x-path->points[i].p.x)*tm*tm+ 71 2.0 * (path->points[i+2].p.x-path->points[i+1].p.x) * tm*tp+ 72 (path->points[i+3].p.x-path->points[i+2].p.x) * tp*tp); 73} 74 75/*! Calculates the derivation of the bezier curve starting 76 at point \c i of \c path at curve parameter \c t. 77*/ 78double hpgs_bezier_path_delta_y(const hpgs_paint_path *path, int i, double t) 79{ 80 double tp,tm; 81 82 tp = t + 0.5; 83 tm = 0.5 - t; 84 85 return 3.0* ((path->points[i+1].p.y-path->points[i].p.y)*tm*tm+ 86 2.0 * (path->points[i+2].p.y-path->points[i+1].p.y) * tm*tp+ 87 (path->points[i+3].p.y-path->points[i+2].p.y) * tp*tp); 88} 89 90/*! Calculates the tangent to the bezier curve starting 91 at point \c i of \c path at curve parameter \c t. 92 93 If \c orientation is \c 1, the tangent is calculated with increasing 94 curve parameter, if we meet a singularity (right tangent). 95 96 If \c orientation is \c -1, the tangent is calculated with decreasing 97 curve parameter, if we meet a singularity (left tangent). 98 99 If \c orientation is \c 0, the tangent is calculated with decreased 100 and increased curve parameter, if we meet a singularity (mid tangent). 101*/ 102 103void hpgs_bezier_path_tangent(const hpgs_paint_path *path, int i, 104 double t, int orientation, 105 double ytol, hpgs_point *p) 106{ 107 hpgs_bezier_path_delta(path,i,t,p); 108 109 double l = hypot(p->x,p->y); 110 111 if (l < ytol) 112 { 113 hpgs_point p1; 114 115 if (orientation <= 0) 116 hpgs_bezier_path_point(path,i,t-1.0e-4,&p1); 117 else 118 hpgs_bezier_path_point(path,i,t,&p1); 119 120 if (orientation >= 0) 121 hpgs_bezier_path_point(path,i,t+1.0e-4,p); 122 else 123 hpgs_bezier_path_point(path,i,t,p); 124 125 p->x -= p1.x; 126 p->y -= p1.y; 127 128 l = hypot(p->x,p->y); 129 if (l < ytol*1.0e-4) return; 130 } 131 132 p->x /= l; 133 p->y /= l; 134} 135 136static int quad_roots (double a, double b, double c, 137 double t0, double t1, double *t) 138{ 139 if (fabs(a) < (fabs(b)>fabs(c) ? fabs(b) : fabs(c)) * 1.0e-8 || fabs(a) < 1.0e-16) 140 { 141 if (fabs(b) < fabs(c) * 1.0e-8) 142 return 0; 143 144 t[0] = -c/b; 145 return t[0] > t0 && t[0] < t1; 146 } 147 148 double p = b/a; 149 double q = c/a; 150 151 double det = 0.25*p-q; 152 153 if (det < 0) return 0; 154 155 det = sqrt(det); 156 157 if (p > 0.0) 158 t[0] = -0.5*p - det; 159 else 160 t[0] = -0.5*p + det; 161 162 if (fabs(t[0]) < 1.0e-8) 163 return t[0] > t0 && t[0] < t1; 164 165 t[1] = q / t[0]; 166 167 if (t[0] <= t0 || t[0] >= t1) 168 { 169 t[0] = t[1]; 170 return t[0] > t0 && t[0] < t1; 171 } 172 173 if (t[1] <= t0 || t[1] >= t1) 174 return 1; 175 176 if (t[0] > t[1]) 177 { 178 double tmp = t[1]; 179 t[1] = t[0]; 180 t[0] = tmp; 181 } 182 183 return fabs(t[0]-t[1]) < 1.0e-8 ? 1 : 2; 184} 185 186 187/*! 188 Calculates the singularities of the bezier curve in the 189 parameter interval from \c t0 to \c t1. 190 191 We calculate points, where the cross product of the first derivative 192 and the second derivative vanishes. These points include points, where 193 the first or second derivative vanishes. Additionally, turning points 194 of the curve are also included. 195 196 \c tx is a pointer to an array of at least two double values, 197 in which the position of the singularities is returned. 198 */ 199void hpgs_bezier_path_singularities(const hpgs_paint_path *path, int i, 200 double t0, double t1, 201 int *nx, double *tx) 202{ 203 // The second order spline p'(t) cross p''(t) 204 double k0 = 205 (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+2].p.y-path->points[i+1].p.y) - 206 (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+2].p.x-path->points[i+1].p.x); 207 208 double k1 = 209 (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) - 210 (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x); 211 212 double k2 = 213 (path->points[i+2].p.x-path->points[i+1].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) - 214 (path->points[i+2].p.y-path->points[i+1].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x); 215 216 217 // coefficients for the quadratic equations. 218 double a = k0 - k1 + k2; 219 220 double b = k2 - k0; 221 222 double c = 0.25 * (k0 + k1 + k2); 223 224 *nx = quad_roots (a,b,c,t0,t1,tx); 225} 226 227static void add_quad (const hpgs_point *p1, 228 const hpgs_point *d1, 229 const hpgs_point *p2, 230 const hpgs_point *d2, 231 int *nx, hpgs_point *points) 232{ 233 double det = d1->y * d2->x - d1->x * d2->y; 234 235 if (fabs(det) < 1.0e-8) 236 { 237 points[*nx].x = 0.5 * (p1->x + p2->x); 238 points[*nx].y = 0.5 * (p1->y + p2->y); 239 ++*nx; 240 } 241 else 242 { 243 double s = (d2->x*(p2->y-p1->y)- 244 d2->y*(p2->x-p1->x) )/det; 245 246 points[*nx].x = p1->x + s*d1->x; 247 points[*nx].y = p1->y + s*d1->y; 248 ++*nx; 249 } 250 251 points[*nx].x = p2->x; 252 points[*nx].y = p2->y; 253 ++*nx; 254} 255 256/*! 257 Approximates the bezier curve segment in the 258 parameter interval from \c t0 to \c t1 with quadratic bezier splines. 259 260 The array \c points is used to return the quadratic bezier segments and 261 must have a dimension of at least 16 points. 262 263 \c nx returns the number of generated points, which is always an even 264 number. The first point is the control point of the first quadratic 265 spline, the second point the endpoint of the first quadratic spline 266 and so on. 267 */ 268void hpgs_bezier_path_to_quadratic(const hpgs_paint_path *path, int i, 269 double t0, double t1, 270 int *nx, hpgs_point *points) 271{ 272 // The second order spline p'(t) cross p''(t) 273 double k0 = 274 (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+2].p.y-path->points[i+1].p.y) - 275 (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+2].p.x-path->points[i+1].p.x); 276 277 double k1 = 278 (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) - 279 (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x); 280 281 double k2 = 282 (path->points[i+2].p.x-path->points[i+1].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) - 283 (path->points[i+2].p.y-path->points[i+1].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x); 284 285 286 double xmax,xmin,ymin,ymax; 287 288 // coefficients for the quadratic equations. 289 double a = k0 - k1 + k2; 290 291 double b = k2 - k0; 292 293 double c = 0.25 * (k0 + k1 + k2); 294 295 // at most three partition points: 296 // extrem, roots of the curvature spline 297 double tpart[5]; 298 299 // get the roots. 300 int ipart; 301 int npart = quad_roots (a,b,c,t0,t1,tpart+1) + 1; 302 303 // add the extreme point. 304 if (fabs(a) > fabs(b)) 305 { 306 double text = -0.5 * b/a; 307 308 if (npart == 3) 309 { 310 tpart[3] = tpart[2]; 311 tpart[2] = text; 312 npart = 3; 313 } 314 else if (text > t0 && text < t1) 315 { 316 if (npart == 2 && text < tpart[1]) 317 { 318 tpart[2] = tpart[1]; 319 tpart[1] = text; 320 } 321 else 322 tpart[npart] = text; 323 324 ++npart; 325 } 326 } 327 328 tpart[0] = t0; 329 tpart[npart] = t1; 330 331 xmax = xmin = path->points[i].p.x; 332 if (path->points[i+1].p.x > xmax) xmax = path->points[i+1].p.x; 333 if (path->points[i+2].p.x > xmax) xmax = path->points[i+2].p.x; 334 if (path->points[i+3].p.x > xmax) xmax = path->points[i+3].p.x; 335 if (path->points[i+1].p.x < xmin) xmin = path->points[i+1].p.x; 336 if (path->points[i+2].p.x < xmin) xmin = path->points[i+2].p.x; 337 if (path->points[i+3].p.x < xmin) xmin = path->points[i+3].p.x; 338 339 ymax = ymin = path->points[i].p.y; 340 if (path->points[i+1].p.y > ymax) ymax = path->points[i+1].p.y; 341 if (path->points[i+2].p.y > ymax) ymax = path->points[i+2].p.y; 342 if (path->points[i+3].p.y > ymax) ymax = path->points[i+3].p.y; 343 if (path->points[i+1].p.y < ymin) ymin = path->points[i+1].p.y; 344 if (path->points[i+2].p.y < ymin) ymin = path->points[i+2].p.y; 345 if (path->points[i+3].p.y < ymin) ymin = path->points[i+3].p.y; 346 347 double curv0 = 348 k0*(0.5-tpart[0])*(0.5-tpart[0]) + 349 k1*(0.5-tpart[0])*(0.5+tpart[0]) + 350 k2*(0.5+tpart[0])*(0.5+tpart[0]); 351 352 *nx = 0; 353 354 for (ipart=0; ipart<npart; ++ipart) 355 { 356 double curv1 = 357 k0*(0.5-tpart[ipart+1])*(0.5-tpart[ipart+1]) + 358 k1*(0.5-tpart[ipart+1])*(0.5+tpart[ipart+1]) + 359 k2*(0.5+tpart[ipart+1])*(0.5+tpart[ipart+1]); 360 361 int n_splines = 1; 362 363 if (fabs(xmax-xmin)<1.0e-8 && fabs(ymax-ymin) < 1.0e-8) 364 n_splines = (int)(16.0 * fabs(curv0-curv1)/((xmax-xmin)*(ymax-ymin))); 365 366 int j; 367 368 if (n_splines < 1) n_splines = 1; 369 else if (n_splines > 8/npart) n_splines = 8/npart; 370 371 hpgs_point d1; 372 hpgs_point p1; 373 374 hpgs_bezier_path_point(path,i,tpart[ipart],&p1); 375 hpgs_bezier_path_tangent(path,i,tpart[ipart],1,1.0e-3,&d1); 376 377 for (j=1;j<=n_splines;++j) 378 { 379 double t = tpart[ipart] + (tpart[ipart+1]-tpart[ipart]) * j / (double)n_splines; 380 381 hpgs_point d2; 382 hpgs_point p2; 383 384 hpgs_bezier_path_point(path,i,t,&p2); 385 hpgs_bezier_path_tangent(path,i,t,j<n_splines ? 0 : -1,1.0e-3,&d2); 386 387 add_quad(&p1,&d1,&p2,&d2,nx,points); 388 389 p1 = p2; 390 d1 = d2; 391 } 392 393 curv0 = curv1; 394 } 395} 396 397/* Internal: Calculates the derivation of the curve length 398 of the bezier curve \c b at curve parameter \c t. 399*/ 400static double bezier_dl(const hpgs_paint_path *path, int i, double t) 401{ 402 double tp,tm; 403 404 tp = t + 0.5; 405 tm = 0.5 - t; 406 407 return 3.0*hypot((path->points[i+1].p.x-path->points[i].p.x)*tm*tm+ 408 2.0 * (path->points[i+2].p.x-path->points[i+1].p.x) * tm*tp+ 409 (path->points[i+3].p.x-path->points[i+2].p.x) * tp*tp, 410 411 (path->points[i+1].p.y-path->points[i].p.y)*tm*tm+ 412 2.0 * (path->points[i+2].p.y-path->points[i+1].p.y) * tm*tp+ 413 (path->points[i+3].p.y-path->points[i+2].p.y) * tp*tp 414 ); 415} 416 417/*! Calculates the curve lengths at an equidistant grid 418 of curve parmeters by na numerical quadrature. 419 420 You must call this subroutine before using 421 \c hpgs_bezier_length_param. 422*/ 423void hpgs_bezier_length_init(hpgs_bezier_length *b, 424 const hpgs_paint_path *path, int i) 425{ 426 int ip; 427 428 b->dl[0] = bezier_dl(path,i,-0.5); 429 430 b->l[0] = 0.0; 431 432 for (ip=1;ip<=16;++ip) 433 { 434 double dlm0,dlm1; 435 436 // Lobatto quadrature with order 4. 437 b->dl[ip] = bezier_dl(path,i,(ip-8.0)*0.0625); 438 439 dlm0 = bezier_dl(path,i,(ip-(8.5+0.22360679774997896964))*0.0625); 440 dlm1 = bezier_dl(path,i,(ip-(8.5-0.22360679774997896964))*0.0625); 441 442 b->l[ip] = b->l[ip-1] + 443 (b->dl[ip-1] + 5.0 * (dlm0 + dlm1) + b->dl[ip])/(16.0*12.0); 444 } 445} 446 447/*! Returns the curve parameter at the bezier curve \c b 448 corresponding to a curve length \c l. 449 450 You must have called \c hpgs_bezier_length_init before using 451 this function. 452*/ 453double hpgs_bezier_length_param(const hpgs_bezier_length *b, double l) 454{ 455 // binary search. 456 457 int i = 0; 458 459 if (l >= b->l[i+8]) i+=8; 460 if (l >= b->l[i+4]) i+=4; 461 if (l >= b->l[i+2]) i+=2; 462 if (l >= b->l[i+1]) ++i; 463 464 double dl = b->l[i+1]-b->l[i]; 465 double ll = (l-b->l[i])/dl; 466 467 return 468 (i - 8.0 + ll*ll*(3.0-2.0*ll)) * 0.0625 + 469 (ll*(ll*(ll-2.0)+1.0) / b->dl[i] + ll*ll*(ll-1.0) / b->dl[i+1] ) * dl; 470} 471 472/*! Calculates the minimal x component of all four control points of 473 the bezier spline. */ 474double hpgs_bezier_path_xmin (const hpgs_paint_path *path, int i) 475{ 476 double x0 = 477 path->points[i].p.x < path->points[i+1].p.x ? 478 path->points[i].p.x : path->points[i+1].p.x; 479 480 double x1 = 481 path->points[i+2].p.x < path->points[i+3].p.x ? 482 path->points[i+2].p.x : path->points[i+3].p.x; 483 484 return x0 < x1 ? x0 : x1; 485} 486 487/*! Calculates the maximal x component of all four control points of 488 the bezier spline. */ 489double hpgs_bezier_path_xmax (const hpgs_paint_path *path, int i) 490{ 491 double x0 = 492 path->points[i].p.x > path->points[i+1].p.x ? 493 path->points[i].p.x : path->points[i+1].p.x; 494 495 double x1 = 496 path->points[i+2].p.x > path->points[i+3].p.x ? 497 path->points[i+2].p.x : path->points[i+3].p.x; 498 499 return x0 > x1 ? x0 : x1; 500} 501 502/*! Calculates the minimal y component of all four control points of 503 the bezier spline. */ 504double hpgs_bezier_path_ymin (const hpgs_paint_path *path, int i) 505{ 506 double y0 = 507 path->points[i].p.y < path->points[i+1].p.y ? 508 path->points[i].p.y : path->points[i+1].p.y; 509 510 double y1 = 511 path->points[i+2].p.y < path->points[i+3].p.y ? 512 path->points[i+2].p.y : path->points[i+3].p.y; 513 514 return y0 < y1 ? y0 : y1; 515} 516 517/*! Calculates the maximal y component of all four control points of 518 the bezier spline. */ 519double hpgs_bezier_path_ymax (const hpgs_paint_path *path, int i) 520{ 521 double y0 = 522 path->points[i].p.y > path->points[i+1].p.y ? 523 path->points[i].p.y : path->points[i+1].p.y; 524 525 double y1 = 526 path->points[i+2].p.y > path->points[i+3].p.y ? 527 path->points[i+2].p.y : path->points[i+3].p.y; 528 529 return y0 > y1 ? y0 : y1; 530} 531 532/*! Calculates the x component of a point on the bezier spline. */ 533double hpgs_bezier_path_x (const hpgs_paint_path *path, int i, double t) 534{ 535 double tp,tm; 536 537 tp = t + 0.5; 538 tm = 0.5 - t; 539 540 return 541 path->points[i].p.x *tm*tm*tm + 542 3.0 * tm*tp * (path->points[i+1].p.x * tm + path->points[i+2].p.x * tp) + 543 path->points[i+3].p.x * tp*tp*tp; 544} 545 546/*! Calculates the y component of a point on the bezier spline. */ 547double hpgs_bezier_path_y (const hpgs_paint_path *path, int i, double t) 548{ 549 double tp,tm; 550 551 tp = t + 0.5; 552 tm = 0.5 - t; 553 554 return 555 path->points[i].p.y *tm*tm*tm + 556 3.0 * tm*tp * (path->points[i+1].p.y * tm + path->points[i+2].p.y * tp) + 557 path->points[i+3].p.y * tp*tp*tp; 558} 559 560/*! Calculates the point on the bezier curve starting 561 at point \c i of \c path at curve parameter \c t. 562*/ 563void hpgs_bezier_path_point(const hpgs_paint_path *path, int i, 564 double t, hpgs_point *p) 565{ 566 double tp,tm; 567 568 tp = t + 0.5; 569 tm = 0.5 - t; 570 571 p->x = 572 path->points[i].p.x *tm*tm*tm + 573 3.0 * tm*tp * (path->points[i+1].p.x * tm + path->points[i+2].p.x * tp) + 574 path->points[i+3].p.x * tp*tp*tp; 575 576 p->y = 577 path->points[i].p.y *tm*tm*tm + 578 3.0 * tm*tp * (path->points[i+1].p.y * tm + path->points[i+2].p.y * tp) + 579 path->points[i+3].p.y * tp*tp*tp; 580} 581