1 /*
2 * FIG : Facility for Interactive Generation of figures
3 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
4 * Parts Copyright (c) 1989-2007 by Brian V. Smith
5 * Parts Copyright (c) 1991 by Paul King
6 *
7 * Any party obtaining a copy of these files is granted, free of charge, a
8 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
9 * nonexclusive right and license to deal in this software and documentation
10 * files (the "Software"), including without limitation the rights to use,
11 * copy, modify, merge, publish, distribute, sublicense and/or sell copies of
12 * the Software, and to permit persons who receive copies from any such
13 * party to do so, with the only requirement being that the above copyright
14 * and this permission notice remain intact.
15 *
16 */
17
18 /*
19 * Routines dealing with geometry under the following headings:
20 * COMPUTE NORMAL, CLOSE TO VECTOR, COMPUTE ARC CENTER,
21 * COMPUTE ANGLE, COMPUTE DIRECTION, LATEX LINE ROUTINES.
22 */
23
24 #include "fig.h"
25 #include "resources.h"
26 #include "object.h"
27 #include "u_geom.h"
28
29 /*************************** ROTATE VECTOR **********************
30 In fact, rotate & scale ;-) */
31
32
33 int compute_arccenter (F_pos p1, F_pos p2, F_pos p3, float *x, float *y);
34 int gcd (int a, int b);
35
36 static void
rotate_vector(double * vx,double * vy,double c,double s)37 rotate_vector(double *vx, double *vy, double c, double s)
38 {
39 double wx, wy;
40 wx = c * (*vx) - s * (*vy);
41 wy = s * (*vx) + c * (*vy);
42 *vx = wx;
43 *vy = wy;
44 }
45
46
47 /*************************** COMPUTE NORMAL **********************
48
49 Input arguments :
50 (x1,y1)(x2,y2) : the vector
51 direction : direction of the normal vector to (x1,y1)(x2,y2)
52 Output arguments :
53 (*x,*y)(x2,y2) : a normal vector.
54 Return value : none
55
56 ******************************************************************/
57
compute_normal(float x1,float y1,int x2,int y2,int direction,int * x,int * y)58 void compute_normal(float x1, float y1, int x2, int y2, int direction, int *x, int *y)
59 {
60 if (direction) { /* counter clockwise */
61 *x = round(x2 - (y2 - y1));
62 *y = round(y2 - (x1 - x2));
63 } else {
64 *x = round(x2 + (y2 - y1));
65 *y = round(y2 + (x1 - x2));
66 }
67 }
68
69 /******************** CLOSE TO VECTOR **************************
70
71 Input arguments:
72 (x1,y1)(x2,y2) : the vector
73 (xp,yp) : the point
74 d : tolerance (max. allowable distance from the point to the vector)
75 dd : d * d
76 Output arguments:
77 (*px,*py) : a point on the vector which is not far from (xp,yp)
78 by more than d. Normally the vector (*px,*py)(xp,yp)
79 is normal to vector (x1,y1)(x2,y2) except when (xp,yp)
80 is within d from (x1,y1) or (x2,y2), in which cases,
81 (*px,*py) = (x1,y1) or (x2,y2) respectively.
82 Return value :
83 0 : No point on the vector is within d from (xp, yp)
84 1 : (*px, *py) is such a point.
85
86 ******************************************************************/
87
close_to_vector(int x1,int y1,int x2,int y2,int xp,int yp,int d,float dd,int * px,int * py)88 int close_to_vector(int x1, int y1, int x2, int y2, int xp, int yp, int d, float dd, int *px, int *py)
89 {
90 int xmin, ymin, xmax, ymax;
91 float x, y, slope, D2, dx, dy;
92
93 if (abs(xp - x1) <= d && abs(yp - y1) <= d) {
94 *px = x1;
95 *py = y1;
96 return (1);
97 }
98 if (abs(xp - x2) <= d && abs(yp - y2) <= d) {
99 *px = x2;
100 *py = y2;
101 return (1);
102 }
103 if (x1 < x2) {
104 xmin = x1 - d;
105 xmax = x2 + d;
106 } else {
107 xmin = x2 - d;
108 xmax = x1 + d;
109 }
110 if (xp < xmin || xmax < xp)
111 return (0);
112
113 if (y1 < y2) {
114 ymin = y1 - d;
115 ymax = y2 + d;
116 } else {
117 ymin = y2 - d;
118 ymax = y1 + d;
119 }
120 if (yp < ymin || ymax < yp)
121 return (0);
122
123 if (x2 == x1) {
124 x = x1;
125 y = yp;
126 } else if (y1 == y2) {
127 x = xp;
128 y = y1;
129 } else {
130 slope = ((float) (x2 - x1)) / ((float) (y2 - y1));
131 y = (slope * (xp - x1 + slope * y1) + yp) / (1 + slope * slope);
132 x = ((float) x1) + slope * (y - y1);
133 }
134 dx = ((float) xp) - x;
135 dy = ((float) yp) - y;
136 D2 = dx * dx + dy * dy;
137 if (D2 < dd) {
138 *px = round(x);
139 *py = round(y);
140 return (1);
141 }
142 return (0);
143 }
144
145 /********************* COMPUTE ARC RADIUS ******************
146 Input arguments :
147 (x1, y1), (x2, y2), (x3, y3): 3 points on arc
148 Output arguments :
149 (*r) : Radius of the arc
150 Return value :
151 0 : if p1, p2, p3 are co-linear.
152 1 : if they are not.
153 *************************************************************/
154
155 int
compute_arcradius(int x1,int y1,int x2,int y2,int x3,int y3,float * r)156 compute_arcradius(int x1, int y1, int x2, int y2, int x3, int y3, float *r)
157 {
158 F_pos p1, p2, p3;
159 float cx, cy, dx, dy;
160
161 p1.x = x1; p1.y = y1;
162 p2.x = x2; p2.y = y2;
163 p3.x = x3; p3.y = y3;
164
165 if (!compute_arccenter(p1, p2, p3, &cx, &cy))
166 return 0;
167 dx = x2 - cx;
168 dy = y2 - cy;
169 if (dx == 0.0 && dy == 0.0) {
170 *r = 0.0;
171 return 0;
172 }
173 else {
174 *r = sqrt(dx * dx + dy * dy);
175 return 1;
176 }
177 }
178
179
180 /********************* COMPUTE ARC CENTER ******************
181
182 New routine 12/11/00 from Thomas Henlich - better at finding center
183
184 Input arguments :
185 p1, p2, p3 : 3 points on the arc
186 Output arguments :
187 (*x,*y) : Center of the arc
188 Return value :
189 0 : if p1, p2, p3 are co-linear.
190 1 : if they are not.
191
192 *************************************************************/
193
194 int
compute_arccenter(F_pos p1,F_pos p2,F_pos p3,float * x,float * y)195 compute_arccenter(F_pos p1, F_pos p2, F_pos p3, float *x, float *y)
196 {
197 double a, b, c, d, e, f, g, h, i, j;
198 double resx, resy;
199
200 if ((p1.x == p3.x && p1.y == p3.y) ||
201 (p1.x == p2.x && p1.y == p2.y) ||
202 (p2.x == p3.x && p2.y == p3.y)) {
203 return 0;
204 }
205 a = (double)p1.x - (double)p2.x;
206 b = (double)p1.x + (double)p2.x;
207 c = (double)p1.y - (double)p2.y;
208 d = (double)p1.y + (double)p2.y;
209 e = (a*b + c*d)/2.0;
210
211 f = (double)p2.x - (double)p3.x;
212 g = (double)p2.x + (double)p3.x;
213 h = (double)p2.y - (double)p3.y;
214 i = (double)p2.y + (double)p3.y;
215 j = (f*g + h*i)/2.0;
216
217 if (a*h - c*f != 0.0)
218 resy = (a*j - e*f)/(a*h - c*f);
219 else {
220 return 0;
221 }
222 if (a != 0.0)
223 resx = (e - resy*c)/a;
224 else if (f != 0.0)
225 resx = (j - resy*h)/f;
226 else {
227 return 0;
228 }
229
230 *x = (float) resx;
231 *y = (float) resy;
232
233 return 1;
234 }
235
236 /********************* CLOSE TO ARC ************************
237
238 Input arguments :
239 a: arc object
240 xp, yp: Point
241 d: tolerance
242 Output arguments :
243 (*px,*py) : `near' point on the arc
244 Return value :
245 0 on failure, 1 on success
246
247 *************************************************************/
248
249 int
close_to_arc(F_arc * a,int xp,int yp,int d,float * px,float * py)250 close_to_arc(F_arc *a, int xp, int yp, int d, float *px, float *py)
251 {
252 int i, ok;
253 double ux, uy, vx, vy, wx, wy, dx, dy, rarc, rp, uang, vang, wang, pang;
254 double cphi, sphi;
255
256 ux = a->point[0].x - a->center.x;
257 uy = a->point[0].y - a->center.y;
258 vx = a->point[1].x - a->center.x;
259 vy = a->point[1].y - a->center.y;
260 wx = a->point[2].x - a->center.x;
261 wy = a->point[2].y - a->center.y;
262 if (((ux == 0.0) && (uy == 0.0)) || ((vx == 0.0) && (vy == 0.0)) ||
263 ((wx == 0.0) && (wy == 0.0)))
264 return 0;
265 rarc = sqrt(ux * ux + uy * uy);
266 dx = xp - a->center.x;
267 dy = yp - a->center.y;
268 if ((dx == 0.0) && (dy == 0.0))
269 return 0;
270 rp = sqrt(dx * dx + dy * dy);
271 if (fabs(rarc - rp) > d)
272 return 0;
273 /* radius ok, check special cases first */
274 for (i = 0; i < 3; i++) {
275 if (abs(xp - a->point[i].x) <= d && abs(yp - a->point[i].y <= d)) {
276 *px = (float)(a->point[i].x);
277 *py = (float)(a->point[i].y);
278 return 1;
279 }
280 }
281
282 /* check angles */
283 uang = atan2(uy, ux);
284 vang = atan2(vy, vx);
285 wang = atan2(wy, wx);
286 pang = atan2(dy, dx);
287
288 if (uang >= wang) {
289 if ((uang >= vang) && (vang >= wang))
290 ok = (uang >= pang) && (pang >= wang);
291 else {
292 wang += 2*M_PI;
293 if (uang >= pang)
294 pang += 2*M_PI;
295 ok = (pang <= wang);
296 }
297 }
298 else {
299 if ((uang <= vang) && (vang <= wang))
300 ok = (uang <= pang) && (pang <= wang);
301 else {
302 wang -= 2*M_PI;
303 if (uang <= pang)
304 pang -= 2*M_PI;
305 ok = (pang >= wang);
306 }
307 }
308 if (!ok)
309 return 0;
310 /* generate point */
311 cphi = cos(pang);
312 sphi = sin(pang);
313 dx = rarc;
314 dy = 0;
315 rotate_vector(&dx, &dy, cphi, sphi);
316 *px = (float)(a->center.x + dx);
317 *py = (float)(a->center.y + dy);
318 return 1;
319 }
320
321 /********************* COMPUTE ELLIPSE DISTANCE ***************
322
323 Input arguments :
324 (dx, dy) vector from center to test point
325 (a, b) main axes (`radii') of ellipse
326 Output arguments :
327 *dis distance from center as is
328 *r distance from center as should be
329 *************************************************************/
330 static void
compute_ellipse_distance(double dx,double dy,double a,double b,double * dis,double * r)331 compute_ellipse_distance(double dx, double dy, double a, double b, double *dis, double *r)
332 {
333 if (dx == 0.0 && dy == 0.0)
334 *dis = 0.0; /* prevent core dumps */
335 else
336 *dis = sqrt(dx * dx + dy * dy);
337 if (a * dy == 0.0 && b * dx == 0.0)
338 *r = 0.0; /* prevent core dumps */
339 else
340 *r = a * b * (*dis) / sqrt(1.0 * b * b * dx * dx + 1.0 * a * a * dy * dy);
341 }
342
343 /********************* CLOSE TO ELLIPSE ************************
344
345 Input arguments :
346 e: (possibly rotated!) ellipse object
347 xp, yp: Point
348 d: tolerance
349 Output arguments :
350 (*ex,*ey) : `near' point on the ellipse
351 (*vx,*vy) : tangent vector
352 Return value :
353 0 on failure, 1 on success
354 *************************************************************/
355 int
close_to_ellipse(F_ellipse * e,int xp,int yp,int d,float * ex,float * ey,float * vx,float * vy)356 close_to_ellipse(F_ellipse *e, int xp, int yp, int d, float *ex, float *ey, float *vx, float *vy)
357 {
358 double a, b, dx, dy, phi, cphi, sphi, qx, qy, dvx, dvy;
359 double dis, r;
360 int wx, wy;
361
362 dx = xp - e->center.x;
363 dy = yp - e->center.y;
364 a = e->radiuses.x;
365 b = e->radiuses.y;
366 phi = e->angle;
367
368 if (phi == 0.0) { /* extra case for speed */
369 compute_ellipse_distance(dx, dy, a, b, &dis, &r);
370 if (fabs(dis - r) <= d) {
371 wx = round(r * dx / dis);
372 wy = round(r * dy / dis);
373 *ex = (float)(e->center.x + wx);
374 *ey = (float)(e->center.y + wy);
375 if (a == 0.0) { /* degenerated */
376 dvx = 0.0;
377 dvy = b;
378 }
379 else if (b == 0.0) { /* degenerated */
380 dvx = a;
381 dvy = 0.0;
382 }
383 else {
384 dvx = - ((double)wy) * a /b;
385 dvy = ((double)wx) * b /a;
386 }
387 *vx = (float)dvx;
388 *vy = (float)dvy;
389 return 1;
390 }
391 }
392 else { /* general case: rotation */
393 cphi = cos(phi);
394 sphi = sin(phi);
395 rotate_vector(&dx, &dy, cphi, sphi);
396 compute_ellipse_distance(dx, dy, a, b, &dis, &r);
397 if (fabs(dis - r) <= d) {
398 qx = r * dx / dis;
399 qy = r * dy / dis;
400 if (a == 0.0) {
401 dvx = 0.0;
402 dvy = b;
403 }
404 else if (b == 0.0) {
405 dvx = a;
406 dvy = 0.0;
407 }
408 else {
409 dvx = -qy * a/b;
410 dvy = qx * b/a;
411 }
412 rotate_vector(&qx, &qy, cphi, -sphi);
413 rotate_vector(&dvx, &dvy, cphi, -sphi);
414 *ex = (float)(e->center.x + qx);
415 *ey = (float)(e->center.y + qy);
416 *vx = (float)dvx;
417 *vy = (float)dvy;
418 return 1;
419 }
420 }
421 return 0;
422 }
423
424
425 /********************* CLOSE TO POLYLINE ************************
426
427 Input arguments :
428 p: polyline object
429 xp, yp: Point
430 d: tolerance
431 sd: singular tolerance
432 Output arguments :
433 (*px,*py) : `near' point on the polyline
434 (*lx1, *ly1), (*lx2, *ly2): neighbor polyline corners
435 Return value :
436 0 on failure, 1 on success
437 *************************************************************/
438
439 int
close_to_polyline(F_line * l,int xp,int yp,int d,int sd,int * px,int * py,int * lx1,int * ly1,int * lx2,int * ly2)440 close_to_polyline(F_line *l, int xp, int yp, int d, int sd, int *px, int *py, int *lx1, int *ly1, int *lx2, int *ly2)
441 {
442 F_point *point;
443 int x1, y1, x2, y2;
444 float tol2;
445 tol2 = (float) d*d;
446
447 point = l->points;
448 x1 = point->x;
449 y1 = point->y;
450 if (abs(xp - x1) <= sd && abs(yp - y1) <= sd) {
451 *px = *lx1 = *lx2 = x1;
452 *py = *ly1 = *ly2 = y1;
453 return 1;
454 }
455 for (point = point->next; point != NULL; point = point->next) {
456 x2 = point->x;
457 y2 = point->y;
458 if (abs(xp - x2) <= sd && abs(yp - y2) <= sd) {
459 *px = *lx1 = *lx2 = x2;
460 *py = *ly1 = *ly2 = y2;
461 return 1;
462 }
463 if (close_to_vector(x1, y1, x2, y2, xp, yp, d, tol2, px, py)) {
464 *lx1 = x1; *ly1 = y1;
465 *lx2 = x2; *ly2 = y2;
466 return 1;
467 }
468 x1 = x2;
469 y1 = y2;
470 }
471 return 0;
472 }
473
474 /********************* CLOSE TO SPLINE ************************
475
476 Input arguments :
477 p: spline object
478 xp, yp: Point
479 d: tolerance
480 Output arguments :
481 (*px,*py) : `near' point on the spline
482 (*lx1, *ly1), (*lx2, *ly2): points for tangent
483 Return value :
484 0 on failure, 1 on success
485 *************************************************************/
486
487
488
489 /* Take a look at close_to_vector in u_geom.c ... and then tell me:
490 Are those people mathematicians? */
491
492 static int
close_to_float_vector(float x1,float y1,float x2,float y2,float xp,float yp,float d,float * px,float * py,float * lambda)493 close_to_float_vector(float x1, float y1, float x2, float y2, float xp, float yp, float d, float *px, float *py, float *lambda)
494 {
495 /* always stores in (*px,*py) the orthogonal projection of (xp,yp)
496 onto the line through (x1,y1) and (x2,y2)
497 returns 1 if (x,y) is closer than d to the *segment* between these two,
498 otherwise 0 */
499 float rx, ry, radius, radius2, nx, ny, distance;
500 int ok;
501
502 rx = x2 - x1;
503 ry = y2 - y1; /* direction vector of the straight line */
504 radius2 = rx * rx + ry * ry;
505 if (radius2 > 0.0)
506 radius = sqrt(radius2); /* its length */
507 else
508 return(0); /* singular case */
509 nx = -ry/radius;
510 ny = rx/radius; /* unit normal vector */
511 distance = (xp - x1) * nx + (yp - y1) * ny; /* distance (Hesse) */
512 *px = xp - distance * nx;
513 *py = yp - distance * ny; /* the projection */
514 if (fabs(distance) > d)
515 ok = 0;
516 else if ((xp-x1) * (xp-x1) + (yp-y1) * (yp-y1) <= d * d) {
517 *lambda = 0.0;
518 ok = 1; /* close to endpoint */
519 }
520 else if ((xp-x2) * (xp-x2) + (yp-y2) * (yp-y2) <= d * d) {
521 *lambda = 1.0;
522 ok = 1; /* close to endpoint */
523 }
524 else {
525 *lambda = ((xp-x1)*rx + (yp-y1)*ry)/radius2; /* the parameter (0=point1, 1=point2) */
526 if ((0.0 <= *lambda) && (*lambda <= 1.0))
527 ok = 1;
528 else
529 ok = 0;
530 }
531 return ok;
532 }
533
534 static int isfirst, prevx, prevy, firstx, firsty;
535 static int tx, ty, td, px1, py1, px2, py2, foundx, foundy;
536
537 /* dirty trick: redefine init_point_array to do nothing and
538 add_point to check distance */
539
540 /**********************************/
541 /* include common spline routines */
542 /**********************************/
543
544 static Boolean add_point(int x, int y);
545 static void init_point_array(void);
546 static Boolean add_closepoint(void);
547
548 #include "u_draw_spline.c"
549
550 static void
init_point_array(void)551 init_point_array(void)
552 {
553 }
554
555 static Boolean
add_point(int x,int y)556 add_point(int x, int y)
557 {
558 float ux, uy, lambda;
559
560 if (isfirst) {
561 firstx = x;
562 firsty = y;
563 isfirst = False;
564 } else {
565 if (close_to_float_vector((float)prevx, (float)prevy, (float)x, (float)y,
566 (float)tx, (float)ty, (float)td, &ux, &uy, &lambda)) {
567 foundx = round(ux);
568 foundy = round(uy);
569 px1 = prevx;
570 py1 = prevy;
571 px2 = x;
572 py2 = y;
573 DONE = True;
574 }
575 }
576 prevx = x;
577 prevy = y;
578
579 return True;
580 }
581
582 static Boolean
add_closepoint(void)583 add_closepoint(void) {
584 return add_point(firstx, firsty);
585 }
586
587 int
close_to_spline(F_spline * spline,int xp,int yp,int d,int * px,int * py,int * lx1,int * ly1,int * lx2,int * ly2)588 close_to_spline(F_spline *spline, int xp, int yp, int d, int *px, int *py, int *lx1, int *ly1, int *lx2, int *ly2)
589 {
590 float precision;
591
592 tx = xp; ty = yp; /* test point */
593 td = d; /* tolerance */
594 isfirst = True;
595 precision = HIGH_PRECISION;
596
597 if (open_spline(spline))
598 compute_open_spline(spline, precision);
599 else
600 compute_closed_spline(spline, precision);
601 if (DONE) {
602 *px = foundx;
603 *py = foundy;
604 *lx1 = px1;
605 *ly1 = py1;
606 *lx2 = px2;
607 *ly2 = py2;
608 }
609 return DONE;
610 }
611
612 /********************* COMPUTE ANGLE ************************
613
614 Input arguments :
615 (dx,dy) : the vector (0,0)(dx,dy)
616 Output arguments : none
617 Return value : the angle of the vector in the range [0, 2PI)
618
619 *************************************************************/
620
621 double
compute_angle(double dx,double dy)622 compute_angle(double dx, double dy) /* compute the angle between 0 to 2PI */
623 {
624 double alpha;
625
626 if (dx == 0) {
627 if (dy > 0)
628 alpha = M_PI_2;
629 else
630 alpha = 3 * M_PI_2;
631 } else if (dy == 0) {
632 if (dx > 0)
633 alpha = 0;
634 else
635 alpha = M_PI;
636 } else {
637 alpha = atan(dy / dx); /* range = -PI/2 to PI/2 */
638 if (dx < 0)
639 alpha += M_PI;
640 else if (dy < 0)
641 alpha += M_2PI;
642 }
643 return (alpha);
644 }
645
646
647 /********************* COMPUTE DIRECTION ********************
648
649 Input arguments :
650 p1, p2, p3 : 3 points of an arc with p1 the first and p3 the last.
651 Output arguments : none
652 Return value :
653 0 : if the arc passes p1, p2 and p3 (in that order) in
654 clockwise direction
655 1 : if direction is counterclockwise
656
657 *************************************************************/
658
659 int
compute_direction(F_pos p1,F_pos p2,F_pos p3)660 compute_direction(F_pos p1, F_pos p2, F_pos p3)
661 {
662 double diff, dx, dy, alpha, theta;
663
664 dx = p2.x - p1.x;
665 dy = p1.y - p2.y; /* because origin of the screen is on the
666 * upper left corner */
667
668 alpha = compute_angle(dx, dy);
669
670 dx = p3.x - p2.x;
671 dy = p2.y - p3.y;
672 theta = compute_angle(dx, dy);
673
674 diff = theta - alpha;
675 if ((0 < diff && diff < M_PI) || diff < -M_PI) {
676 return (1); /* counterclockwise */
677 }
678 return (0); /* clockwise */
679 }
680
681 /********************* COMPUTE ANGLE BY 3 POINTS ************
682 Input arguments:
683 p1, p2, p3 : Points
684 Output arguments : angle in the range [0, 2PI)
685 Return value : 0 on failure, 1 on success
686 *************************************************************/
687
688 int
compute_3p_angle(F_point * p1,F_point * p2,F_point * p3,double * alpha)689 compute_3p_angle(F_point *p1, F_point *p2, F_point *p3, double *alpha)
690 {
691 double vx, vy, wx, wy;
692 double dang;
693
694 vx = (double)(p1->x - p2->x);
695 vy = (double)(p1->y - p2->y);
696 wx = (double)(p3->x - p2->x);
697 wy = (double)(p3->y - p2->y);
698 if (((vx == 0.0) && (vy == 0.0)) || ((wx == 0.0) && (wy == 0.0))) {
699 return(0);
700 }
701 dang = -atan2(wy, wx) + atan2(vy, vx); /* difference angle */
702 if (dang < 0.0)
703 dang += 2*M_PI;
704 *alpha = dang;
705 return(1);
706 }
707
708
709 /********************* COMPUTE LINE ANGLE *******************
710
711 Input arguments :
712 l : the line object
713 p : point
714 Output arguments : angle in the range [0, 2PI)
715 Return value : 0 on failure, 1 on success
716 *************************************************************/
717
718 int
compute_line_angle(F_line * l,F_point * p,double * alpha)719 compute_line_angle(F_line *l, F_point *p, double *alpha)
720 {
721 F_point *q, *vorq, *nachq;
722
723 q = l->points;
724 vorq = NULL;
725 nachq = l->points->next;
726
727 for (; q != NULL; vorq = q, q = nachq, nachq = (nachq) ? nachq->next : NULL) {
728 if ((q->x == p->x) && (q->y == p->y) && (vorq)) {
729 if (nachq) {
730 if (compute_3p_angle(vorq, q, nachq, alpha))
731 return(1);
732 }
733 else if (l->type == T_POLYGON) {
734 if (compute_3p_angle(vorq, q, l->points->next, alpha))
735 return(1);
736 }
737 }
738 }
739 return(0);
740 }
741
742 /********************* COMPUTE ARC ANGLE *******************
743
744 Input arguments :
745 a : the arc
746 Output arguments : orientated angle
747 Return value : 0 on failure, 1 on success
748 *************************************************************/
749
750 int
compute_arc_angle(F_arc * a,double * alpha)751 compute_arc_angle(F_arc *a, double *alpha)
752 {
753 double ux, uy, vx, vy, wx, wy;
754 double uang, vang, wang, dang;
755
756 ux = a->point[0].x - a->center.x;
757 uy = a->point[0].y - a->center.y;
758 vx = a->point[1].x - a->center.x;
759 vy = a->point[1].y - a->center.y;
760 wx = a->point[2].x - a->center.x;
761 wy = a->point[2].y - a->center.y;
762 if (((ux == 0.0) && (uy == 0.0)) || ((vx == 0.0) && (vy == 0.0)) ||
763 ((wx == 0.0) && (wy == 0.0))) {
764 return(0);
765 }
766 uang = atan2(uy, ux);
767 vang = atan2(vy, vx);
768 wang = atan2(wy, wx);
769
770 if (uang >= wang) {
771 if ((uang >= vang) && (vang >= wang))
772 dang = uang - wang;
773 else
774 dang = -2*M_PI + uang - wang;
775 }
776 else {
777 if ((uang <= vang) && (vang <= wang))
778 dang = uang - wang;
779 else
780 dang = 2*M_PI + uang - wang;
781 }
782 *alpha = dang;
783 return(1);
784 }
785
786 /********************* COMPUTE POLY LENGTH *******************
787
788 Input arguments :
789 l : the polyline
790 Output arguments :
791 *lp : computed length
792 Return value : 0 on failure, 1 on success
793 *************************************************************/
794 int
compute_poly_length(F_line * l,float * lp)795 compute_poly_length(F_line *l, float *lp)
796 {
797 double length;
798 F_point *q, *prevq;
799 double dx, dy, dl;
800
801 length = 0.0;
802 prevq = NULL;
803
804 for (q = l->points; q != NULL; prevq = q, q = q->next) {
805 if (prevq) {
806 dx = (double)(q->x - prevq->x);
807 dy = (double)(q->y - prevq->y);
808 if ((dx == 0.0) && (dy == 0.0))
809 dl = 0.0;
810 else
811 dl = sqrt(dx * dx + dy * dy);
812 length += dl;
813 }
814 }
815 *lp = (float)length;
816 return 1;
817 }
818
819 /********************* COMPUTE ARC LENGTH *******************
820
821 Input arguments :
822 a : the arc
823 Output arguments :
824 *lp : computed length
825 Return value : 0 on failure, 1 on success
826 *************************************************************/
827 int
compute_arc_length(F_arc * a,float * lp)828 compute_arc_length(F_arc *a, float *lp)
829 {
830 double alpha, ux, uy, arcradius;
831
832 if (!compute_arc_angle(a, &alpha))
833 return 0;
834 ux = a->point[0].x - a->center.x;
835 uy = a->point[0].y - a->center.y;
836 if ((ux == 0.0) && (uy == 0.0))
837 arcradius = 0.0;
838 else
839 arcradius = sqrt(ux*ux + uy*uy);
840 *lp = (float)(arcradius * fabs(alpha));
841 return 1;
842 }
843
844 /********************* COMPUTE POLYGON AREA *******************
845
846 Input arguments :
847 l : the polygon
848 Output arguments :
849 *lp : computed area
850 Return value : 0 on failure, 1 on success
851 *************************************************************/
852
853 void
compute_poly_area(F_line * l,float * ap)854 compute_poly_area(F_line *l, float *ap)
855 {
856 float area;
857 F_point *q, *prevq;
858
859 area = 0.0;
860 prevq = NULL;
861 for (q = l->points; q != NULL; prevq = q, q = q->next) {
862 if (prevq)
863 area += (float)-prevq->x * (float)q->y + (float)prevq->y * (float)q->x;
864 }
865 area /= 2.0;
866 *ap = area;
867 }
868
869 /********************* COMPUTE ARC AREA *******************
870
871 Input arguments :
872 a : the arc
873 Output arguments :
874 *lp : computed area
875 Return value : 0 on failure, 1 on success
876 *************************************************************/
877
878 int
compute_arc_area(F_arc * a,float * ap)879 compute_arc_area(F_arc *a, float *ap)
880 {
881 double area, tri_area, alpha, ux, uy, arcradius;
882
883 if (!compute_arc_angle(a, &alpha))
884 return 0;
885 ux = a->point[0].x - a->center.x;
886 uy = a->point[0].y - a->center.y;
887 if ((ux == 0.0) && (uy == 0.0))
888 arcradius = 0.0;
889 else
890 arcradius = sqrt(ux*ux + uy*uy);
891 /* sector (`pie-shape')
892 A = (pi * r * r) * (alpha / 2pi) */
893 area = 0.5 * alpha * arcradius * arcradius;
894 /* open arc: subtract triangular area */
895 if (a->type == T_OPEN_ARC) {
896 tri_area = 0.5 * arcradius * arcradius * fabs(sin(alpha));
897 if (area > 0)
898 area -= tri_area;
899 else
900 area += tri_area;
901 }
902 *ap = (float)area;
903 return 1;
904 }
905
906 /********************* COMPUTE ELLIPSE AREA *******************
907
908 Input arguments :
909 e : the ellipse
910 Output arguments :
911 *lp : computed area
912 Return value : 0 on failure, 1 on success
913 *************************************************************/
914
915 int
compute_ellipse_area(F_ellipse * e,float * ap)916 compute_ellipse_area(F_ellipse *e, float *ap)
917 {
918 *ap = (float)(M_PI * e->radiuses.x * e->radiuses.y);
919 return 1;
920 }
921
922
923 /*********************** LATEX LINE ROUTINES ***************************/
924
925 /*
926 * compute greatest common divisor, assuming 0 < a <= b
927 */
928
929 int
pgcd(int a,int b)930 pgcd(int a, int b)
931 {
932 b = b % a;
933 return (b) ? gcd(b, a) : a;
934 }
935
936 /*
937 * compute greatest common divisor
938 */
939
940 int
gcd(int a,int b)941 gcd(int a, int b)
942 {
943 if (a < 0)
944 a = -a;
945 if (b < 0)
946 b = -b;
947 return (a <= b) ? pgcd(a, b) : pgcd(b, a);
948 }
949
950 /*
951 * compute least common multiple
952 */
953
954 int
lcm(int a,int b)955 lcm(int a, int b)
956 {
957 return abs(a * b) / gcd(a, b);
958 }
959
960
961 double rad2deg = 57.295779513082320877;
962
963 struct angle_table {
964 int x, y;
965 double angle;
966 };
967
968 struct angle_table line_angles[25] =
969 {{0, 1, 90.0},
970 {1, 0, 0.0},
971 {1, 1, 45.0},
972 {1, 2, 63.434948822922010648},
973 {1, 3, 71.565051177077989351},
974 {1, 4, 75.963756532073521417},
975 {1, 5, 78.690067525979786913},
976 {1, 6, 80.537677791974382609},
977 {2, 1, 26.565051177077989351},
978 {2, 3, 56.309932474020213086},
979 {2, 5, 68.198590513648188229},
980 {3, 1, 18.434948822922010648},
981 {3, 2, 33.690067525979786913},
982 {3, 4, 53.130102354155978703},
983 {3, 5, 59.036243467926478582},
984 {4, 1, 14.036243467926478588},
985 {4, 3, 36.869897645844021297},
986 {4, 5, 51.340191745909909396},
987 {5, 1, 11.309932474020213086},
988 {5, 2, 21.801409486351811770},
989 {5, 3, 30.963756532073521417},
990 {5, 4, 38.659808254090090604},
991 {5, 6, 50.194428907734805993},
992 {6, 1, 9.4623222080256173906},
993 {6, 5, 39.805571092265194006}
994 };
995
996 struct angle_table arrow_angles[13] =
997 {{0, 1, 90.0},
998 {1, 0, 0.0},
999 {1, 1, 45.0},
1000 {1, 2, 63.434948822922010648},
1001 {1, 3, 71.565051177077989351},
1002 {1, 4, 75.963756532073521417},
1003 {2, 1, 26.565051177077989351},
1004 {2, 3, 56.309932474020213086},
1005 {3, 1, 18.434948822922010648},
1006 {3, 2, 33.690067525979786913},
1007 {3, 4, 53.130102354155978703},
1008 {4, 1, 14.036243467926478588},
1009 {4, 3, 36.869897645844021297},
1010 };
1011
get_slope(int dx,int dy,int * sxp,int * syp,int arrow)1012 void get_slope(int dx, int dy, int *sxp, int *syp, int arrow)
1013 {
1014 double angle;
1015 int i, s, max;
1016 double d, d1;
1017 struct angle_table *st;
1018
1019 if (dx == 0) {
1020 *sxp = 0;
1021 *syp = signof(dy);
1022 return;
1023 }
1024 angle = atan((double) abs(dy) / (double) abs(dx)) * rad2deg;
1025 if (arrow) {
1026 st = arrow_angles;
1027 max = 13;
1028 } else {
1029 st = line_angles;
1030 max = 25;
1031 }
1032 s = 0;
1033 d = 9.9e9;
1034 for (i = 0; i < max; i++) {
1035 d1 = fabs(angle - st[i].angle);
1036 if (d1 < d) {
1037 s = i;
1038 d = d1;
1039 }
1040 }
1041 *sxp = st[s].x;
1042 if (dx < 0)
1043 *sxp = -*sxp;
1044 *syp = st[s].y;
1045 if (dy < 0)
1046 *syp = -*syp;
1047 }
1048
latex_endpoint(int x1,int y1,int x2,int y2,int * xout,int * yout,int arrow,int magnet)1049 void latex_endpoint(int x1, int y1, int x2, int y2, int *xout, int *yout, int arrow, int magnet)
1050 {
1051 int dx, dy, sx, sy, ds, dsx, dsy;
1052
1053 dx = x2 - x1;
1054 dy = y2 - y1;
1055 get_slope(dx, dy, &sx, &sy, arrow);
1056 if (abs(sx) >= abs(sy)) {
1057 ds = lcm(sx, magnet * gcd(sx, magnet));
1058 dsx = (2 * abs(dx) / ds + 1) / 2;
1059 dsx = (dx >= 0) ? dsx * ds : -dsx * ds;
1060 *xout = x1 + dsx;
1061 *yout = y1 + dsx * sy / sx;
1062 } else {
1063 ds = lcm(sy, magnet * gcd(sy, magnet));
1064 dsy = (2 * abs(dy) / ds + 1) / 2;
1065 dsy = (dy >= 0) ? dsy * ds : -dsy * ds;
1066 *yout = y1 + dsy;
1067 *xout = x1 + dsy * sx / sy;
1068 }
1069 }
1070