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 */
14static 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 */
25static 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 */
35static 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 */
99static 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/*
143Start of new circular overlap routine that might be faster.
144
145double 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
177typedef struct
178{
179 double x, y;
180} point;
181
182
183typedef struct
184{
185 point p1, p2;
186} intersections;
187
188static void swap(double *a, double *b)
189{
190 double temp;
191 temp = *a;
192 *a = *b;
193 *b = temp;
194}
195
196static 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 */
205static 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 */
215static 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 */
230static 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. */
301static 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. */
327intersections 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. */
365static 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;
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.);
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. */
552static 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-2024 The KDE developers.
Generated on Fri May 17 2024 11:48:26 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.