1 /* graph.c 1.12 84/03/27 2 * 3 * This file contains the functions for producing the graphics 4 * images in the varian/versatec drivers for ditroff. 5 */ 6 7 8 #include <stdio.h> 9 #include <ctype.h> 10 #include <math.h> 11 12 13 #define TRUE 1 14 #define FALSE 0 15 #define FATAL 1 16 /* imports from dver.c */ 17 #define hmot(n) hpos += n; 18 #define vmot(n) vgoto(vpos + n); 19 20 extern int hpos; 21 extern int vpos; 22 extern int output; 23 extern vgoto(); 24 extern point(); 25 26 #define MAXPOINTS 200 /* number of points legal for a curve */ 27 28 #define SOLID -1 /* line styles: these are used as bit masks to */ 29 #define DOTTED 004 /* create the right style lines. */ 30 #define DASHED 020 31 #define DOTDASHED 024 32 #define LONGDASHED 074 33 /* constants... */ 34 #define pi 3.14159265358979324 35 #define log2_10 3.3219280948873623 36 37 38 int linethickness = 3; /* number of pixels wide to make lines */ 39 int linmod = SOLID; /* type of line (SOLID, DOTTED, DASHED...) */ 40 41 42 43 /*---------------------------------------------------------------------------- 44 * Routine: drawline (horizontal_motion, vertical_motion) 45 * 46 * Results: Draws a line of "linethickness" width and "linmod" style 47 * from current (hpos, vpos) to (hpos + dh, vpos + dv). 48 * 49 * Side Efct: Resulting position is at end of line (hpos + dh, vpos + dv) 50 *----------------------------------------------------------------------------*/ 51 52 drawline(dh, dv) 53 register int dh; 54 register int dv; 55 { 56 if (output) HGtline (hpos, vpos, hpos + dh, vpos + dv); 57 hmot (dh); /* new position is at */ 58 vmot (dv); /* the end of the line */ 59 } 60 61 62 /*---------------------------------------------------------------------------- 63 * Routine: drawcirc (diameter) 64 * 65 * Results: Draws a circle with leftmost point at current (hpos, vpos) 66 * with the given diameter d. 67 * 68 * Side Efct: Resulting position is at (hpos + diameter, vpos) 69 *----------------------------------------------------------------------------*/ 70 71 drawcirc(d) 72 register int d; 73 { /* 0.0 is the angle to sweep the arc: = full circle */ 74 if (output) HGArc (hpos + d/2, vpos, hpos, vpos, 0.0); 75 hmot (d); /* new postion is the right of the circle */ 76 } 77 78 79 /*---------------------------------------------------------------------------- 80 * Routine: drawellip (horizontal_diameter, vertical_diameter) 81 * 82 * Results: Draws regular ellipses given the major "diameters." It does 83 * so using a modified circle algorithm (see RoundEnd) that 84 * increments x and y proportionally to their axes. 85 * 86 * Side Efct: Resulting position is at (hpos + hd, vpos). 87 *----------------------------------------------------------------------------*/ 88 89 drawellip(hd, vd) 90 int hd; 91 int vd; 92 { 93 double xs, ys, xepsilon, yepsilon; 94 register int thick; 95 register int basex; 96 register int basey; 97 register int x; 98 register int y; 99 100 101 basex = hpos + (hd >> 1); /* bases == coordinates of center */ 102 hmot(hd); /* horizontal motion (hpos should */ 103 basey = vpos; /* NOT be used after this) */ 104 /* hd and vd are radii, not diameters */ 105 if ((hd = hd >> 1) < 1) hd++; /* neither diameter can be zero. */ 106 if ((vd = vd >> 1) < 1) vd++; /* hd changed!! no hmot after this */ 107 ys = (double) vd; /* start at top of the ellipse */ 108 xs = 0.0; /* (y = 1/2 diameter, x = 0) */ 109 110 if ((thick = vd) > hd) thick = hd; 111 xepsilon = (double) thick / (double) (vd * vd); 112 yepsilon = (double) thick / (double) (hd * hd); 113 114 /* Calculate trajectory of the ellipse for 1/4 */ 115 /* the circumference (while ys is non-negative) */ 116 /* and mirror to get the other three quadrants. */ 117 118 if (linethickness > 1) { /* more than one pixel thick */ 119 RoundEnd(basex, ((int)(ys + 0.5)) + basey, linethickness, 0); 120 RoundEnd(basex, basey - ((int)(ys + 0.5)), linethickness, 0); 121 while (ys >= 0) { 122 xs += xepsilon * ys; /* generate circumference */ 123 ys -= yepsilon * xs; 124 x = (int)(xs + 0.5); 125 y = (int)(ys + 0.5); 126 RoundEnd(x + basex, y + basey, linethickness, 0); 127 RoundEnd(x + basex, basey - y, linethickness, 0); 128 RoundEnd(basex - x, y + basey, linethickness, 0); 129 RoundEnd(basex - x, basey - y, linethickness, 0); 130 } 131 } else { /* do the perimeter only (no fill) */ 132 point(basex, ((int)(ys + 0.5)) + basey); 133 point(basex, basey - ((int)(ys + 0.5))); 134 while (ys >= 0) { 135 xs += xepsilon * ys; /* generate circumference */ 136 ys -= yepsilon * xs; 137 x = (int)(xs + 0.5); 138 y = (int)(ys + 0.5); 139 point(x + basex, y + basey); 140 point(x + basex, basey - y); 141 point(basex - x, y + basey); 142 point(basex - x, basey - y); 143 } 144 } 145 } 146 147 148 /*---------------------------------------------------------------------------- 149 * Routine: drawarc (xcenter, ycenter, xpoint, ypoint) 150 * 151 * Results: Draws an arc starting at current (hpos, vpos). Center is 152 * at (hpos + cdh, vpos + cdv) and the terminating point is 153 * at <center> + (pdh, pdv). The angle between the lines 154 * formed by the starting, ending, and center points is figured 155 * first, then the points and angle are sent to HGArc for the 156 * drawing. 157 * 158 * Side Efct: Resulting position is at the last point of the arc. 159 *----------------------------------------------------------------------------*/ 160 161 drawarc (cdh, cdv, pdh, pdv) 162 register int cdh; 163 register int cdv; 164 register int pdh; 165 register int pdv; 166 { 167 register double angle; 168 /* figure angle from the three points...*/ 169 /* and convert (and round) to degrees */ 170 angle = (atan2((double) pdh, (double) pdv) 171 - atan2 ((double) -cdh, (double) -cdv)) * 180.0 / pi; 172 /* "normalize" and round */ 173 angle += (angle < 0.0) ? 360.5 : 0.5; 174 175 if (output) HGArc(hpos + cdh, vpos + cdv, hpos, vpos, (int) angle); 176 hmot(cdh + pdh); 177 vmot(cdv + pdv); 178 } 179 180 181 /*---------------------------------------------------------------------------- 182 * Routine: drawwig (character_buffer, file_pointer, type_flag) 183 * 184 * Results: Given the starting position, the motion list in buf, and any 185 * extra characters from fp (terminated by a \n), drawwig sets 186 * up a point list to make a spline from. If "pic" is set picurve 187 * is called to draw the curve in pic style; else it calls HGCurve 188 * for the gremlin-style curve. 189 * 190 * Side Efct: Resulting position is reached from adding successive motions 191 * to the current position. 192 *----------------------------------------------------------------------------*/ 193 194 drawwig (buf, fp, pic) 195 char *buf; 196 FILE *fp; 197 int pic; 198 { 199 register int len = strlen(buf); /* length of the string in "buf" */ 200 register int npts = 2; /* point list index */ 201 register char *ptr = buf; /* "walking" pointer into buf */ 202 int x[MAXPOINTS], y[MAXPOINTS]; /* point list */ 203 204 while (*ptr == ' ') ptr++; /* skip any leading spaces */ 205 x[1] = hpos; /* the curve starts at the */ 206 y[1] = vpos; /* current position */ 207 208 while (*ptr != '\n') { /* curve commands end with a '\n' */ 209 hmot(atoi(ptr)); /* convert motion to curve points */ 210 x[npts] = hpos; /* and remember them */ 211 while (isdigit(*++ptr)); /* skip number*/ 212 while (*++ptr == ' '); /* skip spaces 'tween numbers */ 213 vmot(atoi(ptr)); 214 y[npts] = vpos; 215 while (isdigit(*++ptr)); 216 while (*ptr == ' ') ptr++; 217 /* if the amount we read wasn't the */ 218 /* whole thing, read some more in */ 219 if (len - (ptr - buf) < 15 && *(buf + len - 1) != '\n') { 220 char *cop = buf; 221 222 while (*cop++ = *ptr++); /* copy what's left to the beginning */ 223 if (fgets ((cop - 1), len - (cop - buf), fp) == NULL) 224 error(FATAL, "unexpected end of input"); 225 ptr = buf; 226 } 227 if (npts < MAXPOINTS - 1) /* if too many points, forget some */ 228 npts++; 229 } 230 npts--; /* npts must point to the last coordinate in x and y */ 231 /* now, actually DO the curve */ 232 if (output) { 233 if (pic > 0) 234 picurve(x, y, npts); 235 else if (pic < 0) 236 polygon(x, y, npts); 237 else 238 HGCurve(x, y, npts); 239 } 240 } 241 242 243 /*----------------------------------------------------------------------------* 244 | Routine: drawthick (thickness) 245 | 246 | Results: sets the variable "linethickness" to the given size. 247 | NO motion is involved. 248 *----------------------------------------------------------------------------*/ 249 250 drawthick(s) 251 int s; 252 { 253 linethickness = s; 254 } 255 256 257 /*----------------------------------------------------------------------------* 258 | Routine: drawstyle (style_bit_map) 259 | 260 | Results: sets the variable "linmod" to the given bit map. 261 | NO motion is involved. 262 *----------------------------------------------------------------------------*/ 263 264 drawstyle(s) 265 int s; 266 { 267 linmod = s; 268 } 269 270 271 /*---------------------------------------------------------------------------- 272 * Routine: picurve (xpoints, ypoints, num_of_points) 273 * 274 * Results: Draws a curve delimited by (not through) the line segments 275 * traced by (xpoints, ypoints) point list. This is the "Pic" 276 * style curve. 277 *----------------------------------------------------------------------------*/ 278 279 picurve (x, y, npts) 280 int x[MAXPOINTS]; 281 int y[MAXPOINTS]; 282 int npts; 283 { 284 register int i; /* line segment traverser */ 285 register int nseg; /* effective resolution for each curve */ 286 register float w; /* position factor */ 287 register int xp; /* current point (and intermediary) */ 288 register int yp; 289 int pxp, pyp; /* previous point (to make lines from) */ 290 float t1, t2, t3; /* calculation temps */ 291 int j; /* inner curve segment traverser */ 292 293 294 if (x[1] == x[npts] && y[1] == y[npts]) { 295 x[0] = x[npts - 1]; /* if the lines' ends meet, make */ 296 y[0] = y[npts - 1]; /* sure the curve meets */ 297 x[npts + 1] = x[2]; 298 y[npts + 1] = y[2]; 299 } else { /* otherwise, make the ends of the */ 300 x[0] = x[1]; /* curve touch the ending points of */ 301 y[0] = y[1]; /* the line segments */ 302 x[npts + 1] = x[npts]; 303 y[npts + 1] = y[npts]; 304 } 305 306 pxp = (x[0] + x[1]) / 2; /* make the last point pointers */ 307 pyp = (y[0] + y[1]) / 2; /* point to the start of the 1st line */ 308 309 for (i = 0; i < npts; i++) { /* traverse the line segments */ 310 xp = x[i] - x[i+1]; 311 yp = y[i] - y[i+1]; 312 nseg = (int) sqrt((double)(xp * xp + yp * yp)); 313 xp = x[i+1] - x[i+2]; 314 yp = y[i+1] - y[i+2]; /* "nseg" is the number of line */ 315 /* segments that will be drawn for */ 316 /* each curve segment. ">> 3" is */ 317 /* dropping the resolution ( == / 8) */ 318 nseg = (nseg + (int) sqrt((double)(xp * xp + yp * yp))) >> 3; 319 320 for (j = 1; j < nseg; j++) { 321 w = (float) j / (float) nseg; 322 t1 = 0.5 * w * w; 323 w -= 0.5; 324 t2 = 0.75 - w * w ; 325 w -= 0.5; 326 t3 = 0.5 * w * w; 327 xp = t1 * x[i+2] + t2 * x[i+1] + t3 * x[i] + 0.5; 328 yp = t1 * y[i+2] + t2 * y[i+1] + t3 * y[i] + 0.5; 329 330 HGtline(pxp, pyp, xp, yp); 331 pxp = xp; 332 pyp = yp; 333 } 334 } 335 } 336 337 338 /*---------------------------------------------------------------------------- 339 * Routine: line (xstart, ystart, xend, yend) 340 * 341 * Results: Draws a one-pixel wide line from (x0, y0) to (x1, y1) 342 * using point(x,y) to place the dots. Line style 343 * is done in this routine by masking off unwanted dots. 344 *----------------------------------------------------------------------------*/ 345 346 line(x0, y0, x1, y1) 347 register int x0; 348 register int y0; 349 int x1; 350 int y1; 351 { 352 register int dx; 353 register int dy; 354 register int res; 355 register int slope; 356 int xinc; 357 int yinc; 358 359 xinc = 1; 360 yinc = 1; 361 if ((dx = x1-x0) < 0) { 362 xinc = -1; 363 dx = -dx; 364 } 365 if ((dy = y1-y0) < 0) { 366 yinc = -1; 367 dy = -dy; 368 } 369 slope = xinc*yinc; 370 if (dx >= dy) { 371 res = (dy >> 1) - (dx >> 1); 372 while (x0 != x1) { 373 if ((x0+slope*y0)&linmod) point(x0, y0); 374 if (res >= 0) { 375 res -= dx; 376 y0 += yinc; 377 } 378 res += dy; 379 x0 += xinc; 380 } 381 } else { 382 res = (dx >> 1) - (dy >> 1); 383 while (y0 != y1) { 384 if ((x0+slope*y0)&linmod) point(x0, y0); 385 if (res >= 0) { 386 res -= dy; 387 x0 += xinc; 388 } 389 res += dx; 390 y0 += yinc; 391 } 392 } 393 if((x1+slope*y1)&linmod) point(x1, y1); 394 } 395 396 397 /*---------------------------------------------------------------------------- 398 * Routine: HGArc (xcenter, ycenter, xstart, ystart, angle) 399 * 400 * Results: This routine plots an arc centered about (cx, cy) counter 401 * clockwise starting from the point (px, py) through 'angle' 402 * degrees. If angle is 0, a full circle is drawn. 403 * It does so by calling RoundEnd (fat point maker) for points 404 * along the circle with density depending on the circle's size. 405 *----------------------------------------------------------------------------*/ 406 407 HGArc(cx,cy,px,py,angle) 408 register int cx; 409 register int cy; 410 int px; 411 int py; 412 int angle; 413 { 414 double xs, ys, resolution, epsilon, fullcircle; 415 register int nx; 416 register int ny; 417 register int extent; 418 419 xs = px - cx; 420 ys = py - cy; 421 422 /* calculate drawing parameters */ 423 424 resolution = log10(sqrt( xs * xs + ys * ys)) * log2_10; 425 resolution = ceil(resolution); 426 resolution = pow(2.0, resolution); 427 epsilon = 1.0 / resolution; 428 fullcircle = 2 * pi * resolution; 429 fullcircle = ceil(fullcircle); 430 431 if (angle == 0) 432 extent = fullcircle; 433 else 434 extent = angle * fullcircle / 360.0; 435 436 while (extent-- > 0) { 437 xs += epsilon * ys; 438 nx = cx + (int) (xs + 0.5); 439 ys -= epsilon * xs; 440 ny = cy + (int) (ys + 0.5); 441 RoundEnd(nx, ny, linethickness, FALSE); 442 } /* end for */ 443 } /* end HGArc */ 444 445 446 /*---------------------------------------------------------------------------- 447 * Routine: RoundEnd (x, y, diameter, filled_flag) 448 * 449 * Results: Plots a filled (if requested) circle of the specified diameter 450 * centered about (x, y). 451 *----------------------------------------------------------------------------*/ 452 453 RoundEnd(x, y, diameter, filled) 454 register int x; 455 register int y; 456 int diameter; 457 int filled; 458 { 459 double xs, ys; /* floating point distance form center of circle */ 460 double epsilon; /* "resolution" of the step around circle */ 461 register int cy; /* to index up from center of circle */ 462 register int nx; /* integer distance from center */ 463 register int ny; 464 register int arc; /* counts how far around the circle to go */ 465 466 467 if (diameter < 2) { /* too small to notice */ 468 point(x, y); 469 return; 470 } 471 472 xs = 0; 473 ys = (double) (diameter - 1) / 2.0; 474 epsilon = 1.0 / ys; 475 arc = (pi / 2.0) * ys; 476 if (arc < 4) { /* if too small, make it bigger */ 477 arc += arc; /* to try and fill in more. */ 478 epsilon /= 2.0; 479 } 480 481 /* Calculate the trajectory of the circle for 1/4 the circumference 482 * and mirror appropriately to get the other three quadrants. 483 */ 484 485 nx = 0; /* must start out the x and y for first */ 486 ny = (int) (ys + 0.5); /* painting in while loop */ 487 488 while (arc-- >= 0) { 489 if (filled) { /* fill from center */ 490 cy = 0; 491 } else { /* fill from perimeter only (no fill) */ 492 cy = ny; 493 } 494 while (cy <= ny) { 495 point(nx + x, cy + y); 496 point(nx + x, y - cy); 497 point(x - nx, cy + y); 498 point(x - nx, y - cy); 499 cy++; 500 } /* end while cy */ 501 /* generate circumference */ 502 nx = (int) ((xs += epsilon * ys) + 0.5); 503 ny = (int) ((ys -= epsilon * xs) + 0.5); 504 } /* end while arc */; 505 } /* end RoundEnd */; 506 507 508 /*---------------------------------------------------------------------------- 509 * Routine: Paramaterize (xpoints, ypoints, hparams, num_points) 510 * 511 * Results: This routine calculates parameteric values for use in 512 * calculating curves. The parametric values are returned 513 * in the array h. The values are an approximation of 514 * cumulative arc lengths of the curve (uses cord length). 515 * For additional information, see paper cited below. 516 *----------------------------------------------------------------------------*/ 517 518 static Paramaterize(x, y, h, n) 519 int x[MAXPOINTS]; 520 int y[MAXPOINTS]; 521 float h[MAXPOINTS]; 522 int n; 523 { 524 register int dx; 525 register int dy; 526 register int i; 527 register int j; 528 float u[MAXPOINTS]; 529 530 531 for (i=1; i<=n; ++i) { 532 u[i] = 0; 533 for (j=1; j<i; j++) { 534 dx = x[j+1] - x[j]; 535 dy = y[j+1] - y[j]; 536 u[i] += sqrt ((double) (dx * dx + dy * dy)); 537 } 538 } 539 for (i=1; i<n; ++i) h[i] = u[i+1] - u[i]; 540 } /* end Paramaterize */ 541 542 543 /*---------------------------------------------------------------------------- 544 * Routine: PeriodicSpline (h, z, dz, d2z, d3z, npoints) 545 * 546 * Results: This routine solves for the cubic polynomial to fit a 547 * spline curve to the the points specified by the list 548 * of values. The Curve generated is periodic. The algorithms 549 * for this curve are from the "Spline Curve Techniques" paper 550 * cited below. 551 *----------------------------------------------------------------------------*/ 552 553 static PeriodicSpline(h, z, dz, d2z, d3z, npoints) 554 float h[MAXPOINTS]; /* paramaterization */ 555 int z[MAXPOINTS]; /* point list */ 556 float dz[MAXPOINTS]; /* to return the 1st derivative */ 557 float d2z[MAXPOINTS], d3z[MAXPOINTS]; /* 2nd and 3rd derivatives */ 558 int npoints; /* number of valid points */ 559 { 560 float d[MAXPOINTS]; 561 float deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS]; 562 float c[MAXPOINTS], r[MAXPOINTS], s[MAXPOINTS]; 563 int i; 564 565 /* step 1 */ 566 for (i=1; i<npoints; ++i) { 567 deltaz[i] = h[i] ? ((double) (z[i+1] - z[i])) / h[i] : 0; 568 } 569 h[0] = h[npoints-1]; 570 deltaz[0] = deltaz[npoints-1]; 571 572 /* step 2 */ 573 for (i=1; i<npoints-1; ++i) { 574 d[i] = deltaz[i+1] - deltaz[i]; 575 } 576 d[0] = deltaz[1] - deltaz[0]; 577 578 /* step 3a */ 579 a[1] = 2 * (h[0] + h[1]); 580 b[1] = d[0]; 581 c[1] = h[0]; 582 for (i=2; i<npoints-1; ++i) { 583 a[i] = 2*(h[i-1]+h[i]) - pow ((double) h[i-1],(double)2.0) / a[i-1]; 584 b[i] = d[i-1] - h[i-1] * b[i-1]/a[i-1]; 585 c[i] = -h[i-1] * c[i-1]/a[i-1]; 586 } 587 588 /* step 3b */ 589 r[npoints-1] = 1; 590 s[npoints-1] = 0; 591 for (i=npoints-2; i>0; --i) { 592 r[i] = -(h[i] * r[i+1] + c[i])/a[i]; 593 s[i] = (6 * b[i] - h[i] * s[i+1])/a[i]; 594 } 595 596 /* step 4 */ 597 d2z[npoints-1] = (6 * d[npoints-2] - h[0] * s[1] 598 - h[npoints-1] * s[npoints-2]) 599 / (h[0] * r[1] + h[npoints-1] * r[npoints-2] 600 + 2 * (h[npoints-2] + h[0])); 601 for (i=1; i<npoints-1; ++i) { 602 d2z[i] = r[i] * d2z[npoints-1] + s[i]; 603 } 604 d2z[npoints] = d2z[1]; 605 606 /* step 5 */ 607 for (i=1; i<npoints; ++i) { 608 dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i+1])/6; 609 d3z[i] = h[i] ? (d2z[i+1] - d2z[i])/h[i] : 0; 610 } 611 } /* end PeriodicSpline */ 612 613 614 /*---------------------------------------------------------------------------- 615 * Routine: NaturalEndSpline (h, z, dz, d2z, d3z, npoints) 616 * 617 * Results: This routine solves for the cubic polynomial to fit a 618 * spline curve the the points specified by the list of 619 * values. The alogrithms for this curve are from the 620 * "Spline Curve Techniques" paper cited below. 621 *----------------------------------------------------------------------------*/ 622 623 static NaturalEndSpline(h, z, dz, d2z, d3z, npoints) 624 float h[MAXPOINTS]; /* parameterization */ 625 int z[MAXPOINTS]; /* Point list */ 626 float dz[MAXPOINTS]; /* to return the 1st derivative */ 627 float d2z[MAXPOINTS], d3z[MAXPOINTS]; /* 2nd and 3rd derivatives */ 628 int npoints; /* number of valid points */ 629 { 630 float d[MAXPOINTS]; 631 float deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS]; 632 int i; 633 634 /* step 1 */ 635 for (i=1; i<npoints; ++i) { 636 deltaz[i] = h[i] ? ((double) (z[i+1] - z[i])) / h[i] : 0; 637 } 638 deltaz[0] = deltaz[npoints-1]; 639 640 /* step 2 */ 641 for (i=1; i<npoints-1; ++i) { 642 d[i] = deltaz[i+1] - deltaz[i]; 643 } 644 d[0] = deltaz[1] - deltaz[0]; 645 646 /* step 3 */ 647 a[0] = 2 * (h[2] + h[1]); 648 b[0] = d[1]; 649 for (i=1; i<npoints-2; ++i) { 650 a[i] = 2*(h[i+1]+h[i+2]) - pow((double) h[i+1],(double) 2.0)/a[i-1]; 651 b[i] = d[i+1] - h[i+1] * b[i-1]/a[i-1]; 652 } 653 654 /* step 4 */ 655 d2z[npoints] = d2z[1] = 0; 656 for (i=npoints-1; i>1; --i) { 657 d2z[i] = (6 * b[i-2] - h[i] *d2z[i+1])/a[i-2]; 658 } 659 660 /* step 5 */ 661 for (i=1; i<npoints; ++i) { 662 dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i+1])/6; 663 d3z[i] = h[i] ? (d2z[i+1] - d2z[i])/h[i] : 0; 664 } 665 } /* end NaturalEndSpline */ 666 667 668 /*---------------------------------------------------------------------------- 669 * Routine: HGCurve(xpoints, ypoints, num_points) 670 * 671 * Results: This routine generates a smooth curve through a set of points. 672 * The method used is the parametric spline curve on unit knot 673 * mesh described in "Spline Curve Techniques" by Patrick 674 * Baudelaire, Robert Flegal, and Robert Sproull -- Xerox Parc. 675 *----------------------------------------------------------------------------*/ 676 677 #define PointsPerInterval 32 678 679 HGCurve(x, y, numpoints) 680 int x[MAXPOINTS]; 681 int y[MAXPOINTS]; 682 int numpoints; 683 { 684 float h[MAXPOINTS], dx[MAXPOINTS], dy[MAXPOINTS]; 685 float d2x[MAXPOINTS], d2y[MAXPOINTS], d3x[MAXPOINTS], d3y[MAXPOINTS]; 686 float t, t2, t3; 687 register int j; 688 register int k; 689 register int nx; 690 register int ny; 691 int lx, ly; 692 693 694 lx = x[1]; 695 ly = y[1]; 696 697 /* Solve for derivatives of the curve at each point 698 * separately for x and y (parametric). 699 */ 700 Paramaterize(x, y, h, numpoints); 701 /* closed curve */ 702 if ((x[1] == x[numpoints]) && (y[1] == y[numpoints])) { 703 PeriodicSpline(h, x, dx, d2x, d3x, numpoints); 704 PeriodicSpline(h, y, dy, d2y, d3y, numpoints); 705 } else { 706 NaturalEndSpline(h, x, dx, d2x, d3x, numpoints); 707 NaturalEndSpline(h, y, dy, d2y, d3y, numpoints); 708 } 709 710 /* generate the curve using the above information and 711 * PointsPerInterval vectors between each specified knot. 712 */ 713 714 for (j=1; j<numpoints; ++j) { 715 if ((x[j] == x[j+1]) && (y[j] == y[j+1])) continue; 716 for (k=0; k<=PointsPerInterval; ++k) { 717 t = (float) k * h[j] / (float) PointsPerInterval; 718 t2 = t * t; 719 t3 = t * t * t; 720 nx = x[j] + (int) (t * dx[j] + t2 * d2x[j]/2 + t3 * d3x[j]/6); 721 ny = y[j] + (int) (t * dy[j] + t2 * d2y[j]/2 + t3 * d3y[j]/6); 722 HGtline(lx, ly, nx, ny); 723 lx = nx; 724 ly = ny; 725 } /* end for k */ 726 } /* end for j */ 727 } /* end HGCurve */ 728 729 730 /*---------------------------------------------------------------------------- 731 * Routine: HGtline(xstart, ystart, xend, yend) 732 * 733 * Results: Draws a line of proper thickness by calling "line" numerous 734 * times until the desired thickness is reached. 735 *----------------------------------------------------------------------------*/ 736 737 HGtline(x0, y0, x1, y1) 738 register int x0; 739 register int y0; 740 int x1; 741 int y1; 742 { 743 register int xs; 744 register int xe; 745 register int ys; 746 register int ye; 747 double theta, wx, wy, xx, xy; 748 int addln, j, xdir, ydir, dx, dy; 749 750 751 xdir = ydir = 1; 752 dx = x1 - x0; /* calculate direction to move to */ 753 dy = y1 - y0; /* move to draw additional lines if needed */ 754 if (dx < 0 ) { /* for extra thickness */ 755 dx = -dx; 756 xdir = -1; 757 } 758 if (dy < 0 ) { 759 dy = -dy; 760 ydir = -1; 761 } 762 763 addln = linethickness / 2; 764 RoundEnd (x0, y0, linethickness, TRUE); /* add rounded end */ 765 766 for (j=(-addln); j<=addln; ++j) { 767 if (dy == 0) { 768 xs = x0; 769 xe = x1; 770 ys = ye = y0 + j; 771 } 772 if (dx == 0) { 773 ys = y0; 774 ye = y1; 775 xs = xe = x0 + j; 776 } 777 if ((dx != 0) && (dy != 0)) { 778 theta = pi / 2.0 - atan( ((double) dx)/((double) dy) ); 779 wx = j * sin(theta); 780 wy = j * cos(theta); 781 xs = x0 + (int) (wx * xdir + 0.4); 782 ys = y0 - (int) (wy * ydir + 0.4); 783 xe = x1 + (int) (wx * xdir + 0.4); 784 ye = y1 - (int) (wy * ydir + 0.4); 785 } 786 line(xs, ys, xe, ye); 787 } /* end for */ 788 789 RoundEnd(x1, y1, linethickness, TRUE); /* add rounded end */ 790 } /* end HGtline */ 791