Я понимаю, что на это ответили некоторое время назад, но я решаю ту же проблему, и я не смог найти готового решения, которое можно было бы использовать. Обратите внимание, что мои поля выровнены по оси , это не совсем определено OP. Решение, приведенное ниже, является полностью общим и будет работать для любого количества пересечений (не только для двух). Обратите внимание, что если ваши прямоугольники не выровнены по оси (но все же прямоугольники с прямыми углами, а не общие quads ), вы можете воспользоваться круглыми кругами, повернуть координаты всего так, чтобы закончился прямоугольник выровнен по оси и , затем используйте этот код.
Я хочу использовать интеграцию - это хорошая идея. Начнем с написания очевидной формулы для построения круга:
x = center.x + cos(theta) * radius
y = center.y + sin(theta) * radius
^
|
|**### **
| #* # * * x
|# * # * # y
|# * # *
+-----------------------> theta
* # * #
* # * #
* #* #
***###
Это хорошо, но я не могу интегрировать область этого круга по x или y; это разные количества. Я могу интегрировать только по углу theta, получая участки ломтиков пиццы. Не то, что я хочу. Попробуем изменить аргументы:
(x - center.x) / radius = cos(theta) // the 1st equation
theta = acos((x - center.x) / radius) // value of theta from the 1st equation
y = center.y + sin(acos((x - center.x) / radius)) * radius // substitute to the 2nd equation
Это больше похоже на это. Теперь, учитывая диапазон x, я могу интегрировать по y, чтобы получить площадь верхней половины круга. Это справедливо только для x в [center.x - radius, center.x + radius] (другие значения приведут к мнимым результатам), но мы знаем, что область вне этого диапазона равна нулю, так что это легко обрабатывается. Давайте для простоты предположим единичный круг, мы всегда можем подключить центр и радиус обратно:
y = sin(acos(x)) // x in [-1, 1]
y = sqrt(1 - x * x) // the same thing, arguably faster to compute http://www.wolframalpha.com/input/?i=sin%28acos%28x%29%29+
^ y
|
***|*** <- 1
**** | ****
** | **
* | *
* | *
----|----------+----------|-----> x
-1 1
Эта функция действительно имеет интеграл pi/2, поскольку она является верхней половиной единичного круга (площадь полукруга равна pi * r^2 / 2, и мы уже сказали единица , что означает r = 1). Теперь мы можем вычислить площадь пересечения полукруга и бесконечно высокого прямоугольника, стоящего на оси x (центр круга также лежит на оси x), интегрируя по y:
f(x): integral(sqrt(1 - x * x) * dx) = (sqrt(1 - x * x) * x + asin(x)) / 2 + C // http://www.wolframalpha.com/input/?i=sqrt%281+-+x*x%29
area(x0, x1) = f(max(-1, min(1, x1))) - f(max(-1, min(1, x0))) // the integral is not defined outside [-1, 1] but we want it to be zero out there
~ ~
| ^ |
| | |
| ***|*** <- 1
****###|##|****
**|######|##| **
* |######|##| *
* |######|##| *
----|---|------+--|-------|-----> x
-1 x0 x1 1
Это может быть не очень полезно, поскольку бесконечно высокие коробки - это не то, что нам нужно. Нам нужно добавить еще один параметр, чтобы освободить нижний край бесконечно высокого поля:
g(x, h): integral((sqrt(1 - x * x) - h) * dx) = (sqrt(1 - x * x) * x + asin(x) - 2 * h * x) / 2 + C // http://www.wolframalpha.com/input/?i=sqrt%281+-+x*x%29+-+h
area(x0, x1, h) = g(min(section(h), max(-section(h), x1))) - g(min(section(h), max(-section(h), x0)))
~ ~
| ^ |
| | |
| ***|*** <- 1
****###|##|****
**|######|##| **
* +------+--+ * <- h
* | *
----|---|------+--|-------|-----> x
-1 x0 x1 1
Где h - (положительное) расстояние нижнего края нашего бесконечного прямоугольника от оси x. Функция section вычисляет (положительное) положение пересечения единичного круга с горизонтальной линией, заданной y = h, и мы можем определить ее, решив:
sqrt(1 - x * x) = h // http://www.wolframalpha.com/input/?i=sqrt%281+-+x+*+x%29+%3D+h
section(h): (h < 1)? sqrt(1 - h * h) : 0 // if h is 1 or above, then the section is an empty interval and we want the area integral to be zero
^ y
|
***|*** <- 1
**** | ****
** | **
-----*---------+---------*------- y = h
* | *
----||---------+---------||-----> x
-1| |1
-section(h) section(h)
Теперь мы можем начать действовать. Итак, как рассчитать площадь пересечения конечного прямоугольника, пересекающего единичную окружность над осью x:
area(x0, x1, y0, y1) = area(x0, x1, y0) - area(x0, x1, y1) // where x0 <= x1 and y0 <= y1
~ ~ ~ ~
| ^ | | ^ |
| | | | | |
| ***|*** | ***|***
****###|##|**** ****---+--+**** <- y1
**|######|##| ** ** | **
* +------+--+ * <- y0 * | *
* | * * | *
----|---|------+--|-------|-----> x ----|---|------+--|-------|-----> x
x0 x1 x0 x1
^
|
***|***
****---+--+**** <- y1
**|######|##| **
* +------+--+ * <- y0
* | *
----|---|------+--|-------|-----> x
x0 x1
Это хорошо. Так как насчет коробки, которая не выше оси х? Я бы сказал, что не все коробки. Возникают три простых случая:
- поле выше оси x (используйте приведенное выше уравнение)
- поле находится ниже оси x (переверните знак координат y и используйте приведенное выше уравнение)
- прямоугольник пересекает ось x (разделите прямоугольник на верхнюю и нижнюю половину, рассчитайте площадь обоих, используя приведенное выше, и суммируйте их)
Достаточно легко? Давайте напишем некоторый код:
float section(float h, float r = 1) // returns the positive root of intersection of line y = h with circle centered at the origin and radius r
{
assert(r >= 0); // assume r is positive, leads to some simplifications in the formula below (can factor out r from the square root)
return (h < r)? sqrt(r * r - h * h) : 0; // http://www.wolframalpha.com/input/?i=r+*+sin%28acos%28x+%2F+r%29%29+%3D+h
}
float g(float x, float h, float r = 1) // indefinite integral of circle segment
{
return .5f * (sqrt(1 - x * x / (r * r)) * x * r + r * r * asin(x / r) - 2 * h * x); // http://www.wolframalpha.com/input/?i=r+*+sin%28acos%28x+%2F+r%29%29+-+h
}
float area(float x0, float x1, float h, float r) // area of intersection of an infinitely tall box with left edge at x0, right edge at x1, bottom edge at h and top edge at infinity, with circle centered at the origin with radius r
{
if(x0 > x1)
std::swap(x0, x1); // this must be sorted otherwise we get negative area
float s = section(h, r);
return g(max(-s, min(s, x1)), h, r) - g(max(-s, min(s, x0)), h, r); // integrate the area
}
float area(float x0, float x1, float y0, float y1, float r) // area of the intersection of a finite box with a circle centered at the origin with radius r
{
if(y0 > y1)
std::swap(y0, y1); // this will simplify the reasoning
if(y0 < 0) {
if(y1 < 0)
return area(x0, x1, -y0, -y1, r); // the box is completely under, just flip it above and try again
else
return area(x0, x1, 0, -y0, r) + area(x0, x1, 0, y1, r); // the box is both above and below, divide it to two boxes and go again
} else {
assert(y1 >= 0); // y0 >= 0, which means that y1 >= 0 also (y1 >= y0) because of the swap at the beginning
return area(x0, x1, y0, r) - area(x0, x1, y1, r); // area of the lower box minus area of the higher box
}
}
float area(float x0, float x1, float y0, float y1, float cx, float cy, float r) // area of the intersection of a general box with a general circle
{
x0 -= cx; x1 -= cx;
y0 -= cy; y1 -= cy;
// get rid of the circle center
return area(x0, x1, y0, y1, r);
}
И некоторые базовые юнит-тесты:
printf("%f\n", area(-10, 10, -10, 10, 0, 0, 1)); // unit circle completely inside a huge box, area of intersection is pi
printf("%f\n", area(-10, 0, -10, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(0, 10, -10, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(-10, 10, -10, 0, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(-10, 10, 0, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2
printf("%f\n", area(0, 1, 0, 1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(0, -1, 0, 1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(0, -1, 0, -1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(0, 1, 0, -1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4
printf("%f\n", area(-.5f, .5f, -.5f, .5f, 0, 0, 10)); // unit box completely inside a huge circle, area of intersection is 1
printf("%f\n", area(-20, -10, -10, 10, 0, 0, 1)); // huge box completely outside a circle (left), area of intersection is 0
printf("%f\n", area(10, 20, -10, 10, 0, 0, 1)); // huge box completely outside a circle (right), area of intersection is 0
printf("%f\n", area(-10, 10, -20, -10, 0, 0, 1)); // huge box completely outside a circle (below), area of intersection is 0
printf("%f\n", area(-10, 10, 10, 20, 0, 0, 1)); // huge box completely outside a circle (above), area of intersection is 0
Вывод этого:
3.141593
1.570796
1.570796
1.570796
1.570796
0.785398
0.785398
0.785398
0.785398
1.000000
-0.000000
0.000000
0.000000
0.000000
Что мне кажется правильным. Вы можете захотеть встроить функции, если не доверяете своему компилятору сделать это за вас.
Используется 6 sqrt, 4 asin, 8 div, 16 mul и 17 add для боксов, которые не пересекают ось x, и удвоение этого значения (и еще 1 add) для боксов, которые это делают. Обратите внимание, что деления по радиусу и могут быть уменьшены до двух делений и нескольких умножений. Если рассматриваемый прямоугольник пересекает ось x, но не пересекает ось y, вы можете повернуть все на pi/2 и выполнить расчет по первоначальной стоимости.
Я использую его как ссылку для отладки растрового сглаженного круга с точностью до субпикселя. Это чертовски медленно :), мне нужно вычислить площадь пересечения каждого пикселя в ограничительной рамке круга с кругом, чтобы получить альфа. Вы можете видеть, что это работает, и никакие числовые артефакты, кажется, не появляются. На рисунке ниже изображен ряд кругов с радиусом, увеличивающимся от 0,3 пикселей до 6 пикселей, расположенных по спирали.
