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