src/preproc/grn/hgraph.cpp

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