src/preproc/grn/hgraph.cc

   1 /* Last non-groff version: hgraph.c  1.14 (Berkeley) 84/11/27
   2  *
   3  * This file contains the graphics routines for converting gremlin pictures
   4  * to troff input.
   5  */
   6
   7 #include "gprint.h"
   8
   9
  10 #define MAXVECT 40
  11 #define MAXPOINTS       200
  12 #define LINELENGTH      1
  13 #define PointsPerInterval 64
  14 #define pi              3.14159265358979324
  15 #define twopi           (2.0 * pi)
  16 #define len(a, b)       hypot((double)(b.x-a.x), (double)(b.y-a.y))
  17
  18
  19 extern int dotshifter;          /* for the length of dotted curves */
  20
  21 extern int style[];             /* line and character styles */
  22 extern double thick[];
  23 extern char *tfont[];
  24 extern int tsize[];
  25 extern int stipple_index[];     /* stipple font index for stipples 0 - 16 */
  26 extern char *stipple;           /* stipple type (cf or ug) */
  27
  28
  29 extern double troffscale;       /* imports from main.c */
  30 extern double linethickness;
  31 extern int linmod;
  32 extern int lastx;
  33 extern int lasty;
  34 extern int lastyline;
  35 extern int ytop;
  36 extern int ybottom;
  37 extern int xleft;
  38 extern int xright;
  39 extern enum {
  40   OUTLINE, FILL, BOTH
  41 } polyfill;
  42
  43 extern double adj1;
  44 extern double adj2;
  45 extern double adj3;
  46 extern double adj4;
  47 extern int res;
  48
  49 void HGSetFont(int font, int size);
  50 void HGPutText(int justify, POINT pnt, register char *string);
  51 void HGSetBrush(int mode);
  52 void tmove2(int px, int py);
  53 void doarc(POINT cp, POINT sp, int angle);
  54 void tmove(POINT * ptr);
  55 void cr(void);
  56 void drawwig(POINT * ptr);
  57 void HGtline(int x1, int y1);
  58 void dx(double x);
  59 void dy(double y);
  60 void HGArc(register int cx, register int cy, int px, int py, int angle);
  61 void picurve(register int *x, register int *y, int npts);
  62 void Paramaterize(int x[], int y[], float h[], int n);
  63 void PeriodicSpline(float h[], int z[],
  64                     float dz[], float d2z[], float d3z[],
  65                     int npoints);
  66 void NaturalEndSpline(float h[], int z[],
  67                       float dz[], float d2z[], float d3z[],
  68                       int npoints);
  69
  70
  71
  72 /*----------------------------------------------------------------------------*
  73  | Routine:     HGPrintElt (element_pointer, baseline)
  74  |
  75  | Results:     Examines a picture element and calls the appropriate
  76  |              routine(s) to print them according to their type.  After the
  77  |              picture is drawn, current position is (lastx, lasty).
  78  *----------------------------------------------------------------------------*/
  79
  80 void
  81 HGPrintElt(ELT *element,
  82            int baseline)
  83 {
  84   register POINT *p1;
  85   register POINT *p2;
  86   register int length;
  87   register int graylevel;
  88
  89   if (!DBNullelt(element) && !Nullpoint((p1 = element->ptlist))) {
  90     /* p1 always has first point */
  91     if (TEXT(element->type)) {
  92       HGSetFont(element->brushf, element->size);
  93       switch (element->size) {
  94       case 1:
  95         p1->y += adj1;
  96         break;
  97       case 2:
  98         p1->y += adj2;
  99         break;
 100       case 3:
 101         p1->y += adj3;
 102         break;
 103       case 4:
 104         p1->y += adj4;
 105         break;
 106       default:
 107         break;
 108       }
 109       HGPutText(element->type, *p1, element->textpt);
 110     } else {
 111       if (element->brushf)              /* if there is a brush, the */
 112         HGSetBrush(element->brushf);    /* graphics need it set     */
 113
 114       switch (element->type) {
 115
 116       case ARC:
 117         p2 = PTNextPoint(p1);
 118         tmove(p2);
 119         doarc(*p1, *p2, element->size);
 120         cr();
 121         break;
 122
 123       case CURVE:
 124         length = 0;             /* keep track of line length */
 125         drawwig(p1);
 126         cr();
 127         break;
 128
 129       case VECTOR:
 130         length = 0;             /* keep track of line length so */
 131         tmove(p1);              /* single lines don't get long  */
 132         while (!Nullpoint((p1 = PTNextPoint(p1)))) {
 133           HGtline((int) (p1->x * troffscale),
 134                   (int) (p1->y * troffscale));
 135           if (length++ > LINELENGTH) {
 136             length = 0;
 137             printf("\\\n");
 138           }
 139         }                       /* end while */
 140         cr();
 141         break;
 142
 143       case POLYGON:
 144         {
 145           /* brushf = style of outline; size = color of fill:
 146            * on first pass (polyfill=FILL), do the interior using 'P'
 147            *    unless size=0
 148            * on second pass (polyfill=OUTLINE), do the outline using a series
 149            *    of vectors. It might make more sense to use \D'p ...',
 150            *    but there is no uniform way to specify a 'fill character'
 151            *    that prints as 'no fill' on all output devices (and
 152            *    stipple fonts).
 153            * If polyfill=BOTH, just use the \D'p ...' command.
 154            */
 155           float firstx = p1->x;
 156           float firsty = p1->y;
 157
 158           length = 0;           /* keep track of line length so */
 159                                 /* single lines don't get long  */
 160
 161           if (polyfill == FILL || polyfill == BOTH) {
 162             /* do the interior */
 163             char command = (polyfill == BOTH && element->brushf) ? 'p' : 'P';
 164
 165             /* include outline, if there is one and */
 166             /* the -p flag was set                  */
 167
 168             /* switch based on what gremlin gives */
 169             switch (element->size) {
 170             case 1:
 171               graylevel = 1;
 172               break;
 173             case 3:
 174               graylevel = 2;
 175               break;
 176             case 12:
 177               graylevel = 3;
 178               break;
 179             case 14:
 180               graylevel = 4;
 181               break;
 182             case 16:
 183               graylevel = 5;
 184               break;
 185             case 19:
 186               graylevel = 6;
 187               break;
 188             case 21:
 189               graylevel = 7;
 190               break;
 191             case 23:
 192               graylevel = 8;
 193               break;
 194             default:            /* who's giving something else? */
 195               graylevel = NSTIPPLES;
 196               break;
 197             }
 198             /* int graylevel = element->size; */
 199
 200             if (graylevel < 0)
 201               break;
 202             if (graylevel > NSTIPPLES)
 203               graylevel = NSTIPPLES;
 204             printf("\\h'-%du'\\D'f %du'",
 205                    stipple_index[graylevel],
 206                    stipple_index[graylevel]);
 207             cr();
 208             tmove(p1);
 209             printf("\\D'%c", command);
 210
 211             while (!Nullpoint((PTNextPoint(p1)))) {
 212               p1 = PTNextPoint(p1);
 213               dx((double) p1->x);
 214               dy((double) p1->y);
 215               if (length++ > LINELENGTH) {
 216                 length = 0;
 217                 printf("\\\n");
 218               }
 219             } /* end while */
 220
 221             /* close polygon if not done so by user */
 222             if ((firstx != p1->x) || (firsty != p1->y)) {
 223               dx((double) firstx);
 224               dy((double) firsty);
 225             }
 226             putchar('\'');
 227             cr();
 228             break;
 229           }
 230           /* else polyfill == OUTLINE; only draw the outline */
 231           if (!(element->brushf))
 232             break;
 233           length = 0;           /* keep track of line length */
 234           tmove(p1);
 235
 236           while (!Nullpoint((PTNextPoint(p1)))) {
 237             p1 = PTNextPoint(p1);
 238             HGtline((int) (p1->x * troffscale),
 239                     (int) (p1->y * troffscale));
 240             if (length++ > LINELENGTH) {
 241               length = 0;
 242               printf("\\\n");
 243             }
 244           }                     /* end while */
 245
 246           /* close polygon if not done so by user */
 247           if ((firstx != p1->x) || (firsty != p1->y)) {
 248             HGtline((int) (firstx * troffscale),
 249                     (int) (firsty * troffscale));
 250           }
 251           cr();
 252           break;
 253         }                       /* end case POLYGON */
 254       }                         /* end switch */
 255     }                           /* end else Text */
 256   }                             /* end if */
 257 }                               /* end PrintElt */
 258
 259
 260 /*----------------------------------------------------------------------------*
 261  | Routine:     HGPutText (justification, position_point, string)
 262  |
 263  | Results:     Given the justification, a point to position with, and a
 264  |              string to put, HGPutText first sends the string into a
 265  |              diversion, moves to the positioning point, then outputs
 266  |              local vertical and horizontal motions as needed to justify
 267  |              the text.  After all motions are done, the diversion is
 268  |              printed out.
 269  *----------------------------------------------------------------------------*/
 270
 271 void
 272 HGPutText(int justify,
 273           POINT pnt,
 274           register char *string)
 275 {
 276   int savelasty = lasty;        /* vertical motion for text is to be */
 277                                 /* ignored.  Save current y here     */
 278
 279   printf(".nr g8 \\n(.d\n");    /* save current vertical position. */
 280   printf(".ds g9 \"");          /* define string containing the text. */
 281   while (*string) {             /* put out the string */
 282     if (*string == '\\' &&
 283         *(string + 1) == '\\') {        /* one character at a */
 284       printf("\\\\\\");                 /* time replacing //  */
 285       string++;                         /* by //// to prevent */
 286     }                                   /* interpretation at  */
 287     printf("%c", *(string++));          /* printout time      */
 288   }
 289   printf("\n");
 290
 291   tmove(&pnt);                  /* move to positioning point */
 292
 293   switch (justify) {
 294     /* local vertical motions                                            */
 295     /* (the numbers here are used to be somewhat compatible with gprint) */
 296   case CENTLEFT:
 297   case CENTCENT:
 298   case CENTRIGHT:
 299     printf("\\v'0.85n'");       /* down half */
 300     break;
 301
 302   case TOPLEFT:
 303   case TOPCENT:
 304   case TOPRIGHT:
 305     printf("\\v'1.7n'");        /* down whole */
 306   }
 307
 308   switch (justify) {
 309     /* local horizontal motions */
 310   case BOTCENT:
 311   case CENTCENT:
 312   case TOPCENT:
 313     printf("\\h'-\\w'\\*(g9'u/2u'");    /* back half */
 314     break;
 315
 316   case BOTRIGHT:
 317   case CENTRIGHT:
 318   case TOPRIGHT:
 319     printf("\\h'-\\w'\\*(g9'u'");       /* back whole */
 320   }
 321
 322   printf("\\&\\*(g9\n");        /* now print the text. */
 323   printf(".sp |\\n(g8u\n");     /* restore vertical position */
 324   lasty = savelasty;            /* vertical position restored to where it */
 325   lastx = xleft;                /* was before text, also horizontal is at */
 326                                 /* left                                   */
 327 }                               /* end HGPutText */
 328
 329
 330 /*----------------------------------------------------------------------------*
 331  | Routine:     doarc (center_point, start_point, angle)
 332  |
 333  | Results:     Produces either drawarc command or a drawcircle command
 334  |              depending on the angle needed to draw through.
 335  *----------------------------------------------------------------------------*/
 336
 337 void
 338 doarc(POINT cp,
 339       POINT sp,
 340       int angle)
 341 {
 342   if (angle)                    /* arc with angle */
 343     HGArc((int) (cp.x * troffscale), (int) (cp.y * troffscale),
 344           (int) (sp.x * troffscale), (int) (sp.y * troffscale), angle);
 345   else                          /* a full circle (angle == 0) */
 346     HGArc((int) (cp.x * troffscale), (int) (cp.y * troffscale),
 347           (int) (sp.x * troffscale), (int) (sp.y * troffscale), 0);
 348 }
 349
 350
 351 /*----------------------------------------------------------------------------*
 352  | Routine:     HGSetFont (font_number, Point_size)
 353  |
 354  | Results:     ALWAYS outputs a .ft and .ps directive to troff.  This is
 355  |              done because someone may change stuff inside a text string.
 356  |              Changes thickness back to default thickness.  Default
 357  |              thickness depends on font and pointsize.
 358  *----------------------------------------------------------------------------*/
 359
 360 void
 361 HGSetFont(int font,
 362           int size)
 363 {
 364   printf(".ft %s\n"
 365          ".ps %d\n", tfont[font - 1], tsize[size - 1]);
 366   linethickness = DEFTHICK;
 367 }
 368
 369
 370 /*----------------------------------------------------------------------------*
 371  | Routine:     HGSetBrush (line_mode)
 372  |
 373  | Results:     Generates the troff commands to set up the line width and
 374  |              style of subsequent lines.  Does nothing if no change is
 375  |              needed.
 376  |
 377  | Side Efct:   Sets `linmode' and `linethicknes'.
 378  *----------------------------------------------------------------------------*/
 379
 380 void
 381 HGSetBrush(int mode)
 382 {
 383   register int printed = 0;
 384
 385   if (linmod != style[--mode]) {
 386     /* Groff doesn't understand \Ds, so we take it out */
 387     /* printf ("\\D's %du'", linmod = style[mode]); */
 388     linmod = style[mode];
 389     printed = 1;
 390   }
 391   if (linethickness != thick[mode]) {
 392     linethickness = thick[mode];
 393     printf("\\h'-%.2lfp'\\D't %.2lfp'", linethickness, linethickness);
 394     printed = 1;
 395   }
 396   if (printed)
 397     cr();
 398 }
 399
 400
 401 /*----------------------------------------------------------------------------*
 402  | Routine:     dx (x_destination)
 403  |
 404  | Results:     Scales and outputs a number for delta x (with a leading
 405  |              space) given `lastx' and x_destination.
 406  |
 407  | Side Efct:   Resets `lastx' to x_destination.
 408  *----------------------------------------------------------------------------*/
 409
 410 void
 411 dx(double x)
 412 {
 413   register int ix = (int) (x * troffscale);
 414
 415   printf(" %du", ix - lastx);
 416   lastx = ix;
 417 }
 418
 419
 420 /*----------------------------------------------------------------------------*
 421  | Routine:     dy (y_destination)
 422  |
 423  | Results:     Scales and outputs a number for delta y (with a leading
 424  |              space) given `lastyline' and y_destination.
 425  |
 426  | Side Efct:   Resets `lastyline' to y_destination.  Since `line' vertical
 427  |              motions don't affect `page' ones, `lasty' isn't updated.
 428  *----------------------------------------------------------------------------*/
 429
 430 void
 431 dy(double y)
 432 {
 433   register int iy = (int) (y * troffscale);
 434
 435   printf(" %du", iy - lastyline);
 436   lastyline = iy;
 437 }
 438
 439
 440 /*----------------------------------------------------------------------------*
 441  | Routine:     tmove2 (px, py)
 442  |
 443  | Results:     Produces horizontal and vertical moves for troff given the
 444  |              pair of points to move to and knowing the current position.
 445  |              Also puts out a horizontal move to start the line.  This is
 446  |              a variation without the .sp command.
 447  *----------------------------------------------------------------------------*/
 448
 449 void
 450 tmove2(int px,
 451        int py)
 452 {
 453   register int dx;
 454   register int dy;
 455
 456   if (dy = py - lasty) {
 457     printf("\\v'%du'", dy);
 458   }
 459   lastyline = lasty = py;       /* lasty is always set to current */
 460   if (dx = px - lastx) {
 461     printf("\\h'%du'", dx);
 462     lastx = px;
 463   }
 464 }
 465
 466
 467 /*----------------------------------------------------------------------------*
 468  | Routine:     tmove (point_pointer)
 469  |
 470  | Results:     Produces horizontal and vertical moves for troff given the
 471  |              pointer of a point to move to and knowing the current
 472  |              position.  Also puts out a horizontal move to start the
 473  |              line.
 474  *----------------------------------------------------------------------------*/
 475
 476 void
 477 tmove(POINT * ptr)
 478 {
 479   register int ix = (int) (ptr->x * troffscale);
 480   register int iy = (int) (ptr->y * troffscale);
 481   register int dx;
 482   register int dy;
 483
 484   if (dy = iy - lasty) {
 485     printf(".sp %du\n", dy);
 486   }
 487   lastyline = lasty = iy;       /* lasty is always set to current */
 488   if (dx = ix - lastx) {
 489     printf("\\h'%du'", dx);
 490     lastx = ix;
 491   }
 492 }
 493
 494
 495 /*----------------------------------------------------------------------------*
 496  | Routine:     cr ( )
 497  |
 498  | Results:     Ends off an input line.  `.sp -1' is also added to counteract
 499  |              the vertical move done at the end of text lines.
 500  |
 501  | Side Efct:   Sets `lastx' to `xleft' for troff's return to left margin.
 502  *----------------------------------------------------------------------------*/
 503
 504 void
 505 cr(void)
 506 {
 507   printf("\n.sp -1\n");
 508   lastx = xleft;
 509 }
 510
 511
 512 /*----------------------------------------------------------------------------*
 513  | Routine:     line ( )
 514  |
 515  | Results:     Draws a single solid line to (x,y).
 516  *----------------------------------------------------------------------------*/
 517
 518 void
 519 line(int px,
 520      int py)
 521 {
 522   printf("\\D'l");
 523   printf(" %du", px - lastx);
 524   printf(" %du'", py - lastyline);
 525   lastx = px;
 526   lastyline = lasty = py;
 527 }
 528
 529
 530 /*----------------------------------------------------------------------------
 531  | Routine:     drawwig (ptr)
 532  |
 533  | Results:     The point sequence found in the structure pointed by ptr is
 534  |              placed in integer arrays for further manipulation by the
 535  |              existing routing.  With the proper parameters, HGCurve is
 536  |              called.
 537  *----------------------------------------------------------------------------*/
 538
 539 void
 540 drawwig(POINT * ptr)
 541 {
 542   register int npts;                    /* point list index */
 543   int x[MAXPOINTS], y[MAXPOINTS];       /* point list */
 544
 545   for (npts = 1; !Nullpoint(ptr); ptr = PTNextPoint(ptr), npts++) {
 546     x[npts] = (int) (ptr->x * troffscale);
 547     y[npts] = (int) (ptr->y * troffscale);
 548   }
 549   if (--npts) {
 550     /* HGCurve(&x[0], &y[0], npts); */  /*Gremlin looks different, so... */
 551     picurve(&x[0], &y[0], npts);
 552   }
 553 }
 554
 555
 556 /*----------------------------------------------------------------------------
 557  | Routine:     HGArc (xcenter, ycenter, xstart, ystart, angle)
 558  |
 559  | Results:     This routine plots an arc centered about (cx, cy) counter
 560  |              clockwise starting from the point (px, py) through `angle'
 561  |              degrees.  If angle is 0, a full circle is drawn.  It does so
 562  |              by creating a draw-path around the arc whose density of
 563  |              points depends on the size of the arc.
 564  *----------------------------------------------------------------------------*/
 565
 566 void
 567 HGArc(register int cx,
 568       register int cy,
 569       int px,
 570       int py,
 571       int angle)
 572 {
 573   double xs, ys, resolution, fullcircle;
 574   int m;
 575   register int mask;
 576   register int extent;
 577   register int nx;
 578   register int ny;
 579   register int length;
 580   register double epsilon;
 581
 582   xs = px - cx;
 583   ys = py - cy;
 584
 585   length = 0;
 586
 587   resolution = (1.0 + hypot(xs, ys) / res) * PointsPerInterval;
 588   /* mask = (1 << (int) log10(resolution + 1.0)) - 1; */
 589   (void) frexp(resolution, &m);         /* A bit more elegant than log10 */
 590   for (mask = 1; mask < m; mask = mask << 1);
 591   mask -= 1;
 592
 593   epsilon = 1.0 / resolution;
 594   fullcircle = (2.0 * pi) * resolution;
 595   if (angle == 0)
 596     extent = (int) fullcircle;
 597   else
 598     extent = (int) (angle * fullcircle / 360.0);
 599
 600   HGtline(px, py);
 601   while (--extent >= 0) {
 602     xs += epsilon * ys;
 603     nx = cx + (int) (xs + 0.5);
 604     ys -= epsilon * xs;
 605     ny = cy + (int) (ys + 0.5);
 606     if (!(extent & mask)) {
 607       HGtline(nx, ny);          /* put out a point on circle */
 608       if (length++ > LINELENGTH) {
 609         length = 0;
 610         printf("\\\n");
 611       }
 612     }
 613   }                             /* end for */
 614 }                               /* end HGArc */
 615
 616
 617 /*----------------------------------------------------------------------------
 618  | Routine:     picurve (xpoints, ypoints, num_of_points)
 619  |
 620  | Results:     Draws a curve delimited by (not through) the line segments
 621  |              traced by (xpoints, ypoints) point list.  This is the `Pic'
 622  |              style curve.
 623  *----------------------------------------------------------------------------*/
 624
 625 void
 626 picurve(register int *x,
 627         register int *y,
 628         int npts)
 629 {
 630   register int nseg;            /* effective resolution for each curve */
 631   register int xp;              /* current point (and temporary) */
 632   register int yp;
 633   int pxp, pyp;                 /* previous point (to make lines from) */
 634   int i;                        /* inner curve segment traverser */
 635   int length = 0;
 636   double w;                     /* position factor */
 637   double t1, t2, t3;            /* calculation temps */
 638
 639   if (x[1] == x[npts] && y[1] == y[npts]) {
 640     x[0] = x[npts - 1];         /* if the lines' ends meet, make */
 641     y[0] = y[npts - 1];         /* sure the curve meets          */
 642     x[npts + 1] = x[2];
 643     y[npts + 1] = y[2];
 644   } else {                      /* otherwise, make the ends of the  */
 645     x[0] = x[1];                /* curve touch the ending points of */
 646     y[0] = y[1];                /* the line segments                */
 647     x[npts + 1] = x[npts];
 648     y[npts + 1] = y[npts];
 649   }
 650
 651   pxp = (x[0] + x[1]) / 2;      /* make the last point pointers       */
 652   pyp = (y[0] + y[1]) / 2;      /* point to the start of the 1st line */
 653   tmove2(pxp, pyp);
 654
 655   for (; npts--; x++, y++) {    /* traverse the line segments */
 656     xp = x[0] - x[1];
 657     yp = y[0] - y[1];
 658     nseg = (int) hypot((double) xp, (double) yp);
 659     xp = x[1] - x[2];
 660     yp = y[1] - y[2];
 661                                 /* `nseg' is the number of line    */
 662                                 /* segments that will be drawn for */
 663                                 /* each curve segment.             */
 664     nseg = (int) ((double) (nseg + (int) hypot((double) xp, (double) yp)) /
 665                   res * PointsPerInterval);
 666
 667     for (i = 1; i < nseg; i++) {
 668       w = (double) i / (double) nseg;
 669       t1 = w * w;
 670       t3 = t1 + 1.0 - (w + w);
 671       t2 = 2.0 - (t3 + t1);
 672       xp = (((int) (t1 * x[2] + t2 * x[1] + t3 * x[0])) + 1) / 2;
 673       yp = (((int) (t1 * y[2] + t2 * y[1] + t3 * y[0])) + 1) / 2;
 674
 675       HGtline(xp, yp);
 676       if (length++ > LINELENGTH) {
 677         length = 0;
 678         printf("\\\n");
 679       }
 680     }
 681   }
 682 }
 683
 684
 685 /*----------------------------------------------------------------------------
 686  | Routine:     HGCurve(xpoints, ypoints, num_points)
 687  |
 688  | Results:     This routine generates a smooth curve through a set of
 689  |              points.  The method used is the parametric spline curve on
 690  |              unit knot mesh described in `Spline Curve Techniques' by
 691  |              Patrick Baudelaire, Robert Flegal, and Robert Sproull --
 692  |              Xerox Parc.
 693  *----------------------------------------------------------------------------*/
 694
 695 void
 696 HGCurve(int *x,
 697         int *y,
 698         int numpoints)
 699 {
 700   float h[MAXPOINTS], dx[MAXPOINTS], dy[MAXPOINTS];
 701   float d2x[MAXPOINTS], d2y[MAXPOINTS], d3x[MAXPOINTS], d3y[MAXPOINTS];
 702   float t, t2, t3;
 703   register int j;
 704   register int k;
 705   register int nx;
 706   register int ny;
 707   int lx, ly;
 708   int length = 0;
 709
 710   lx = x[1];
 711   ly = y[1];
 712   tmove2(lx, ly);
 713
 714   /*
 715    * Solve for derivatives of the curve at each point separately for x and y
 716    * (parametric).
 717    */
 718   Paramaterize(x, y, h, numpoints);
 719
 720   /* closed curve */
 721   if ((x[1] == x[numpoints]) && (y[1] == y[numpoints])) {
 722     PeriodicSpline(h, x, dx, d2x, d3x, numpoints);
 723     PeriodicSpline(h, y, dy, d2y, d3y, numpoints);
 724   } else {
 725     NaturalEndSpline(h, x, dx, d2x, d3x, numpoints);
 726     NaturalEndSpline(h, y, dy, d2y, d3y, numpoints);
 727   }
 728
 729   /*
 730    * generate the curve using the above information and PointsPerInterval
 731    * vectors between each specified knot.
 732    */
 733
 734   for (j = 1; j < numpoints; ++j) {
 735     if ((x[j] == x[j + 1]) && (y[j] == y[j + 1]))
 736       continue;
 737     for (k = 0; k <= PointsPerInterval; ++k) {
 738       t = (float) k *h[j] / (float) PointsPerInterval;
 739       t2 = t * t;
 740       t3 = t * t * t;
 741       nx = x[j] + (int) (t * dx[j] + t2 * d2x[j] / 2 + t3 * d3x[j] / 6);
 742       ny = y[j] + (int) (t * dy[j] + t2 * d2y[j] / 2 + t3 * d3y[j] / 6);
 743       HGtline(nx, ny);
 744       if (length++ > LINELENGTH) {
 745         length = 0;
 746         printf("\\\n");
 747       }
 748     }                           /* end for k */
 749   }                             /* end for j */
 750 }                               /* end HGCurve */
 751
 752
 753 /*----------------------------------------------------------------------------
 754  | Routine:     Paramaterize (xpoints, ypoints, hparams, num_points)
 755  |
 756  | Results:     This routine calculates parameteric values for use in
 757  |              calculating curves.  The parametric values are returned
 758  |              in the array h.  The values are an approximation of
 759  |              cumulative arc lengths of the curve (uses cord length).
 760  |              For additional information, see paper cited below.
 761  *----------------------------------------------------------------------------*/
 762
 763 void
 764 Paramaterize(int x[],
 765              int y[],
 766              float h[],
 767              int n)
 768 {
 769   register int dx;
 770   register int dy;
 771   register int i;
 772   register int j;
 773   float u[MAXPOINTS];
 774
 775   for (i = 1; i <= n; ++i) {
 776     u[i] = 0;
 777     for (j = 1; j < i; j++) {
 778       dx = x[j + 1] - x[j];
 779       dy = y[j + 1] - y[j];
 780       /* Here was overflowing, so I changed it.       */
 781       /* u[i] += sqrt ((double) (dx * dx + dy * dy)); */
 782       u[i] += hypot((double) dx, (double) dy);
 783     }
 784   }
 785   for (i = 1; i < n; ++i)
 786     h[i] = u[i + 1] - u[i];
 787 }                               /* end Paramaterize */
 788
 789
 790 /*----------------------------------------------------------------------------
 791  | Routine:     PeriodicSpline (h, z, dz, d2z, d3z, npoints)
 792  |
 793  | Results:     This routine solves for the cubic polynomial to fit a spline
 794  |              curve to the the points specified by the list of values.
 795  |              The Curve generated is periodic.  The algorithms for this
 796  |              curve are from the `Spline Curve Techniques' paper cited
 797  |              above.
 798  *----------------------------------------------------------------------------*/
 799
 800 void
 801 PeriodicSpline(float h[],       /* paramaterization  */
 802                int z[],         /* point list */
 803                float dz[],      /* to return the 1st derivative */
 804                float d2z[],     /* 2nd derivative */
 805                float d3z[],     /* 3rd derivative */
 806                int npoints)     /* number of valid points */
 807 {
 808   float d[MAXPOINTS];
 809   float deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS];
 810   float c[MAXPOINTS], r[MAXPOINTS], s[MAXPOINTS];
 811   int i;
 812
 813   /* step 1 */
 814   for (i = 1; i < npoints; ++i) {
 815     deltaz[i] = h[i] ? ((double) (z[i + 1] - z[i])) / h[i] : 0;
 816   }
 817   h[0] = h[npoints - 1];
 818   deltaz[0] = deltaz[npoints - 1];
 819
 820   /* step 2 */
 821   for (i = 1; i < npoints - 1; ++i) {
 822     d[i] = deltaz[i + 1] - deltaz[i];
 823   }
 824   d[0] = deltaz[1] - deltaz[0];
 825
 826   /* step 3a */
 827   a[1] = 2 * (h[0] + h[1]);
 828   b[1] = d[0];
 829   c[1] = h[0];
 830   for (i = 2; i < npoints - 1; ++i) {
 831     a[i] = 2 * (h[i - 1] + h[i]) -
 832            pow((double) h[i - 1], (double) 2.0) / a[i - 1];
 833     b[i] = d[i - 1] - h[i - 1] * b[i - 1] / a[i - 1];
 834     c[i] = -h[i - 1] * c[i - 1] / a[i - 1];
 835   }
 836
 837   /* step 3b */
 838   r[npoints - 1] = 1;
 839   s[npoints - 1] = 0;
 840   for (i = npoints - 2; i > 0; --i) {
 841     r[i] = -(h[i] * r[i + 1] + c[i]) / a[i];
 842     s[i] = (6 * b[i] - h[i] * s[i + 1]) / a[i];
 843   }
 844
 845   /* step 4 */
 846   d2z[npoints - 1] = (6 * d[npoints - 2] - h[0] * s[1]
 847                       - h[npoints - 1] * s[npoints - 2])
 848                      / (h[0] * r[1] + h[npoints - 1] * r[npoints - 2]
 849                       + 2 * (h[npoints - 2] + h[0]));
 850   for (i = 1; i < npoints - 1; ++i) {
 851     d2z[i] = r[i] * d2z[npoints - 1] + s[i];
 852   }
 853   d2z[npoints] = d2z[1];
 854
 855   /* step 5 */
 856   for (i = 1; i < npoints; ++i) {
 857     dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i + 1]) / 6;
 858     d3z[i] = h[i] ? (d2z[i + 1] - d2z[i]) / h[i] : 0;
 859   }
 860 }                               /* end PeriodicSpline */
 861
 862
 863 /*----------------------------------------------------------------------------
 864  | Routine:     NaturalEndSpline (h, z, dz, d2z, d3z, npoints)
 865  |
 866  | Results:     This routine solves for the cubic polynomial to fit a spline
 867  |              curve the the points specified by the list of values.  The
 868  |              alogrithms for this curve are from the `Spline Curve
 869  |              Techniques' paper cited above.
 870  *----------------------------------------------------------------------------*/
 871
 872 void
 873 NaturalEndSpline(float h[],     /* parameterization */
 874                  int z[],       /* Point list */
 875                  float dz[],    /* to return the 1st derivative */
 876                  float d2z[],   /* 2nd derivative */
 877                  float d3z[],   /* 3rd derivative */
 878                  int npoints)   /* number of valid points */
 879 {
 880   float d[MAXPOINTS];
 881   float deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS];
 882   int i;
 883
 884   /* step 1 */
 885   for (i = 1; i < npoints; ++i) {
 886     deltaz[i] = h[i] ? ((double) (z[i + 1] - z[i])) / h[i] : 0;
 887   }
 888   deltaz[0] = deltaz[npoints - 1];
 889
 890   /* step 2 */
 891   for (i = 1; i < npoints - 1; ++i) {
 892     d[i] = deltaz[i + 1] - deltaz[i];
 893   }
 894   d[0] = deltaz[1] - deltaz[0];
 895
 896   /* step 3 */
 897   a[0] = 2 * (h[2] + h[1]);
 898   b[0] = d[1];
 899   for (i = 1; i < npoints - 2; ++i) {
 900     a[i] = 2 * (h[i + 1] + h[i + 2]) -
 901             pow((double) h[i + 1], (double) 2.0) / a[i - 1];
 902     b[i] = d[i + 1] - h[i + 1] * b[i - 1] / a[i - 1];
 903   }
 904
 905   /* step 4 */
 906   d2z[npoints] = d2z[1] = 0;
 907   for (i = npoints - 1; i > 1; --i) {
 908     d2z[i] = (6 * b[i - 2] - h[i] * d2z[i + 1]) / a[i - 2];
 909   }
 910
 911   /* step 5 */
 912   for (i = 1; i < npoints; ++i) {
 913     dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i + 1]) / 6;
 914     d3z[i] = h[i] ? (d2z[i + 1] - d2z[i]) / h[i] : 0;
 915   }
 916 }                               /* end NaturalEndSpline */
 917
 918
 919 /*----------------------------------------------------------------------------*
 920  | Routine:     change (x_position, y_position, visible_flag)
 921  |
 922  | Results:     As HGtline passes from the invisible to visible (or vice
 923  |              versa) portion of a line, change is called to either draw
 924  |              the line, or initialize the beginning of the next one.
 925  |              Change calls line to draw segments if visible_flag is set
 926  |              (which means we're leaving a visible area).
 927  *----------------------------------------------------------------------------*/
 928
 929 void
 930 change(register int x,
 931        register int y,
 932        register int vis)
 933 {
 934   static int xorg;
 935   static int yorg;
 936   static int length = 0;
 937
 938   if (vis) {                    /* leaving a visible area, draw it. */
 939     line(x, y);
 940     if (length++ > LINELENGTH) {
 941       length = 0;
 942       printf("\\\n");
 943     }
 944   } else {                      /* otherwise, we're entering one, remember */
 945                                 /* beginning                               */
 946     tmove2(x, y);
 947   }
 948 }
 949
 950
 951 /*----------------------------------------------------------------------------
 952  | Routine:     HGtline (xstart, ystart, xend, yend)
 953  |
 954  | Results:     Draws a line from current position to (x1,y1) using line(x1,
 955  |              y1) to place individual segments of dotted or dashed lines.
 956  *----------------------------------------------------------------------------*/
 957
 958 void
 959 HGtline(int x1,
 960         int y1)
 961 {
 962   register int x0 = lastx;
 963   register int y0 = lasty;
 964   register int dx;
 965   register int dy;
 966   register int oldcoord;
 967   register int res1;
 968   register int visible;
 969   register int res2;
 970   register int xinc;
 971   register int yinc;
 972   register int dotcounter;
 973
 974   if (linmod == SOLID) {
 975     line(x1, y1);
 976     return;
 977   }
 978
 979   /* for handling different resolutions */
 980   dotcounter = linmod << dotshifter;
 981
 982   xinc = 1;
 983   yinc = 1;
 984   if ((dx = x1 - x0) < 0) {
 985     xinc = -xinc;
 986     dx = -dx;
 987   }
 988   if ((dy = y1 - y0) < 0) {
 989     yinc = -yinc;
 990     dy = -dy;
 991   }
 992   res1 = 0;
 993   res2 = 0;
 994   visible = 0;
 995   if (dx >= dy) {
 996     oldcoord = y0;
 997     while (x0 != x1) {
 998       if ((x0 & dotcounter) && !visible) {
 999         change(x0, y0, 0);
1000         visible = 1;
1001       } else if (visible && !(x0 & dotcounter)) {
1002         change(x0 - xinc, oldcoord, 1);
1003         visible = 0;
1004       }
1005       if (res1 > res2) {
1006         oldcoord = y0;
1007         res2 += dx - res1;
1008         res1 = 0;
1009         y0 += yinc;
1010       }
1011       res1 += dy;
1012       x0 += xinc;
1013     }
1014   } else {
1015     oldcoord = x0;
1016     while (y0 != y1) {
1017       if ((y0 & dotcounter) && !visible) {
1018         change(x0, y0, 0);
1019         visible = 1;
1020       } else if (visible && !(y0 & dotcounter)) {
1021         change(oldcoord, y0 - yinc, 1);
1022         visible = 0;
1023       }
1024       if (res1 > res2) {
1025         oldcoord = x0;
1026         res2 += dy - res1;
1027         res1 = 0;
1028         x0 += xinc;
1029       }
1030       res1 += dx;
1031       y0 += yinc;
1032     }
1033   }
1034   if (visible)
1035     change(x1, y1, 1);
1036   else
1037     change(x1, y1, 0);
1038 }
1039
1040 /* EOF */