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