1 /* Dia -- an diagram creation/manipulation program
2  * Copyright (C) 1998 Alexander Larsson
3  *
4  * Matrix code from GIMP:app/transform_tool.c
5  * Copyright (C) 1995 Spencer Kimball and Peter Mattis
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20  */
21 
22 /* include glib.h first, so we don't pick up its inline functions as well */
23 #include <config.h>
24 
25 #include <glib.h>
26 /* include normal versions of the inlined functions here ... */
27 #undef G_INLINE_FUNC
28 #define G_INLINE_FUNC extern
29 #define G_CAN_INLINE 1
30 #include "geometry.h"
31 
32 #include "object.h"
33 #include "units.h"
34 
35 void
rectangle_union(Rectangle * r1,const Rectangle * r2)36 rectangle_union(Rectangle *r1, const Rectangle *r2)
37 {
38   r1->top = MIN( r1->top, r2->top );
39   r1->bottom = MAX( r1->bottom, r2->bottom );
40   r1->left = MIN( r1->left, r2->left );
41   r1->right = MAX( r1->right, r2->right );
42 }
43 
44 void
int_rectangle_union(IntRectangle * r1,const IntRectangle * r2)45 int_rectangle_union(IntRectangle *r1, const IntRectangle *r2)
46 {
47   r1->top = MIN( r1->top, r2->top );
48   r1->bottom = MAX( r1->bottom, r2->bottom );
49   r1->left = MIN( r1->left, r2->left );
50   r1->right = MAX( r1->right, r2->right );
51 }
52 
53 void
rectangle_intersection(Rectangle * r1,const Rectangle * r2)54 rectangle_intersection(Rectangle *r1, const Rectangle *r2)
55 {
56   r1->top = MAX( r1->top, r2->top );
57   r1->bottom = MIN( r1->bottom, r2->bottom );
58   r1->left = MAX( r1->left, r2->left );
59   r1->right = MIN( r1->right, r2->right );
60 
61   /* Is rectangle empty? */
62   if ( (r1->top >= r1->bottom) ||
63        (r1->left >= r1->right) ) {
64     r1->top = r1->bottom = r1->left = r1->right = 0.0;
65   }
66 }
67 
68 int
rectangle_intersects(const Rectangle * r1,const Rectangle * r2)69 rectangle_intersects(const Rectangle *r1, const Rectangle *r2)
70 {
71   if ( (r1->right < r2->left) ||
72        (r1->left > r2->right) ||
73        (r1->top > r2->bottom) ||
74        (r1->bottom < r2->top) )
75     return FALSE;
76 
77   return TRUE;
78 }
79 
80 int
point_in_rectangle(const Rectangle * r,const Point * p)81 point_in_rectangle(const Rectangle* r, const Point *p)
82 {
83   if ( (p->x < r->left) ||
84        (p->x > r->right) ||
85        (p->y > r->bottom) ||
86        (p->y < r->top) )
87     return FALSE;
88 
89   return TRUE;
90 }
91 
92 int
rectangle_in_rectangle(const Rectangle * outer,const Rectangle * inner)93 rectangle_in_rectangle(const Rectangle* outer, const Rectangle *inner)
94 {
95   if ( (inner->left < outer->left) ||
96        (inner->right > outer->right) ||
97        (inner->top < outer->top) ||
98        (inner->bottom > outer->bottom))
99     return FALSE;
100 
101   return TRUE;
102 }
103 
104 void
rectangle_add_point(Rectangle * r,const Point * p)105 rectangle_add_point(Rectangle *r, const Point *p)
106 {
107   if (p->x < r->left)
108     r->left = p->x;
109   else if (p->x > r->right)
110     r->right = p->x;
111 
112   if (p->y < r->top)
113     r->top = p->y;
114   else if (p->y > r->bottom)
115     r->bottom = p->y;
116 }
117 
118 
119 /*
120  * This function estimates the distance from a point to a rectangle.
121  * If the point is in the rectangle, 0.0 is returned. Otherwise the
122  * distance in a manhattan metric from the point to the nearest point
123  * on the rectangle is returned.
124  */
125 real
distance_rectangle_point(const Rectangle * rect,const Point * point)126 distance_rectangle_point(const Rectangle *rect, const Point *point)
127 {
128   real dx = 0.0;
129   real dy = 0.0;
130 
131   if (point->x < rect->left ) {
132     dx = rect->left - point->x;
133   } else if (point->x > rect->right ) {
134     dx = point->x - rect->right;
135   }
136 
137   if (point->y < rect->top ) {
138     dy = rect->top - point->y;
139   } else if (point->y > rect->bottom ) {
140     dy = point->y - rect->bottom;
141   }
142 
143   return dx + dy;
144 }
145 
146 /*
147  * This function estimates the distance from a point to a line segment
148  * specified by two endpoints.
149  * If the point is on the line segment, 0.0 is returned. Otherwise the
150  * distance in the R^2 metric from the point to the nearest point
151  * on the line segment is returned. Does one sqrt per call.
152  * Philosophical bug: line_width is ignored iff point is beyond
153  * end of line segment.
154  */
155 real
distance_line_point(const Point * line_start,const Point * line_end,real line_width,const Point * point)156 distance_line_point(const Point *line_start, const Point *line_end,
157 		    real line_width, const Point *point)
158 {
159   Point v1, v2;
160   real v1_lensq;
161   real projlen;
162   real perp_dist;
163 
164   v1 = *line_end;
165   point_sub(&v1,line_start);
166 
167   v2 = *point;
168   point_sub(&v2, line_start);
169 
170   v1_lensq = point_dot(&v1,&v1);
171 
172   if (v1_lensq<0.000001) {
173     return sqrt(point_dot(&v2,&v2));
174   }
175 
176   projlen = point_dot(&v1,&v2) / v1_lensq;
177 
178   if (projlen<0.0) {
179     return sqrt(point_dot(&v2,&v2));
180   }
181 
182   if (projlen>1.0) {
183     Point v3 = *point;
184     point_sub(&v3,line_end);
185     return sqrt(point_dot(&v3,&v3));
186   }
187 
188   point_scale(&v1, projlen);
189   point_sub(&v1, &v2);
190 
191   perp_dist = sqrt(point_dot(&v1,&v1));
192 
193   perp_dist -= line_width / 2.0;
194   if (perp_dist < 0.0) {
195     perp_dist = 0.0;
196   }
197 
198   return perp_dist;
199 }
200 
201 /* returns 1 if the line crosses the ray from (-Inf, rayend.y) to rayend */
202 static int
line_crosses_ray(const Point * line_start,const Point * line_end,const Point * rayend)203 line_crosses_ray(const Point *line_start,
204                  const Point *line_end, const Point *rayend)
205 {
206   if ((line_start->y <= rayend->y && line_end->y > rayend->y) || /* upward crossing */
207       (line_start->y > rayend->y && line_end->y <= rayend->y)) { /* downward crossing */
208     real vt = (rayend->y - line_start->y) / (line_end->y - line_start->y);
209     if (rayend->x < line_start->x + vt * (line_end->x - line_start->x)) /* intersect */
210       return 1;
211   }
212   return 0;
213 }
214 
215 real
distance_polygon_point(const Point * poly,guint npoints,real line_width,const Point * point)216 distance_polygon_point(const Point *poly, guint npoints, real line_width,
217 		       const Point *point)
218 {
219   guint i, last = npoints - 1;
220   real line_dist = G_MAXFLOAT;
221   guint crossings = 0;
222 
223   /* calculate ray crossings and line distances */
224   for (i = 0; i < npoints; i++) {
225     real dist;
226 
227     crossings += line_crosses_ray(&poly[last], &poly[i], point);
228     dist = distance_line_point(&poly[last], &poly[i], line_width, point);
229     line_dist = MIN(line_dist, dist);
230     last = i;
231   }
232   /* If there is an odd number of ray crossings, we are inside the polygon.
233    * Otherwise, return the minium distance from a line segment */
234   if (crossings % 2 == 1)
235     return 0.0;
236   else
237     return line_dist;
238 }
239 
240 /* number of segments to use in bezier curve approximation */
241 #define NBEZ_SEGS 10
242 
243 /* if cross is not NULL, it will be incremented for each ray crossing */
244 static real
bez_point_distance_and_ray_crosses(const Point * b1,const Point * b2,const Point * b3,const Point * b4,real line_width,const Point * point,guint * cross)245 bez_point_distance_and_ray_crosses(const Point *b1,
246                                    const Point *b2, const Point *b3,
247                                    const Point *b4,
248 				   real line_width, const Point *point,
249                                    guint *cross)
250 {
251   static gboolean calculated_coeff = FALSE;
252   static real coeff[NBEZ_SEGS+1][4];
253   int i;
254   real line_dist = G_MAXFLOAT;
255   Point prev, pt;
256 
257   if (!calculated_coeff) {
258     for (i = 0; i <= NBEZ_SEGS; i++) {
259       real t1 = ((real)i)/NBEZ_SEGS, t2 = t1*t1, t3 = t1*t2;
260       real it1 = 1-t1, it2 = it1*it1, it3 = it1*it2;
261 
262       coeff[i][0] = it3;
263       coeff[i][1] = 3 * t1 * it2;
264       coeff[i][2] = 3 * t2 * it1;
265       coeff[i][3] = t3;
266     }
267   }
268   calculated_coeff = TRUE;
269 
270   prev.x = coeff[0][0] * b1->x + coeff[0][1] * b2->x +
271            coeff[0][2] * b3->x + coeff[0][3] * b4->x;
272   prev.y = coeff[0][0] * b1->y + coeff[0][1] * b2->y +
273            coeff[0][2] * b3->y + coeff[0][3] * b4->y;
274   for (i = 1; i <= NBEZ_SEGS; i++) {
275     real dist;
276 
277     pt.x = coeff[i][0] * b1->x + coeff[i][1] * b2->x +
278            coeff[i][2] * b3->x + coeff[i][3] * b4->x;
279     pt.y = coeff[i][0] * b1->y + coeff[i][1] * b2->y +
280            coeff[i][2] * b3->y + coeff[i][3] * b4->y;
281 
282     dist = distance_line_point(&prev, &pt, line_width, point);
283     line_dist = MIN(line_dist, dist);
284     if (cross)
285       *cross += line_crosses_ray(&prev, &pt, point);
286 
287     prev = pt;
288   }
289   return line_dist;
290 }
291 
292 real
distance_bez_seg_point(const Point * b1,const Point * b2,const Point * b3,const Point * b4,real line_width,const Point * point)293 distance_bez_seg_point(const Point *b1, const Point *b2,
294                        const Point *b3, const Point *b4,
295 		       real line_width, const Point *point)
296 {
297   return bez_point_distance_and_ray_crosses(b1, b2, b3, b4,
298 					    line_width, point, NULL);
299 }
300 
301 real
distance_bez_line_point(const BezPoint * b,guint npoints,real line_width,const Point * point)302 distance_bez_line_point(const BezPoint *b, guint npoints,
303 			real line_width, const Point *point)
304 {
305   Point last;
306   guint i;
307   real line_dist = G_MAXFLOAT;
308 
309   g_return_val_if_fail(b[0].type == BEZ_MOVE_TO, -1);
310 
311   last = b[0].p1;
312 
313   for (i = 1; i < npoints; i++) {
314     real dist;
315 
316     switch (b[i].type) {
317     case BEZ_MOVE_TO:
318       /*g_assert_not_reached();*/
319       g_warning("BEZ_MOVE_TO found half way through a bezier line");
320       break;
321     case BEZ_LINE_TO:
322       dist = distance_line_point(&last, &b[i].p1, line_width, point);
323       line_dist = MIN(line_dist, dist);
324       last = b[i].p1;
325       break;
326     case BEZ_CURVE_TO:
327       dist = bez_point_distance_and_ray_crosses(&last, &b[i].p1, &b[i].p2,
328 						&b[i].p3, line_width, point,
329 						NULL);
330       line_dist = MIN(line_dist, dist);
331       last = b[i].p3;
332       break;
333     }
334   }
335   return line_dist;
336 }
337 
338 real
distance_bez_shape_point(const BezPoint * b,guint npoints,real line_width,const Point * point)339 distance_bez_shape_point(const BezPoint *b, guint npoints,
340                          real line_width, const Point *point)
341 {
342   Point last;
343   guint i;
344   real line_dist = G_MAXFLOAT;
345   guint crossings = 0;
346 
347   g_return_val_if_fail(b[0].type == BEZ_MOVE_TO, -1);
348 
349   last = b[0].p1;
350 
351   for (i = 1; i < npoints; i++) {
352     real dist;
353 
354     switch (b[i].type) {
355     case BEZ_MOVE_TO:
356       /*g_assert_not_reached();*/
357       g_warning("BEZ_MOVE_TO found half way through a bezier shape");
358       break;
359     case BEZ_LINE_TO:
360       dist = distance_line_point(&last, &b[i].p1, line_width, point);
361       crossings += line_crosses_ray(&last, &b[i].p1, point);
362       line_dist = MIN(line_dist, dist);
363       last = b[i].p1;
364       break;
365     case BEZ_CURVE_TO:
366       dist = bez_point_distance_and_ray_crosses(&last, &b[i].p1, &b[i].p2,
367 						&b[i].p3, line_width, point,
368 						&crossings);
369       line_dist = MIN(line_dist, dist);
370       last = b[i].p3;
371       break;
372     }
373   }
374   /* If there is an odd number of ray crossings, we are inside the polygon.
375    * Otherwise, return the minium distance from a line segment */
376   if (crossings % 2 == 1)
377     return 0.0;
378   else
379     return line_dist;
380 }
381 
382 real
distance_ellipse_point(const Point * centre,real width,real height,real line_width,const Point * point)383 distance_ellipse_point(const Point *centre, real width, real height,
384 		       real line_width, const Point *point)
385 {
386   /* A faster intersection method would be to scaled the ellipse and the
387    * point to where the ellipse is a circle, intersect with the circle,
388    * then scale back.
389    */
390   real w2 = width*width, h2 = height*height;
391   real scale, rad, dist;
392   Point pt;
393 
394   /* find the radius of the ellipse in the appropriate direction */
395 
396   /* find the point of intersection between line (x=cx+(px-cx)t; y=cy+(py-cy)t)
397    * and ellipse ((x-cx)^2)/(w/2)^2 + ((y-cy)^2)/(h/2)^2 = 1 */
398   /* radius along ray is sqrt((px-cx)^2 * t^2 + (py-cy)^2 * t^2) */
399 
400   pt = *point;
401   point_sub(&pt, centre);
402   pt.x *= pt.x;
403   pt.y *= pt.y;  /* pt = (point - centre).^2 */
404 
405   scale = w2 * h2 / (4*h2*pt.x + 4*w2*pt.y);
406   rad = sqrt((pt.x + pt.y)*scale) + line_width/2;
407 
408   dist = sqrt(pt.x + pt.y);
409 
410   if (dist <= rad)
411     return 0.0;
412   else
413     return dist - rad;
414 }
415 
416 
417 void
transform_point(Matrix m,Point * src,Point * dest)418 transform_point (Matrix m, Point *src, Point *dest)
419 {
420   real xx, yy, ww;
421 
422   xx = m[0][0] * src->x + m[0][1] * src->y + m[0][2];
423   yy = m[1][0] * src->x + m[1][1] * src->y + m[1][2];
424   ww = m[2][0] * src->x + m[2][1] * src->y + m[2][2];
425 
426   if (!ww)
427     ww = 1.0;
428 
429   dest->x = xx / ww;
430   dest->y = yy / ww;
431 }
432 
433 void
mult_matrix(Matrix m1,Matrix m2)434 mult_matrix (Matrix m1, Matrix m2)
435 {
436   Matrix result;
437   int i, j, k;
438 
439   for (i = 0; i < 3; i++)
440     for (j = 0; j < 3; j++)
441       {
442 	result [i][j] = 0.0;
443 	for (k = 0; k < 3; k++)
444 	  result [i][j] += m1 [i][k] * m2[k][j];
445       }
446 
447   /*  copy the result into matrix 2  */
448   for (i = 0; i < 3; i++)
449     for (j = 0; j < 3; j++)
450       m2 [i][j] = result [i][j];
451 }
452 
453 void
identity_matrix(Matrix m)454 identity_matrix (Matrix m)
455 {
456   int i, j;
457 
458   for (i = 0; i < 3; i++)
459     for (j = 0; j < 3; j++)
460       m[i][j] = (i == j) ? 1 : 0;
461 
462 }
463 
464 void
translate_matrix(Matrix m,real x,real y)465 translate_matrix (Matrix m, real x, real y)
466 {
467   Matrix trans;
468 
469   identity_matrix (trans);
470   trans[0][2] = x;
471   trans[1][2] = y;
472   mult_matrix (trans, m);
473 }
474 
475 void
scale_matrix(Matrix m,real x,real y)476 scale_matrix (Matrix m, real x, real y)
477 {
478   Matrix scale;
479 
480   identity_matrix (scale);
481   scale[0][0] = x;
482   scale[1][1] = y;
483   mult_matrix (scale, m);
484 }
485 
486 void
rotate_matrix(Matrix m,real theta)487 rotate_matrix (Matrix m, real theta)
488 {
489   Matrix rotate;
490   real cos_theta, sin_theta;
491 
492   cos_theta = cos (theta);
493   sin_theta = sin (theta);
494 
495   identity_matrix (rotate);
496   rotate[0][0] = cos_theta;
497   rotate[0][1] = -sin_theta;
498   rotate[1][0] = sin_theta;
499   rotate[1][1] = cos_theta;
500   mult_matrix (rotate, m);
501 }
502 
503 void
xshear_matrix(Matrix m,real shear)504 xshear_matrix (Matrix m, real shear)
505 {
506   Matrix shear_m;
507 
508   identity_matrix (shear_m);
509   shear_m[0][1] = shear;
510   mult_matrix (shear_m, m);
511 }
512 
513 void
yshear_matrix(Matrix m,real shear)514 yshear_matrix (Matrix m, real shear)
515 {
516   Matrix shear_m;
517 
518   identity_matrix (shear_m);
519   shear_m[1][0] = shear;
520   mult_matrix (shear_m, m);
521 }
522 
523 /*  find the determinate for a 3x3 matrix  */
524 static real
determinate(Matrix m)525 determinate (Matrix m)
526 {
527   int i;
528   double det = 0;
529 
530   for (i = 0; i < 3; i ++)
531     {
532       det += m[0][i] * m[1][(i+1)%3] * m[2][(i+2)%3];
533       det -= m[2][i] * m[1][(i+1)%3] * m[0][(i+2)%3];
534     }
535 
536   return det;
537 }
538 
539 
540 
541 void
point_convex(Point * dst,const Point * src1,const Point * src2,real alpha)542 point_convex(Point *dst, const Point *src1, const Point *src2, real alpha)
543 {
544   /* Make convex combination of src1 and src2:
545      dst = alpha * src1 + (1-alpha) * src2;
546   */
547   point_copy(dst,src1);
548   point_scale(dst,alpha);
549   point_add_scaled(dst,src2,1.0 - alpha);
550 }
551 
552 
553 
554 /* The following functions were originally taken from Graphics Gems III,
555    pg 193 and corresponding code pgs 496-499
556 
557    They were modified to work within the dia coordinate system
558  */
559 /* return angle subtended by two vectors */
dot2(Point * p1,Point * p2)560 real dot2(Point *p1, Point *p2)
561 {
562   real d, t;
563   d = sqrt(((p1->x*p1->x)+(p1->y*p1->y)) *
564            ((p2->x*p2->x)+(p2->y*p2->y)));
565   if ( d != 0.0 )
566   {
567     t = (p1->x*p2->x+p1->y*p2->y)/d;
568     return (dia_acos(t));
569   }
570   else
571   {
572     return 0.0;
573   }
574 }
575 
576 /* Find a,b,c in ax + by + c = 0 for line p1, p2 */
line_coef(real * a,real * b,real * c,Point * p1,Point * p2)577 void line_coef(real *a, real *b, real *c, Point *p1, Point *p2)
578 {
579   *c = ( p2->x * p1->y ) - ( p1->x * p2->y );
580   *a = p2->y - p1->y;
581   *b = p1->x - p2->x;
582 }
583 
584 /* Return signed distance from line
585    ax + by + c = 0 to point p.              */
line_to_point(real a,real b,real c,Point * p)586 real line_to_point(real a, real b , real c, Point *p) {
587   real d, lp;
588   d = sqrt((a*a)+(b*b));
589   if ( d == 0.0 ) {
590     lp = 0.0;
591   } else {
592     lp = (a*p->x+b*p->y+c)/d;
593   }
594   return (lp);
595 }
596 
597 /* Given line l = ax + by + c = 0 and point p,
598    compute x,y so point perp is perpendicular to line l. */
point_perp(Point * p,real a,real b,real c,Point * perp)599 void point_perp(Point *p, real a, real b, real c, Point *perp) {
600   real d, cp;
601   perp->x = 0.0;
602   perp->y = 0.0;
603   d = a*a + b*b;
604   cp = a*p->y-b*p->x;
605   if ( d != 0.0 ) {
606     perp->x = (-a*c-b*cp)/d;
607     perp->y = (a*cp-b*c)/d;
608   }
609   return;
610 }
611 
612 /* Compute a circular arc fillet between lines L1 (p1 to p2)
613    and L2 (p3 to p4) with radius r.
614    The circle center is c.
615    The start angle is pa
616    The end angle is aa,
617    The points p1-p4 will be modified as necessary */
fillet(Point * p1,Point * p2,Point * p3,Point * p4,real r,Point * c,real * pa,real * aa)618 void fillet(Point *p1, Point *p2, Point *p3, Point *p4,
619                    real r, Point *c, real *pa, real *aa) {
620   real a1, b1, c1;  /* Coefficients for L1 */
621   real a2, b2, c2;  /* Coefficients for L2 */
622   real d, d1, d2;
623   real c1p, c2p;
624   real rr;
625   real start_angle, stop_angle;
626   Point mp,gv1,gv2;
627   gboolean righthand = FALSE;
628 
629   line_coef(&a1,&b1,&c1,p1,p2);
630   line_coef(&a2,&b2,&c2,p3,p4);
631 
632   if ( (a1*b2) == (a2*b1) ) /* Parallel or coincident lines */
633   {
634     return;
635   }
636 
637   mp.x = (p3->x + p4->x) / 2.0;          /* Find midpoint of p3 p4 */
638   mp.y = (p3->y + p4->y) / 2.0;
639   d1 = line_to_point(a1, b1, c1, &mp);    /* Find distance p1 p2 to
640                                              midpoint p3 p4 */
641   if ( d1 == 0.0 ) return;                /* p1p2 to p3 */
642 
643   mp.x = (p1->x + p2->x) / 2.0;          /* Find midpoint of p1 p2 */
644   mp.y = (p1->y + p2->y) / 2.0;
645   d2 = line_to_point(a2, b2, c2, &mp);    /* Find distance p3 p4 to
646                                              midpoint p1 p2 */
647   if ( d2 == 0.0 ) return;
648 
649   rr = r;
650   if ( d1 <= 0.0 ) rr = -rr;
651 
652   c1p = c1-rr*sqrt((a1*a1)+(b1*b1));     /* Line parallell l1 at d */
653   rr = r;
654 
655   if ( d2 <= 0.0 ) rr = -rr;
656   c2p = c2-rr*sqrt((a2*a2)+(b2*b2));     /* Line parallell l2 at d */
657   d = a1*b2-a2*b1;
658 
659   c->x = ( c2p*b1-c1p*b2 ) / d;          /* Intersect constructed lines */
660   c->y = ( c1p*a2-c2p*a1 ) / d;          /* to find center of arc */
661 
662   point_perp(c,a1,b1,c1,p2);             /* Clip or extend lines as required */
663   point_perp(c,a2,b2,c2,p3);
664 
665   /* need to negate the y coords to calculate angles correctly */
666   gv1.x = p2->x-c->x; gv1.y = -(p2->y-c->y);
667   gv2.x = p3->x-c->x; gv2.y = -(p3->y-c->y);
668 
669   start_angle = atan2(gv1.y,gv1.x);   /* Beginning angle for arc */
670   stop_angle = dot2(&gv1,&gv2);
671   if ( point_cross(&gv1,&gv2) < 0.0 ) {
672     stop_angle = -stop_angle; /* Angle subtended by arc */
673     /* also this means we need to swap our angles later on */
674     righthand = TRUE;
675   }
676   /* now calculate the actual angles in a form that the draw arc function
677      of the renderer can use */
678   start_angle = start_angle*180.0/G_PI;
679   stop_angle  = start_angle + stop_angle*180.0/G_PI;
680   while (start_angle < 0.0) start_angle += 360.0;
681   while (stop_angle < 0.0)  stop_angle += 360.0;
682   /* swap the start and stop if we had to negate the cross product */
683   if ( righthand ) {
684     real tmp = start_angle;
685     start_angle = stop_angle;
686     stop_angle = tmp;
687   }
688   *pa = start_angle;
689   *aa = stop_angle;
690 }
691 
692 int
three_point_circle(const Point * p1,const Point * p2,const Point * p3,Point * center,real * radius)693 three_point_circle (const Point *p1, const Point *p2, const Point *p3,
694                     Point* center, real* radius)
695 {
696   const real epsilon = 1e-4;
697   real x1 = p1->x;
698   real y1 = p1->y;
699   real x2 = p2->x;
700   real y2 = p2->y;
701   real x3 = p3->x;
702   real y3 = p3->y;
703   real ma, mb;
704   if (fabs(x2 - x1) < epsilon)
705     return 0;
706   if (fabs(x3 - x2) < epsilon)
707     return 0;
708   ma = (y2 - y1) / (x2 - x1);
709   mb = (y3 - y2) / (x3 - x2);
710   if (fabs (mb - ma) < epsilon)
711     return 0;
712   center->x = (ma*mb*(y1-y3)+mb*(x1+x2)-ma*(x2+x3))/(2*(mb-ma));
713   if (fabs(ma)>epsilon)
714     center->y = -1/ma*(center->x - (x1+x2)/2.0) + (y1+y2)/2.0;
715   else if (fabs(mb)>epsilon)
716     center->y = -1/mb*(center->x - (x2+x3)/2.0) + (y2+y3)/2.0;
717   else
718     return 0;
719   *radius = distance_point_point(center, p1);
720   return 1;
721 }
722 
723 
724 /* moved this here since it is being reused by rounded polyline code*/
725 real
point_cross(Point * p1,Point * p2)726 point_cross(Point *p1, Point *p2)
727 {
728   return p1->x * p2->y - p2->x * p1->y;
729 }
730 
731 /* moved this here since it is used by line and bezier */
732 /* Computes one point between end and objmid which
733  * touches the edge of the object */
734 Point
calculate_object_edge(Point * objmid,Point * end,DiaObject * obj)735 calculate_object_edge(Point *objmid, Point *end, DiaObject *obj)
736 {
737 #define MAXITER 25
738 #ifdef TRACE_DIST
739   Point trace[MAXITER];
740   real disttrace[MAXITER];
741 #endif
742   Point mid1, mid2, mid3;
743   real dist;
744   int i = 0;
745 
746   mid1 = *objmid;
747   mid2.x = (objmid->x+end->x)/2;
748   mid2.y = (objmid->y+end->y)/2;
749   mid3 = *end;
750 
751   /* If the other end is inside the object */
752   dist = obj->ops->distance_from(obj, &mid3);
753   if (dist < 0.001) return mid1;
754 
755 
756   do {
757     dist = obj->ops->distance_from(obj, &mid2);
758     if (dist < 0.0000001) {
759       mid1 = mid2;
760     } else {
761       mid3 = mid2;
762     }
763     mid2.x = (mid1.x + mid3.x)/2;
764     mid2.y = (mid1.y + mid3.y)/2;
765 #ifdef TRACE_DIST
766     trace[i] = mid2;
767     disttrace[i] = dist;
768 #endif
769     i++;
770   } while (i < MAXITER && (dist < 0.0000001 || dist > 0.001));
771 
772 #ifdef TRACE_DIST
773   if (i == MAXITER) {
774     for (i = 0; i < MAXITER; i++) {
775       printf("%d: %f, %f: %f\n", i, trace[i].x, trace[i].y, disttrace[i]);
776     }
777     printf("i = %d, dist = %f\n", i, dist);
778   }
779 #endif
780 
781   return mid2;
782 }
783 
784 /*!
785  * \brief asin wrapped to limit to valid result range
786  *
787  * Although clipping the value might hide some miscalculation
788  * elsewhere this should still be used in all of Dia. Continuing
789  * calculation with bogus values might final fall on our foots
790  * because rendering libraries might not be graceful.
791  * See https://bugzilla.gnome.org/show_bug.cgi?id=710818
792  */
793 real
dia_asin(real x)794 dia_asin (real x)
795 {
796   /* clamp to valid range */
797   if (x <= -1.0)
798     return -M_PI/2;
799   if (x >= 1.0)
800     return M_PI/2;
801   return asin (x);
802 }
803 
804 /*!
805  * \brief acos wrapped to limit to valid result range
806  */
807 real
dia_acos(real x)808 dia_acos (real x)
809 {
810   /* clamp to valid range */
811   if (x <= -1.0)
812     return M_PI;
813   if (x >= 1.0)
814     return 0.0;
815   return acos (x);
816 }
817