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