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