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