Я ищу быстрый способ определить область пересечения прямоугольника и круга (мне нужно выполнить миллионы этих вычислений).
Конкретное свойство состоит в том, что во всех случаях круг и прямоугольник всегда имеют 2 точки пересечения.
Площадь пересечения между кругом и прямоугольником
Ответ 1
Учитывая 2 точки пересечения:
0 вершин находится внутри круга: область круговой сегмент
XXXXX -------------------
X X X X Circular segment
X X XX XX
+-X-------X--+ XXXXXXXX
| X X |
| XXXXX |
1 вершина находится внутри круга: сумма площадей кругового сегмента и треугольника.
XXXXX XXXXXXXXX
X X Triangle ->X _-X
X X X _- X
X +--X--+ X _- X <- Circular segment
X | X | X- XXX
XXXXX | XXXX
| |
2 вершины находятся внутри круга: сумма площади двух треугольников и кругового сегмента
XXXXX +------------X
X X | _--'/'X
X +--X--- Triangle->| _-- / X
X | X |_-- /XX <- Circular segment
X +-X---- +-------XX
XXXXX Triangle^
3 вершины находятся внутри круга: площадь прямоугольника минус площадь треугольника плюс площадь кругового сегмента
XXXXX
X +--X+ XXX
X | X -------XXX-----+ <- Triangle outside
X | |X Rect ''. XXX |
X +---+X ''. XX|
X X ''. X <- Circular segment inside
X X ^|X
X X | X
XXXXX
Чтобы вычислить эти области:
-
Большинство пунктов, которые вам нужно использовать, можно найти, найдя пересечение линии и круга
-
Области, которые нужно вычислить, можно найти для вычисления площади кругового сегмента и вычисляет площадь треугольника.
-
Вы можете определить, находится ли вершина внутри круга, вычисляя, если его расстояние от центра меньше радиуса.
Ответ 2
Я надеюсь, что его неплохая форма для публикации ответа на такой старый вопрос. Я просмотрел приведенные выше решения и разработал алгоритм, похожий на первый ответ Дэниелса, но немного лучше.
Короче говоря, предположим, что полная площадь находится в прямоугольнике, вычитайте четыре сегмента во внешних полуплоскостях, затем добавьте любые области в четырех внешних квадрантах, отбросив тривиальные случаи на этом пути.
pseudocde (мой фактический код - всего ~ 12 строк..)
find the signed (negative out) normalized distance from the circle center
to each of the infinitely extended rectangle edge lines,
ie.
d_1=(xcenter-xleft)/r
d_2=(ycenter-ybottom)/r
etc
for convenience order 1,2,3,4 around the edge. If the rectangle is not
aligned with the cartesian coordinates this step is more complicated but
the remainder of the algorithm is the same
If ANY d_i <=- 1 return 0
if ALL d_i >= 1 return Pi r^2
this leave only one remaining fully outside case: circle center in
an external quadrant, and distance to corner greater than circle radius:
for each adjacent i,j (ie. i,j=1,2;2,3;3,4;4,1)
if d_i<=0 and d_j <= 0 and d_i^2+d_j^2 > 1 return 0
now begin with full circle area and subtract any areas in the
four external half planes
Area= Pi r^2
for each d_i>-1
a_i=arcsin( d_i ) #save a_i for next step
Area -= r^2/2 (Pi - 2 a_i - sin(2 a_i))
At this point note we have double counted areas in the four external
quadrants, so add back in:
for each adjacent i,j
if d_i < 1 and d_j < 1 and d_i^2+d_j^2 < 1
Area += r^2/4 (Pi- 2 a_i - 2 a_j -sin(2 a_i) -sin(2 a_j) + 4 sin(a_i) sin(a_j))
return Area
Кстати, последняя формула для области круга, содержащейся в плоском квадранте, легко выводится как сумма кругового сегмента, двух правых треугольников и прямоугольника.
Enjoy.
Ответ 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, выложенным по спирали.
Ответ 4
Ниже описано, как рассчитать перекрывающуюся область между кругом и прямоугольником, где центр круга лежит вне прямоугольника. Другие проблемы могут быть сведены к этой проблеме.
Площадь можно рассчитать, интегрируя уравнение окружности y = sqrt [a ^ 2 - (x-h) ^ 2] + k, где a - радиус, (h, k) - центр круга, чтобы найти область под кривой. Вы можете использовать компьютерную интеграцию, где область разделена на множество маленьких прямоугольников и вычисляет их сумму или просто использует закрытую форму здесь.
И вот источник С#, реализующий концепцию выше. Заметим, что существует специальный случай, когда указанный х лежит за пределами круга. Я просто использую простой способ обхода (который не дает правильных ответов во всех случаях)
public static void RunSnippet()
{
// test code
double a,h,k,x1,x2;
a = 10;
h = 4;
k = 0;
x1 = -100;
x2 = 100;
double r1 = Integrate(x1, a, h, k);
double r2 = Integrate(x2, a, h, k);
Console.WriteLine(r2 - r1);
}
private static double Integrate(double x, double a,double h, double k)
{
double a0 = a*a - (h-x)*(h-x);
if(a0 <= 0.0){
if(k == 0.0)
return Math.PI * a * a / 4.0 * Math.Sign(x);
else
throw new Exception("outside boundaries");
}
double a1 = Math.Sqrt(a*a - (h-x)*(h-x)) * (h-x);
double area = 0.5 * Math.Atan(a1 / ((h-x)*(h-x) - a*a))*a*a - 0.5 * a1 + k * x;
return area;
}
Примечание.. Эта проблема очень похожа на ту, что находится в Ошибка в Google Code Jam 2008: Fly Swatter, Вы также можете щелкнуть ссылку на ссылку, чтобы загрузить исходный код решения.
Ответ 5
Спасибо за ответы,
Я забыл упомянуть, что оценки области были достаточно. Что; почему в конце концов, посмотрев все варианты, я пошел с оценкой monte-carlo, где я генерирую случайные точки в круге и проверяю, находятся ли они в коробке.
В моем случае это, скорее всего, более результативно. (У меня есть сетка, расположенная над кругом, и мне нужно измерить отношение круга, принадлежащего каждой из ячеек сетки.)
Спасибо
Ответ 6
Возможно, вы можете использовать ответ на этот вопрос, где задается область пересечения между кругом и треугольником. Разделите прямоугольник на два треугольника и используйте описанные там алгоритмы.
Другой способ - провести линию между двумя точками пересечения, это разбивает вашу область на многоугольник (3 или 4 ребра) и круглый сегмент, для обоих вы должны иметь возможность находить библиотеки проще или сами выполнять вычисления.
Ответ 7
Вот еще одно решение проблемы:
public static bool IsIntersected(PointF circle, float radius, RectangleF rectangle)
{
var rectangleCenter = new PointF((rectangle.X + rectangle.Width / 2),
(rectangle.Y + rectangle.Height / 2));
var w = rectangle.Width / 2;
var h = rectangle.Height / 2;
var dx = Math.Abs(circle.X - rectangleCenter.X);
var dy = Math.Abs(circle.Y - rectangleCenter.Y);
if (dx > (radius + w) || dy > (radius + h)) return false;
var circleDistance = new PointF
{
X = Math.Abs(circle.X - rectangle.X - w),
Y = Math.Abs(circle.Y - rectangle.Y - h)
};
if (circleDistance.X <= (w))
{
return true;
}
if (circleDistance.Y <= (h))
{
return true;
}
var cornerDistanceSq = Math.Pow(circleDistance.X - w, 2) +
Math.Pow(circleDistance.Y - h, 2);
return (cornerDistanceSq <= (Math.Pow(radius, 2)));
}