Я ищу быстрый способ определить область пересечения прямоугольника и круга (мне нужно выполнить миллионы этих вычислений).
Конкретное свойство состоит в том, что во всех случаях круг и прямоугольник всегда имеют 2 точки пересечения.
Ответ 3
Я понимаю, что это было ответили некоторое время назад, но я решаю ту же проблему, и я не мог найти готовое решение, которое я мог бы использовать. Обратите внимание, что мои поля выровнены по оси, это не совсем определено 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). Теперь мы можем вычислить площадь пересечения полукруга и бесконечно высокий квадрат, стоящий на оси х (центр круга также лежит на оси х), интегрируя по 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)
Теперь мы можем понять, что происходит. Итак, как вычислить площадь пересечения конечной коробки, пересекающей единичный круг над осью х:
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 (используйте приведенное выше уравнение)
- поле находится ниже оси 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 добавок для ящиков, которые не пересекают ось x, и double из этого (и еще 1 add) для ящиков, которые делают. Обратите внимание, что дивизии имеют радиус и могут быть сведены к двум делениям и кучке умножений. Если поле, о котором идет речь, пересекало ось x, но не пересекало ось y, вы можете повернуть все на pi/2 и выполнить расчет в первоначальной стоимости.
Я использую его как ссылку на отладочный субпиксельный точный антиалиасированный круговой растеризатор. Это медленно, как ад:), мне нужно вычислить площадь пересечения каждого пикселя в ограничивающей рамке круга с кругом, чтобы получить альфу. Вы можете видеть, что он работает, и нет никаких числовых артефактов. На рисунке ниже представлен график группы кругов с радиусом, увеличивающимся от 0,3 px до 6 px, выложенным по спирали.
![круги]()
