Kstars

overlap.h
1 /* Licensed under a 3-clause BSD style license.
2 
3  Functions for calculating exact overlap between shapes.
4 
5  Original cython version by Thomas Robitaille. Converted to C by Kyle
6  Barbary.*/
7 
8 #pragma once
9 
10 #include <math.h>
11 
12 /* Return area of a circle arc between (x1, y1) and (x2, y2) with radius r */
13 /* reference: http://mathworld.wolfram.com/CircularSegment.html */
14 static double area_arc(double x1, double y1, double x2, double y2,
15  double r)
16 {
17  double a, theta;
18 
19  a = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
20  theta = 2. * asin(0.5 * a / r);
21  return 0.5 * r * r * (theta - sin(theta));
22 }
23 
24 /* Area of a triangle defined by three vertices */
25 static double area_triangle(double x1, double y1, double x2, double y2,
26  double x3, double y3)
27 {
28  return 0.5 * fabs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2));
29 }
30 
31 /* Core of circular overlap routine.
32  * Assumes that xmax >= xmin >= 0.0, ymax >= ymin >= 0.0.
33  * (can always modify input to conform to this).
34  */
35 static double circoverlap_core(double xmin, double ymin,
36  double xmax, double ymax, double r)
37 {
38  double a, b, x1, x2, y1, y2, r2, xmin2, ymin2, xmax2, ymax2;
39 
40  xmin2 = xmin*xmin;
41  ymin2 = ymin*ymin;
42  r2 = r*r;
43  if (xmin2 + ymin2 > r2)
44  return 0.;
45 
46  xmax2 = xmax*xmax;
47  ymax2 = ymax*ymax;
48  if (xmax2 + ymax2 < r2)
49  return (xmax-xmin)*(ymax-ymin);
50 
51  a = xmax2 + ymin2; /* (corner 1 distance)^2 */
52  b = xmin2 + ymax2; /* (corner 2 distance)^2 */
53 
54  if (a < r2 && b < r2)
55  {
56  x1 = sqrt(r2 - ymax2);
57  y1 = ymax;
58  x2 = xmax;
59  y2 = sqrt(r2 - xmax2);
60  return ((xmax-xmin)*(ymax-ymin) -
61  area_triangle(x1, y1, x2, y2, xmax, ymax) +
62  area_arc(x1, y1, x2, y2, r));
63  }
64 
65  if (a < r2)
66  {
67  x1 = xmin;
68  y1 = sqrt(r2 - xmin2);
69  x2 = xmax;
70  y2 = sqrt(r2 - xmax2);
71  return (area_arc(x1, y1, x2, y2, r) +
72  area_triangle(x1, y1, x1, ymin, xmax, ymin) +
73  area_triangle(x1, y1, x2, ymin, x2, y2));
74  }
75 
76  if (b < r2)
77  {
78  x1 = sqrt(r2 - ymin2);
79  y1 = ymin;
80  x2 = sqrt(r2 - ymax2);
81  y2 = ymax;
82  return (area_arc(x1, y1, x2, y2, r) +
83  area_triangle(x1, y1, xmin, y1, xmin, ymax) +
84  area_triangle(x1, y1, xmin, y2, x2, y2));
85  }
86 
87  /* else */
88  x1 = sqrt(r2 - ymin2);
89  y1 = ymin;
90  x2 = xmin;
91  y2 = sqrt(r2 - xmin2);
92  return (area_arc(x1, y1, x2, y2, r) +
93  area_triangle(x1, y1, x2, y2, xmin, ymin));
94 }
95 
96 
97 
98 /* Area of overlap of a rectangle and a circle */
99 static double circoverlap(double xmin, double ymin, double xmax, double ymax,
100  double r)
101 {
102  /* some subroutines demand that r > 0 */
103  if (r <= 0.)
104  return 0.;
105 
106  if (0. <= xmin)
107  {
108  if (0. <= ymin)
109  return circoverlap_core(xmin, ymin, xmax, ymax, r);
110  else if (0. >= ymax)
111  return circoverlap_core(-ymax, xmin, -ymin, xmax, r);
112  else
113  return (circoverlap(xmin, ymin, xmax, 0., r) +
114  circoverlap(xmin, 0., xmax, ymax, r));
115  }
116  else if (0. >= xmax)
117  {
118  if (0. <= ymin)
119  return circoverlap_core(-xmax, ymin, -xmin, ymax, r);
120  else if (0. >= ymax)
121  return circoverlap_core(-xmax, -ymax, -xmin, -ymin, r);
122  else
123  return (circoverlap(xmin, ymin, xmax, 0., r) +
124  circoverlap(xmin, 0., xmax, ymax, r));
125  }
126  else
127  {
128  if (0. <= ymin)
129  return (circoverlap(xmin, ymin, 0., ymax, r) +
130  circoverlap(0., ymin, xmax, ymax, r));
131  if (0. >= ymax)
132  return (circoverlap(xmin, ymin, 0., ymax, r) +
133  circoverlap(0., ymin, xmax, ymax, r));
134  else
135  return (circoverlap(xmin, ymin, 0., 0., r) +
136  circoverlap(0., ymin, xmax, 0., r) +
137  circoverlap(xmin, 0., 0., ymax, r) +
138  circoverlap(0., 0., xmax, ymax, r));
139  }
140 }
141 
142 /*
143 Start of new circular overlap routine that might be faster.
144 
145 double circoverlap_new(double dx, double dy, double r)
146 {
147  double xmin, xmax, ymin, ymax, xmin2, xmax2, ymin2, ymax2, r2;
148 
149  if (dx < 0.)
150  dx = -dx;
151  if (dy < 0.)
152  dy = -dy;
153  if (dy > dx)
154  {
155  r2 = dy;
156  dy = dx;
157  dx = r2;
158  }
159 
160  xmax = dx + 0.5;
161  ymax = dy + 0.5;
162  xmax2 = xmax*xmax;
163  ymax2 = ymax*ymax;
164  r2 = r*r;
165 
166  if (xmax2 + ymax2 < r2)
167  return 1.;
168 
169  xmin2 = xmin*xmin;
170  if (xmin2 +
171 }
172 */
173 
174 /*****************************************************************************/
175 /* ellipse overlap functions */
176 
177 typedef struct
178 {
179  double x, y;
180 } point;
181 
182 
183 typedef struct
184 {
185  point p1, p2;
186 } intersections;
187 
188 static void swap(double *a, double *b)
189 {
190  double temp;
191  temp = *a;
192  *a = *b;
193  *b = temp;
194 }
195 
196 static void swap_point(point *a, point *b)
197 {
198  point temp;
199  temp = *a;
200  *a = *b;
201  *b = temp;
202 }
203 
204 /* rotate values to the left: a=b, b=c, c=a */
205 static void rotate(double *a, double *b, double *c)
206 {
207  double temp;
208  temp = *a;
209  *a = *b;
210  *b = *c;
211  *c = temp;
212 }
213 
214 /* Check if a point (x,y) is inside a triangle */
215 static int in_triangle(double x, double y, double x1, double y1,
216  double x2, double y2, double x3, double y3)
217 {
218  int c;
219 
220  c = 0;
221  c += ((y1 > y) != (y2 > y)) && (x < (x2-x1) * (y-y1) / (y2-y1) + x1);
222  c += ((y2 > y) != (y3 > y)) && (x < (x3-x2) * (y-y2) / (y3-y2) + x2);
223  c += ((y3 > y) != (y1 > y)) && (x < (x1-x3) * (y-y3) / (y1-y3) + x3);
224 
225  return c % 2 == 1;
226 }
227 
228 
229 /* Intersection of a line defined by two points with a unit circle */
230 static intersections circle_line(double x1, double y1, double x2, double y2)
231 {
232  double a, b, delta, dx, dy;
233  double tol = 1.e-10;
234  intersections inter;
235 
236  dx = x2 - x1;
237  dy = y2 - y1;
238 
239  if (fabs(dx) < tol && fabs(dy) < tol)
240  {
241  inter.p1.x = 2.;
242  inter.p1.y = 2.;
243  inter.p2.x = 2.;
244  inter.p2.y = 2.;
245  }
246  else if (fabs(dx) > fabs(dy))
247  {
248  /* Find the slope and intercept of the line */
249  a = dy / dx;
250  b = y1 - a * x1;
251 
252  /* Find the determinant of the quadratic equation */
253  delta = 1. + a*a - b*b;
254  if (delta > 0.) /* solutions exist */
255  {
256  delta = sqrt(delta);
257 
258  inter.p1.x = (-a*b - delta) / (1. + a*a);
259  inter.p1.y = a * inter.p1.x + b;
260  inter.p2.x = (-a*b + delta) / (1. + a*a);
261  inter.p2.y = a * inter.p2.x + b;
262  }
263  else /* no solution, return values > 1 */
264  {
265  inter.p1.x = 2.;
266  inter.p1.y = 2.;
267  inter.p2.x = 2.;
268  inter.p2.y = 2.;
269  }
270  }
271  else
272  {
273  /* Find the slope and intercept of the line */
274  a = dx / dy;
275  b = x1 - a * y1;
276 
277  /* Find the determinant of the quadratic equation */
278  delta = 1. + a*a - b*b;
279  if (delta > 0.) /* solutions exist */
280  {
281  delta = sqrt(delta);
282  inter.p1.y = (-a*b - delta) / (1. + a*a);
283  inter.p1.x = a * inter.p1.y + b;
284  inter.p2.y = (-a*b + delta) / (1. + a*a);
285  inter.p2.x = a * inter.p2.y + b;
286  }
287  else /* no solution, return values > 1 */
288  {
289  inter.p1.x = 2.;
290  inter.p1.y = 2.;
291  inter.p2.x = 2.;
292  inter.p2.y = 2.;
293  }
294  }
295 
296  return inter;
297 }
298 
299 /* The intersection of a line with the unit circle. The intersection the
300  * closest to (x2, y2) is chosen. */
301 static point circle_segment_single2(double x1, double y1, double x2, double y2)
302 {
303  double dx1, dy1, dx2, dy2;
304  intersections inter;
305  point pt1, pt2, pt;
306 
307  inter = circle_line(x1, y1, x2, y2);
308  pt1 = inter.p1;
309  pt2 = inter.p2;
310 
311  /*Can be optimized, but just checking for correctness right now */
312  dx1 = fabs(pt1.x - x2);
313  dy1 = fabs(pt1.y - y2);
314  dx2 = fabs(pt2.x - x2);
315  dy2 = fabs(pt2.y - y2);
316 
317  if (dx1 > dy1) /* compare based on x-axis */
318  pt = (dx1 > dx2) ? pt2 : pt1;
319  else
320  pt = (dy1 > dy2) ? pt2 : pt1;
321 
322  return pt;
323 }
324 
325 /* Intersection(s) of a segment with the unit circle. Discard any
326  solution not on the segment. */
327 intersections circle_segment(double x1, double y1, double x2, double y2)
328 {
329  intersections inter, inter_new;
330  point pt1, pt2;
331 
332  inter = circle_line(x1, y1, x2, y2);
333  pt1 = inter.p1;
334  pt2 = inter.p2;
335 
336  if ((pt1.x > x1 && pt1.x > x2) || (pt1.x < x1 && pt1.x < x2) ||
337  (pt1.y > y1 && pt1.y > y2) || (pt1.y < y1 && pt1.y < y2))
338  {
339  pt1.x = 2.;
340  pt1.y = 2.;
341  }
342  if ((pt2.x > x1 && pt2.x > x2) || (pt2.x < x1 && pt2.x < x2) ||
343  (pt2.y > y1 && pt2.y > y2) || (pt2.y < y1 && pt2.y < y2))
344  {
345  pt2.x = 2.;
346  pt2.y = 2.;
347  }
348 
349  if (pt1.x > 1. && pt2.x < 2.)
350  {
351  inter_new.p1 = pt1;
352  inter_new.p2 = pt2;
353  }
354  else
355  {
356  inter_new.p1 = pt2;
357  inter_new.p2 = pt1;
358  }
359  return inter_new;
360 }
361 
362 
363 /* Given a triangle defined by three points (x1, y1), (x2, y2), and
364  (x3, y3), find the area of overlap with the unit circle. */
365 static double triangle_unitcircle_overlap(double x1, double y1,
366  double x2, double y2,
367  double x3, double y3)
368 {
369  double d1, d2, d3, area, xp, yp;
370  int in1, in2, in3, on1, on2, on3;
371  int intersect13, intersect23;
372  intersections inter;
373  point pt1, pt2, pt3, pt4, pt5, pt6;
374 
375  /* Find distance of all vertices to circle center */
376  d1 = x1*x1 + y1*y1;
377  d2 = x2*x2 + y2*y2;
378  d3 = x3*x3 + y3*y3;
379 
380  /* Order vertices by distance from origin */
381  if (d1 < d2)
382  {
383  if (d2 < d3)
384  {}
385  else if (d1 < d3)
386  {
387  swap(&x2, &x3);
388  swap(&y2, &y3);
389  swap(&d2, &d3);
390  }
391  else
392  {
393  rotate(&x1, &x3, &x2);
394  rotate(&y1, &y3, &y2);
395  rotate(&d1, &d3, &d2);
396  }
397  }
398  else
399  {
400  if (d1 < d3)
401  {
402  swap(&x1, &x2);
403  swap(&y1, &y2);
404  swap(&d1, &d2);
405  }
406  else if (d2 < d3)
407  {
408  rotate(&x1, &x2, &x3);
409  rotate(&y1, &y2, &y3);
410  rotate(&d1, &d2, &d3);
411  }
412  else
413  {
414  swap(&x1, &x3);
415  swap(&y1, &y3);
416  swap(&d1, &d3);
417  }
418  }
419 
420  /* Determine number of vertices inside circle */
421  in1 = d1 < 1.;
422  in2 = d2 < 1.;
423  in3 = d3 < 1.;
424 
425  /* Determine which vertices are on the circle */
426  on1 = fabs(d1 - 1.) < 1.e-10;
427  on2 = fabs(d2 - 1.) < 1.e-10;
428  on3 = fabs(d3 - 1.) < 1.e-10;
429 
430  if (on3 || in3) /* triangle completely within circle */
431  area = area_triangle(x1, y1, x2, y2, x3, y3);
432 
433  else if (in2 || on2)
434  {
435  /* If vertex 1 or 2 are on the edge of the circle, then we use
436  * the dot product to vertex 3 to determine whether an
437  * intersection takes place. */
438  intersect13 = !on1 || (x1*(x3-x1) + y1*(y3-y1) < 0.);
439  intersect23 = !on2 || (x2*(x3-x2) + y2*(y3-y2) < 0.);
440  if (intersect13 && intersect23)
441  {
442  pt1 = circle_segment_single2(x1, y1, x3, y3);
443  pt2 = circle_segment_single2(x2, y2, x3, y3);
444  area = (area_triangle(x1, y1, x2, y2, pt1.x, pt1.y) +
445  area_triangle(x2, y2, pt1.x, pt1.y, pt2.x, pt2.y) +
446  area_arc(pt1.x, pt1.y, pt2.x, pt2.y, 1.));
447  }
448  else if (intersect13)
449  {
450  pt1 = circle_segment_single2(x1, y1, x3, y3);
451  area = (area_triangle(x1, y1, x2, y2, pt1.x, pt1.y) +
452  area_arc(x2, y2, pt1.x, pt1.y, 1.));
453  }
454  else if (intersect23)
455  {
456  pt2 = circle_segment_single2(x2, y2, x3, y3);
457  area = (area_triangle(x1, y1, x2, y2, pt2.x, pt2.y) +
458  area_arc(x1, y1, pt2.x, pt2.y, 1.));
459  }
460  else
461  area = area_arc(x1, y1, x2, y2, 1.);
462  }
463  else if (in1)
464  {
465  /* Check for intersections of far side with circle */
466  inter = circle_segment(x2, y2, x3, y3);
467  pt1 = inter.p1;
468  pt2 = inter.p2;
469  pt3 = circle_segment_single2(x1, y1, x2, y2);
470  pt4 = circle_segment_single2(x1, y1, x3, y3);
471 
472  if (pt1.x > 1.) /* indicates no intersection */
473  {
474  /* check if the pixel vertex (x1, y2) and the origin are on
475  * different sides of the circle segment. If they are, the
476  * circle segment spans more than pi radians.
477  * We use the formula (y-y1) * (x2-x1) > (y2-y1) * (x-x1)
478  * to determine if (x, y) is on the left of the directed
479  * line segment from (x1, y1) to (x2, y2) */
480  if (((0.-pt3.y) * (pt4.x-pt3.x) > (pt4.y-pt3.y) * (0.-pt3.x)) !=
481  ((y1-pt3.y) * (pt4.x-pt3.x) > (pt4.y-pt3.y) * (x1-pt3.x)))
482  {
483  area = (area_triangle(x1, y1, pt3.x, pt3.y, pt4.x, pt4.y) +
484  PI - area_arc(pt3.x, pt3.y, pt4.x, pt4.y, 1.));
485  }
486  else
487  {
488  area = (area_triangle(x1, y1, pt3.x, pt3.y, pt4.x, pt4.y) +
489  area_arc(pt3.x, pt3.y, pt4.x, pt4.y, 1.));
490  }
491  }
492  else
493  {
494  /* ensure that pt1 is the point closest to (x2, y2) */
495  if (((pt2.x-x2)*(pt2.x-x2) + (pt2.y-y2)*(pt2.y-y2)) <
496  ((pt1.x-x2)*(pt1.x-x2) + (pt1.y-y2)*(pt1.y-y2)))
497  swap_point(&pt1, &pt2);
498 
499  area = (area_triangle(x1, y1, pt3.x, pt3.y, pt1.x, pt1.y) +
500  area_triangle(x1, y1, pt1.x, pt1.y, pt2.x, pt2.y) +
501  area_triangle(x1, y1, pt2.x, pt2.y, pt4.x, pt4.y) +
502  area_arc(pt1.x, pt1.y, pt3.x, pt3.y, 1.) +
503  area_arc(pt2.x, pt2.y, pt4.x, pt4.y, 1.));
504  }
505  }
506  else
507  {
508  inter = circle_segment(x1, y1, x2, y2);
509  pt1 = inter.p1;
510  pt2 = inter.p2;
511  inter = circle_segment(x2, y2, x3, y3);
512  pt3 = inter.p1;
513  pt4 = inter.p2;
514  inter = circle_segment(x3, y3, x1, y1);
515  pt5 = inter.p1;
516  pt6 = inter.p2;
517  if (pt1.x <= 1.)
518  {
519  xp = 0.5 * (pt1.x + pt2.x);
520  yp = 0.5 * (pt1.y + pt2.y);
521  area = (triangle_unitcircle_overlap(x1, y1, x3, y3, xp, yp) +
522  triangle_unitcircle_overlap(x2, y2, x3, y3, xp, yp));
523  }
524  else if (pt3.x <= 1.)
525  {
526  xp = 0.5 * (pt3.x + pt4.x);
527  yp = 0.5 * (pt3.y + pt4.y);
528  area = (triangle_unitcircle_overlap(x3, y3, x1, y1, xp, yp) +
529  triangle_unitcircle_overlap(x2, y2, x1, y1, xp, yp));
530  }
531  else if (pt5.x <= 1.)
532  {
533  xp = 0.5 * (pt5.x + pt6.x);
534  yp = 0.5 * (pt5.y + pt6.y);
535  area = (triangle_unitcircle_overlap(x1, y1, x2, y2, xp, yp) +
536  triangle_unitcircle_overlap(x3, y3, x2, y2, xp, yp));
537  }
538  else /* no intersections */
539  {
540  if (in_triangle(0., 0., x1, y1, x2, y2, x3, y3))
541  return PI;
542  else
543  return 0.;
544  }
545  }
546  return area;
547 }
548 
549 /* exact overlap between a rectangle defined by (xmin, ymin, xmax,
550  ymax) and an ellipse with major and minor axes rx and ry
551  respectively and position angle theta. */
552 static double ellipoverlap(double xmin, double ymin, double xmax, double ymax,
553  double a, double b, double theta)
554 {
555  double cos_m_theta, sin_m_theta, scale;
556  double x1, y1, x2, y2, x3, y3, x4, y4;
557 
558  cos_m_theta = cos(-theta);
559  sin_m_theta = sin(-theta);
560 
561  /* scale by which the areas will be shrunk */
562  scale = a * b;
563 
564  /* Reproject rectangle to a frame in which ellipse is a unit circle */
565  x1 = (xmin * cos_m_theta - ymin * sin_m_theta) / a;
566  y1 = (xmin * sin_m_theta + ymin * cos_m_theta) / b;
567  x2 = (xmax * cos_m_theta - ymin * sin_m_theta) / a;
568  y2 = (xmax * sin_m_theta + ymin * cos_m_theta) / b;
569  x3 = (xmax * cos_m_theta - ymax * sin_m_theta) / a;
570  y3 = (xmax * sin_m_theta + ymax * cos_m_theta) / b;
571  x4 = (xmin * cos_m_theta - ymax * sin_m_theta) / a;
572  y4 = (xmin * sin_m_theta + ymax * cos_m_theta) / b;
573 
574  /* Divide resulting quadrilateral into two triangles and find
575  intersection with unit circle */
576  return scale * (triangle_unitcircle_overlap(x1, y1, x2, y2, x3, y3) +
577  triangle_unitcircle_overlap(x1, y1, x4, y4, x3, y3));
578 }
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Mon Aug 15 2022 04:04:03 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.