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