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