1 /* graph.c	1.7	83/07/09
2  *
3  *	This file contains the functions for producing the graphics
4  *   images in the varian/versatec drivers for ditroff.
5  */
6 
7 
8 #include <stdio.h>
9 #include <ctype.h>
10 #include <math.h>
11 
12 
13 #define TRUE	1
14 #define FALSE	0
15 				/* imports from dver.c */
16 #define  hmot(n)	hpos += n;
17 #define  vmot(n)	vgoto(vpos + n);
18 
19 extern int hpos;
20 extern int vpos;
21 extern vgoto();
22 extern point();
23 
24 #define MAXPOINTS 200	/* number of points legal for a curve */
25 
26 #define SOLID -1	/* line styles:  these are used as bit masks to */
27 #define DOTTED 004	/* create the right style lines. */
28 #define DASHED 020
29 #define DOTDASHED 024
30 #define LONGDASHED 074
31 				/* constants... */
32 #define  pi		3.14159265358979324
33 #define  log2_10	3.3219280948873623
34 
35 
36 int	linethickness = 3;	/* number of pixels wide to make lines */
37 int	linmod = SOLID;		/* type of line (SOLID, DOTTED, DASHED...) */
38 
39 
40 
41 /*----------------------------------------------------------------------------
42  * Routine:	drawline (horizontal_motion, vertical_motion)
43  *
44  * Results:	Draws a line of "linethickness" width and "linmod" style
45  *		from current (hpos, vpos) to (hpos + dh, vpos + dv).
46  *
47  * Side Efct:	Resulting position is at end of line (hpos + dh, vpos + dv)
48  *----------------------------------------------------------------------------*/
49 
50 drawline(dh, dv)
51 register int dh;
52 register int dv;
53 {
54     HGtline (hpos, vpos, hpos + dh, vpos + dv);
55     hmot (dh);					/* new position is at */
56     vmot (dv);					/* the end of the line */
57 }
58 
59 
60 /*----------------------------------------------------------------------------
61  * Routine:	drawcirc (diameter)
62  *
63  * Results:	Draws a circle with leftmost point at current (hpos, vpos)
64  *		with the given diameter d.
65  *
66  * Side Efct:	Resulting position is at (hpos + diameter, vpos)
67  *----------------------------------------------------------------------------*/
68 
69 drawcirc(d)
70 register int d;
71 {			/* 0.0 is the angle to sweep the arc: = full circle */
72     HGArc (hpos + d/2, vpos, hpos, vpos, 0.0);
73     hmot (d);			/* new postion is the right of the circle */
74 }
75 
76 
77 /*----------------------------------------------------------------------------
78  * Routine:	drawellip (horizontal_diameter, vertical_diameter)
79  *
80  * Results:	Draws regular ellipses given the major "diameters."  It does
81  *		so by drawing many small lines, every other pixel.  The ellipse
82  *		formula:  ((x-x0)/hrad)**2 + ((y-y0)/vrad)**2 = 1 is used,
83  *		converting to:  y = y0 +- vrad * sqrt(1 - ((x-x0)/hrad)**2).
84  *		The line segments are duplicated (mirrored) on either side of
85  *		the horizontal "diameter".
86  *
87  * Side Efct:	Resulting position is at (hpos + hd, vpos).
88  *
89  * Bugs:	Odd numbered horizontal axes are rounded up to even numbers.
90  *----------------------------------------------------------------------------*/
91 
92 drawellip(hd, vd)
93 register int hd;
94 int vd;
95 {
96     register int bx;		/* multiplicative x factor */
97     register int x;		/* x position drawing to */
98     register int yk;		/* the square-root term */
99     register int y;		/* y position drawing to */
100     double k1;			/* k? are constants depending on parameters */
101     int k2, oldy1, oldy2;	/* oldy? are last y points drawn to */
102 
103 
104     hd = 2 * ((hd + 1) / 2);	/* don't accept odd diameters */
105 
106     bx = 4 * (hpos + hd);
107     x = hpos;
108     k1 = vd / (2.0 * hd);
109     k2 = hd * hd - 4 * (hpos + hd/2) * (hpos + hd/2);
110     oldy1 = vpos;
111     oldy2 = vpos;
112 
113     hmot (hd);		/* end position is the right-hand side of the ellipse */
114 
115     do {
116 	yk = (int) (k1 * sqrt((double) (k2 + (bx -= 8) * (x += 2))));
117 
118 	HGtline (x-2, oldy1, x, y = vpos + yk);    /* top half of ellipse */
119 	oldy1 = y;
120 	HGtline (x-2, oldy2, x, y = vpos - yk);	  /* bottom half of ellipse */
121 	oldy2 = y;
122 
123     } while (hd -= 2);
124 }
125 
126 
127 /*----------------------------------------------------------------------------
128  * Routine:	drawarc (xcenter, ycenter, xpoint, ypoint)
129  *
130  * Results:	Draws an arc starting at current (hpos, vpos).  Center is
131  *		at (hpos + cdh, vpos + cdv) and the terminating point is
132  *		at <center> + (pdh, pdv).  The angle between the lines
133  *		formed by the starting, ending, and center points is figured
134  *		first, then the points and angle are sent to HGArc for the
135  *		drawing.
136  *
137  * Side Efct:	Resulting position is at the last point of the arc.
138  *----------------------------------------------------------------------------*/
139 
140 drawarc (cdh, cdv, pdh, pdv)
141 register int cdh;
142 register int cdv;
143 register int pdh;
144 register int pdv;
145 {
146     register double angle;
147 				/* figure angle from the three points...*/
148 				/* and convert (and round) to degrees */
149     angle = (atan2((double) pdh, (double) pdv)
150 		- atan2 ((double) -cdh, (double) -cdv)) * 180.0 / pi;
151 				/* "normalize" and round */
152     angle += (angle < 0.0)  ?  360.5 : 0.5;
153 
154     HGArc(hpos + cdh, vpos + cdv, hpos, vpos, (int) angle);
155     hmot(cdh + pdh);
156     vmot(cdv + pdv);
157 }
158 
159 
160 /*----------------------------------------------------------------------------
161  * Routine:	drawwig (character_buffer, file_pointer, type_flag)
162  *
163  * Results:	Given the starting position, the motion list in buf, and any
164  *		extra characters from fp (terminated by a \n), drawwig sets
165  *		up a point list to make a spline from.  If "pic" is set picurve
166  *		is called to draw the curve in pic style; else it calls HGCurve
167  *		for the gremlin-style curve.
168  *
169  * Side Efct:	Resulting position is reached from adding successive motions
170  *		to the current position.
171  *----------------------------------------------------------------------------*/
172 
173 drawwig (buf, fp, pic)
174 char *buf;
175 FILE *fp;
176 int pic;
177 {
178     register int len = strlen(buf);	/* length of the string in "buf" */
179     register int npts = 2;		/* point list index */
180     register char *ptr = buf;		/* "walking" pointer into buf */
181     int x[MAXPOINTS], y[MAXPOINTS];	/* point list */
182 
183     while (*ptr == ' ') ptr++;		/* skip any leading spaces */
184     x[1] = hpos;		/* the curve starts at the */
185     y[1] = vpos;		/* current position */
186 
187     while (*ptr != '\n') {		/* curve commands end with a '\n' */
188 	hmot(atoi(ptr));		/* convert motion to curve points */
189 	x[npts] = hpos;			/* and remember them */
190 	while (isdigit(*++ptr));		/* skip number*/
191 	while (*++ptr == ' ');		/* skip spaces 'tween numbers */
192 	vmot(atoi(ptr));
193 	y[npts] = vpos;
194 	while (isdigit(*++ptr));
195 	while (*ptr == ' ') ptr++;
196 				/* if the amount we read wasn't the */
197 		 		/*    whole thing, read some more in */
198 	if (len - (ptr - buf) < 15 && *(buf + len - 1) != '\n') {
199 	    char *cop = buf;
200 
201 	    while (*cop++ = *ptr++);	/* copy what's left to the beginning */
202 	    fgets ((cop - 1), len - (cop - buf), fp);
203 	    ptr = buf;
204 	}
205 	if (npts < MAXPOINTS - 1)	/* if too many points, forget some */
206 	    npts++;
207     }
208     npts--;	/* npts must point to the last coordinate in x and y */
209 				/* now, actually DO the curve */
210     if (pic)
211 	picurve(x, y, npts);
212     else
213 	HGCurve(x, y, npts);
214 }
215 
216 
217 /*----------------------------------------------------------------------------*
218  | Routine:	drawthick (thickness)
219  |
220  | Results:	sets the variable "linethickness" to the given size.
221  |		NO motion is involved.
222  *----------------------------------------------------------------------------*/
223 
224 drawthick(s)
225 int s;
226 {
227     linethickness = s;
228 }
229 
230 
231 /*----------------------------------------------------------------------------*
232  | Routine:	drawstyle (style_bit_map)
233  |
234  | Results:	sets the variable "linmod" to the given bit map.
235  |		NO motion is involved.
236  *----------------------------------------------------------------------------*/
237 
238 drawstyle(s)
239 int s;
240 {
241     linmod = s;
242 }
243 
244 
245 /*----------------------------------------------------------------------------
246  * Routine:	picurve (xpoints, ypoints, num_of_points)
247  *
248  * Results:	Draws a curve delimited by (not through) the line segments
249  *		traced by (xpoints, ypoints) point list.  This is the "Pic"
250  *		style curve.
251  *
252  * Bugs:	does nothing yet....
253  *----------------------------------------------------------------------------*/
254 
255 picurve (x, y, npts)
256 int x[MAXPOINTS];
257 int y[MAXPOINTS];
258 int npts;
259 {
260     register int i;		/* line segment traverser */
261     register int nseg;		/* effective resolution for each curve */
262     register float w;		/* position factor */
263     register int xp;		/* current point (and intermediary) */
264     register int yp;
265     int pxp, pyp;		/* previous point (to make lines from) */
266     float t1, t2, t3;		/* calculation temps */
267     int j;			/* inner curve segment traverser */
268 
269 
270     if (x[1] == x[npts] && y[1] == y[npts]) {
271 	x[0] = x[npts - 1];		/* if the lines' ends meet, make */
272 	y[0] = y[npts - 1];		/* sure the curve meets */
273 	x[npts + 1] = x[2];
274 	y[npts + 1] = y[2];
275     } else {				/* otherwise, make the ends of the */
276 	x[0] = x[1];			/* curve touch the ending points of */
277 	y[0] = y[1];			/* the line segments */
278 	x[npts + 1] = x[npts];
279 	y[npts + 1] = y[npts];
280     }
281 
282     pxp = (x[0] + x[1]) / 2;		/* make the last point pointers */
283     pyp = (y[0] + y[1]) / 2;		/* point to the start of the 1st line */
284 
285     for (i = 0; i < npts; i++) {	/* traverse the line segments */
286 	xp = x[i] - x[i+1];
287 	yp = y[i] - y[i+1];
288 	nseg = (int) sqrt((double)(xp * xp + yp * yp));
289 	xp = x[i+1] - x[i+2];
290 	yp = y[i+1] - y[i+2];		/* "nseg" is the number of line */
291 					/* segments that will be drawn for */
292 					/* each curve segment.  ">> 3" is */
293 					/* dropping the resolution ( == / 8) */
294 	nseg = (nseg + (int) sqrt((double)(xp * xp + yp * yp))) >> 3;
295 
296 	for (j = 1; j < nseg; j++) {
297 	    w = (float) j / (float) nseg;
298 	    t1 = 0.5 * w * w;
299 	    w -= 0.5;
300 	    t2 = 0.75 - w * w ;
301 	    w -= 0.5;
302 	    t3 = 0.5 * w * w;
303 	    xp = t1 * x[i+2] + t2 * x[i+1] + t3 * x[i] + 0.5;
304 	    yp = t1 * y[i+2] + t2 * y[i+1] + t3 * y[i] + 0.5;
305 
306 	    HGtline(pxp, pyp, xp, yp);
307 	    pxp = xp;
308 	    pyp = yp;
309 	}
310     }
311 }
312 
313 
314 /*----------------------------------------------------------------------------
315  * Routine:	line (xstart, ystart, xend, yend)
316  *
317  * Results:	Draws a one-pixel wide line from (x0, y0) to
318  *		(x1, y1) using point(x,y) to place the dots.
319  *----------------------------------------------------------------------------*/
320 
321 line(x0, y0, x1, y1)
322 register int x0;
323 register int y0;
324 int x1;
325 int y1;
326 {
327     register int dx;
328     register int dy;
329     register int res1;
330     int res2;
331     int xinc;
332     int yinc;
333     int slope;
334 
335     xinc = 1;
336     yinc = 1;
337     if ((dx = x1-x0) < 0) {
338         xinc = -1;
339         dx = -dx;
340     }
341     if ((dy = y1-y0) < 0) {
342         yinc = -1;
343         dy = -dy;
344     }
345     slope = xinc*yinc;
346     res1 = 0;
347     res2 = 0;
348     if (dx >= dy)
349         while (x0 != x1) {
350             if((x0+slope*y0)&linmod) point(x0, y0);
351             if (res1 > res2) {
352                 res2 += dx - res1;
353                 res1 = 0;
354                 y0 += yinc;
355             }
356             res1 += dy;
357             x0 += xinc;
358         }
359     else
360         while (y0 != y1) {
361             if((x0+slope*y0)&linmod) point(x0, y0);
362             if (res1 > res2) {
363                 res2 += dy - res1;
364                 res1 = 0;
365                 x0 += xinc;
366             }
367             res1 += dx;
368             y0 += yinc;
369         }
370     if((x1+slope*y1)&linmod) point(x1, y1);
371 }
372 
373 
374 /*----------------------------------------------------------------------------
375  * Routine:	HGArc (xcenter, ycenter, xstart, ystart, angle)
376  *
377  * Results:	This routine plots an arc centered about (cx, cy) counter
378  *		clockwise starting from the point (px, py) through 'angle'
379  *		degrees.  If angle is 0, a full circle is drawn.
380  *		It does so by calling RoundEnd (fat point maker) for points
381  *		along the circle with density depending on the circle's size.
382  *----------------------------------------------------------------------------*/
383 
384 HGArc(cx,cy,px,py,angle)
385 register int cx;
386 register int cy;
387 register int px;
388 register int py;
389 register int angle;
390 {
391     double xs, ys, resolution, epsilon, degreesperpoint, fullcircle;
392     double t1, t2;
393     int i, extent, nx, ny;
394 
395     xs = px - cx;
396     ys = py - cy;
397 
398 /* calculate drawing parameters */
399 
400     t1 = log10(sqrt( xs * xs + ys * ys)) * log2_10;
401     t1 = ceil(t1);
402     resolution = pow(2.0, t1);
403     epsilon = 1.0 / resolution;
404     fullcircle = 2 * pi * resolution;
405     fullcircle = ceil(fullcircle);
406     degreesperpoint = 360.0 / fullcircle;
407 
408     if (angle == 0) extent = fullcircle;
409     else extent = angle/degreesperpoint;
410 
411     for (i=0; i<extent; ++i) {
412         xs += epsilon * ys;
413         nx = (int) (xs + cx + 0.5);
414         ys -= epsilon * xs;
415         ny = (int) (ys + cy + 0.5);
416         RoundEnd(nx, ny, (int) (linethickness/2), FALSE);
417     }   /* end for */
418 }  /* end HGArc */
419 
420 
421 /*----------------------------------------------------------------------------
422  * Routine:	RoundEnd (x, y, radius, filled_flag)
423  *
424  * Results:	Plots a filled (if requested) circle of the specified radius
425  *		centered about (x, y).
426  *----------------------------------------------------------------------------*/
427 
428 RoundEnd(x, y, radius, filled)
429 int x, y, radius;
430 int filled;
431 
432 {
433     double xs, ys, epsilon;
434     int i, j, k, extent, nx, ny;
435     int cx, cy;
436 
437 
438     if (radius < 1) {	/* too small to notice */
439         point(x, y);
440         return;
441     }
442 
443     xs = 0;
444     ys = radius;
445     epsilon = 1.0 / radius;
446     extent = pi * radius / 2.0;    /* 1/4 the circumference */
447 
448         /* Calculate the trajectory of the circle for 1/4 the circumference
449          * and mirror appropriately to get the other three quadrants.
450          */
451 
452     point(x, y+((int) ys));    /* take care if end of arc missed by */
453     point(x, y-((int) ys));    /* below formulation                 */
454     for (i=0; i<extent; ++i)
455     {
456              /* generate circumference */
457         xs += epsilon * ys;
458         nx = (int) (xs + x + 0.5);
459         if (nx < x) nx = x;  /* 1st quadrant, should be positive */
460         ys -= epsilon * xs;
461         ny = (int) (ys + y + 0.5);
462         if (ny < y) ny = y;  /* 1st quadrant, should be positive */
463 
464         if (filled == TRUE) {	/* fill from center */
465             cx = x;
466             cy = y;
467         } else {		/* fill from perimeter only (no fill) */
468             cx = nx;
469             cy = ny;
470         }
471         for (j=cx; j<=nx; ++j) {
472             for (k=cy; k<=ny; ++k) {
473                 point(j, k);
474                 point(j, 2*y-k);
475                 point(2*x-j, k);
476                 point(2*x-j, 2*y-k);
477             }  /* end for k */
478         }  /* end for j */;
479     }  /* end for i */;
480 }  /* end RoundEnd */;
481 
482 
483 /*----------------------------------------------------------------------------
484  * Routine:	Paramaterize (xpoints, ypoints, hparams, num_points)
485  *
486  * Results:	This routine calculates parameteric values for use in
487  *		calculating curves.  The parametric values are returned
488  *		in the array h.  The values are an approximation of
489  *		cumulative arc lengths of the curve (uses cord length).
490  *		For additional information, see paper cited below.
491  *----------------------------------------------------------------------------*/
492 
493 static Paramaterize(x, y, h, n)
494 int x[MAXPOINTS];
495 int y[MAXPOINTS];
496 float h[MAXPOINTS];
497 int n;
498 {
499 	register int dx;
500 	register int dy;
501 	register int i;
502 	register int j;
503 	float u[MAXPOINTS];
504 
505 
506 	for (i=1; i<=n; ++i) {
507 	    u[i] = 0;
508 	    for (j=1; j<i; j++) {
509 		dx = x[j+1] - x[j];
510 		dy = y[j+1] - y[j];
511 		u[i] += sqrt ((double) (dx * dx + dy * dy));
512 	    }
513 	}
514 	for (i=1; i<n; ++i)  h[i] = u[i+1] - u[i];
515 }  /* end Paramaterize */
516 
517 
518 /*----------------------------------------------------------------------------
519  * Routine:	PeriodicSpline (h, z, dz, d2z, d3z, npoints)
520  *
521  * Results:	This routine solves for the cubic polynomial to fit a
522  *		spline curve to the the points  specified by the list
523  *		of values.  The Curve generated is periodic.  The algorithms
524  *		for this curve are from the "Spline Curve Techniques" paper
525  *		cited below.
526  *----------------------------------------------------------------------------*/
527 
528 static PeriodicSpline(h, z, dz, d2z, d3z, npoints)
529 float h[MAXPOINTS];		/* paramaterization  */
530 int z[MAXPOINTS];		/* point list */
531 float dz[MAXPOINTS];			/* to return the 1st derivative */
532 float d2z[MAXPOINTS], d3z[MAXPOINTS];	/* 2nd and 3rd derivatives */
533 int npoints;				/* number of valid points */
534 {
535 	float d[MAXPOINTS];
536 	float deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS];
537 	float c[MAXPOINTS], r[MAXPOINTS], s[MAXPOINTS];
538 	int i;
539 
540 						/* step 1 */
541 	for (i=1; i<npoints; ++i) {
542 	    deltaz[i] = h[i] ? ((double) (z[i+1] - z[i])) / h[i] : 0;
543 	}
544 	h[0] = h[npoints-1];
545 	deltaz[0] = deltaz[npoints-1];
546 
547 						/* step 2 */
548 	for (i=1; i<npoints-1; ++i) {
549 	    d[i] = deltaz[i+1] - deltaz[i];
550 	}
551 	d[0] = deltaz[1] - deltaz[0];
552 
553 						/* step 3a */
554 	a[1] = 2 * (h[0] + h[1]);
555 	b[1] = d[0];
556 	c[1] = h[0];
557 	for (i=2; i<npoints-1; ++i) {
558 	    a[i] = 2*(h[i-1]+h[i]) - pow ((double) h[i-1],(double)2.0) / a[i-1];
559 	    b[i] = d[i-1] - h[i-1] * b[i-1]/a[i-1];
560 	    c[i] = -h[i-1] * c[i-1]/a[i-1];
561 	}
562 
563 						/* step 3b */
564 	r[npoints-1] = 1;
565 	s[npoints-1] = 0;
566 	for (i=npoints-2; i>0; --i) {
567 	    r[i] = -(h[i] * r[i+1] + c[i])/a[i];
568 	    s[i] = (6 * b[i] - h[i] * s[i+1])/a[i];
569 	}
570 
571 						/* step 4 */
572 	d2z[npoints-1] = (6 * d[npoints-2] - h[0] * s[1]
573 	                   - h[npoints-1] * s[npoints-2])
574 	                 / (h[0] * r[1] + h[npoints-1] * r[npoints-2]
575 	                    + 2 * (h[npoints-2] + h[0]));
576 	for (i=1; i<npoints-1; ++i) {
577 	    d2z[i] = r[i] * d2z[npoints-1] + s[i];
578 	}
579 	d2z[npoints] = d2z[1];
580 
581 						/* step 5 */
582 	for (i=1; i<npoints; ++i) {
583 	    dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i+1])/6;
584 	    d3z[i] = h[i] ? (d2z[i+1] - d2z[i])/h[i] : 0;
585 	}
586 }  /* end PeriodicSpline */
587 
588 
589 /*----------------------------------------------------------------------------
590  * Routine:	NaturalEndSpline (h, z, dz, d2z, d3z, npoints)
591  *
592  * Results:	This routine solves for the cubic polynomial to fit a
593  *		spline curve the the points  specified by the list of
594  *		values.  The alogrithms for this curve are from the
595  *		"Spline Curve Techniques" paper cited below.
596  *----------------------------------------------------------------------------*/
597 
598 static NaturalEndSpline(h, z, dz, d2z, d3z, npoints)
599 float h[MAXPOINTS];		/* parameterization */
600 int z[MAXPOINTS];		/* Point list */
601 float dz[MAXPOINTS];			/* to return the 1st derivative */
602 float d2z[MAXPOINTS], d3z[MAXPOINTS];	/* 2nd and 3rd derivatives */
603 int npoints;				/* number of valid points */
604 {
605 	float d[MAXPOINTS];
606 	float deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS];
607 	int i;
608 
609 						/* step 1 */
610 	for (i=1; i<npoints; ++i) {
611 	    deltaz[i] = h[i] ? ((double) (z[i+1] - z[i])) / h[i] : 0;
612 	}
613 	deltaz[0] = deltaz[npoints-1];
614 
615 						/* step 2 */
616 	for (i=1; i<npoints-1; ++i) {
617 	    d[i] = deltaz[i+1] - deltaz[i];
618 	}
619 	d[0] = deltaz[1] - deltaz[0];
620 
621 						/* step 3 */
622 	a[0] = 2 * (h[2] + h[1]);
623 	b[0] = d[1];
624 	for (i=1; i<npoints-2; ++i) {
625 	    a[i] = 2*(h[i+1]+h[i+2]) - pow((double) h[i+1],(double) 2.0)/a[i-1];
626 	    b[i] = d[i+1] - h[i+1] * b[i-1]/a[i-1];
627 	}
628 
629 						/* step 4 */
630 	d2z[npoints] = d2z[1] = 0;
631 	for (i=npoints-1; i>1; --i) {
632 	    d2z[i] = (6 * b[i-2] - h[i] *d2z[i+1])/a[i-2];
633 	}
634 
635 						/* step 5 */
636 	for (i=1; i<npoints; ++i) {
637 	    dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i+1])/6;
638 	    d3z[i] = h[i] ? (d2z[i+1] - d2z[i])/h[i] : 0;
639 	}
640 }  /* end NaturalEndSpline */
641 
642 
643 /*----------------------------------------------------------------------------
644  * Routine:	HGCurve(xpoints, ypoints, num_points)
645  *
646  * Results:	This routine generates a smooth curve through a set of points.
647  *		The method used is the parametric spline curve on unit knot
648  *		mesh described in "Spline Curve Techniques" by Patrick
649  *		Baudelaire, Robert Flegal, and Robert Sproull -- Xerox Parc.
650  *----------------------------------------------------------------------------*/
651 
652 #define PointsPerInterval 32
653 
654 HGCurve(x, y, numpoints)
655 int x[MAXPOINTS];
656 int y[MAXPOINTS];
657 int numpoints;
658 {
659 	float h[MAXPOINTS], dx[MAXPOINTS], dy[MAXPOINTS];
660 	float d2x[MAXPOINTS], d2y[MAXPOINTS], d3x[MAXPOINTS], d3y[MAXPOINTS];
661 	float t, t2, t3;
662 	register int j;
663 	register int k;
664 	register int nx;
665 	register int ny;
666 	int lx, ly;
667 
668 
669 	lx = x[1];
670 	ly = y[1];
671 
672 	     /* Solve for derivatives of the curve at each point
673               * separately for x and y (parametric).
674 	      */
675 	Paramaterize(x, y, h, numpoints);
676 							/* closed curve */
677 	if ((x[1] == x[numpoints]) && (y[1] == y[numpoints])) {
678 	    PeriodicSpline(h, x, dx, d2x, d3x, numpoints);
679 	    PeriodicSpline(h, y, dy, d2y, d3y, numpoints);
680 	} else {
681 	    NaturalEndSpline(h, x, dx, d2x, d3x, numpoints);
682 	    NaturalEndSpline(h, y, dy, d2y, d3y, numpoints);
683 	}
684 
685 	      /* generate the curve using the above information and
686 	       * PointsPerInterval vectors between each specified knot.
687 	       */
688 
689 	for (j=1; j<numpoints; ++j) {
690 	    if ((x[j] == x[j+1]) && (y[j] == y[j+1])) continue;
691 	    for (k=0; k<=PointsPerInterval; ++k) {
692 		t = (float) k * h[j] / (float) PointsPerInterval;
693 		t2 = t * t;
694 		t3 = t * t * t;
695 		nx = x[j] + (int) (t * dx[j] + t2 * d2x[j]/2 + t3 * d3x[j]/6);
696 		ny = y[j] + (int) (t * dy[j] + t2 * d2y[j]/2 + t3 * d3y[j]/6);
697 		HGtline(lx, ly, nx, ny);
698 		lx = nx;
699 		ly = ny;
700 	    }  /* end for k */
701 	}  /* end for j */
702 }  /* end HGCurve */
703 
704 
705 /*----------------------------------------------------------------------------
706  * Routine:	HGtline(xstart, ystart, xend, yend)
707  *
708  * Results:	Draws a line of proper thickness by calling "line" numerous
709  *		times until the desired thickness is reached.
710  *----------------------------------------------------------------------------*/
711 
712 HGtline(x0, y0, x1, y1)
713 register int x0;
714 register int y0;
715 int x1;
716 int y1;
717 {
718         register int xs;
719 	register int xe;
720 	register int ys;
721 	register int ye;
722         double morelen, theta, wx, wy, xx, xy;
723         int addln, j, xdir, ydir, dx, dy;
724 
725 
726         xdir = ydir = 1;
727         dx = x1 - x0;		/* calculate direction to move to  */
728         dy = y1 - y0;		/* move to draw additional lines if needed */
729         if (dx < 0 ) {		/* for extra thickness */
730             dx = -dx;
731             xdir = -1;
732         }
733         if (dy < 0 ) {
734             dy = -dy;
735             ydir = -1;
736         }
737 
738         morelen = linethickness / 2;
739 	addln = (int) morelen;
740         RoundEnd (x0, y0, (int) morelen, TRUE);    /* add rounded end */
741 
742         for (j=(-addln); j<=addln; ++j) {
743 	    if (dy == 0) {
744 		xs = x0;
745 		xe = x1;
746 		ys = ye = y0 + j;
747 	    }
748 	    if (dx == 0) {
749 		ys = y0;
750 		ye = y1;
751 		xs = xe = x0 + j;
752 	    }
753 	    if ((dx != 0) && (dy != 0)) {
754 		theta =  pi / 2.0 - atan( ((double) dx)/((double) dy) );
755 		wx = j * sin(theta);
756 		wy = j * cos(theta);
757 		xs = x0 + wx * xdir;
758 		ys = y0 - wy * ydir;
759 		xe = x1 + wx * xdir;
760 		ye = y1 - wy * ydir;
761 	    }
762 	    line(xs, ys, xe, ye);
763         }  /* end for */
764 
765         RoundEnd(x1, y1, (int) morelen, TRUE);    /* add rounded end */
766 }  /* end HGtline */
767