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