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