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