Площадь пересечения круга и прямоугольника
https://stackoverflow.com/questions/622287
-
05-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Я ищу быстрый способ определить площадь пересечения прямоугольника и круга (мне нужно выполнить миллионы таких вычислений).
Особым свойством является то, что во всех случаях круг и прямоугольник всегда имеют 2 точки пересечения.
Решение
Учитывая 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
Чтобы вычислить эти площади:
Большинство точек, которые вам понадобятся, можно найти, найдя пересечение линии и окружности
Области, которые необходимо вычислить, можно найти с помощью вычисление площади кругового сегмента и вычисление площади треугольника.
Вы можете определить, находится ли вершина внутри круга, посчитав, меньше ли ее расстояние от центра радиуса.
Другие советы
Я понимаю, что на этот вопрос был дан ответ некоторое время назад, но я решаю ту же проблему и не смог найти готовое работоспособное решение, которое я мог бы использовать.Обратите внимание, что мои коробки ось выровнена, это не совсем указано в ФП.Решение, приведенное ниже, является полностью общим и будет работать для любого количества перекрестков (а не только для двух).Обратите внимание: если ваши блоки не выровнены по осям (но все же блоки имеют прямые углы, а не общие четверные), вы можете воспользоваться тем, что круги круглые, повернуть координаты всего так, чтобы прямоугольник оказался выровнен по оси, и затем используйте этот код.
Я хочу использовать интеграцию — это хорошая идея.Начнем с написания очевидной формулы для построения окружности:
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 (используйте приведенное выше уравнение)
- прямоугольник находится под осью 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, и удвоенное число (и еще 1 добавление) для блоков, которые пересекаются.Обратите внимание, что деления производятся по радиусу и могут быть сокращены до двух делений и набора умножений.Если рассматриваемый прямоугольник пересекал ось X, но не пересекал ось Y, вы могли бы повернуть все на pi/2 и делаем расчет в первоначальной стоимости.
Я использую его как справочник для отладки растеризатора сглаженных кругов с точностью до субпикселя.Это чертовски медленно :), мне нужно вычислить площадь пересечения каждого пикселя в ограничивающей рамке круга с кругом, чтобы получить альфу.Вы можете видеть, что это работает, и никаких числовых артефактов не появляется.На рисунке ниже показан участок группы кругов с радиусом, увеличивающимся от 0,3 пикселей до примерно 6 пикселей, расположенных по спирали.
Надеюсь, это не плохой тон — опубликовать ответ на такой старый вопрос.Я просмотрел приведенные выше решения и разработал алгоритм, который похож на первый ответ Дэниэлса, но немного сложнее.
Короче говоря, предположим, что вся площадь находится в прямоугольнике, вычтите четыре сегмента во внешних полуплоскостях, затем добавьте все площади в четырех внешних квадрантах, отбрасывая при этом тривиальные случаи.
псевдокод (мой фактический код составляет всего ~ 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
Кстати, последняя формула площади круга, содержащегося в плоском квадранте, легко выводится как сумма сегмента окружности, двух прямоугольных треугольников и прямоугольника.
Наслаждаться.
Ниже описано, как вычислить перекрывающуюся область между кругом и прямоугольником, в которой центр круга находится за пределами прямоугольника. Другие случаи могут быть сведены к этой проблеме.
Площадь можно рассчитать путем интегрирования уравнения окружности y = sqrt [a ^ 2 - (xh) ^ 2] + k , где a - радиус, (h, k) - центр круга, найти область под кривой. Вы можете использовать компьютерную интеграцию, где область делится на множество маленьких прямоугольников и вычислять их сумму, или просто использовать здесь закрытую форму.
А вот источник C #, реализующий концепцию выше. Обратите внимание, что существует особый случай, когда указанный х лежит вне границ круга. Я просто использую простой обходной путь (который не дает правильных ответов во всех случаях)
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 года : Мухобойка . Вы также можете нажать на ссылку для оценки, чтобы загрузить исходный код решения.
Спасибо за ответы,
Я забыл упомянуть, что оценки площади было достаточно. Тот; Поэтому, в конце концов, после просмотра всех вариантов, я пошел с оценкой Монте-Карло, где я генерирую случайные точки в круге и проверяю, находятся ли они в рамке.
В моем случае это, вероятно, более производительный. (У меня есть сетка, расположенная над кругом, и я должен измерить отношение круга, принадлежащего каждой из ячеек сетки.)
Спасибо
Возможно, вы сможете использовать ответ для этот вопрос , где задается область пересечения круга и треугольника. Разделите ваш прямоугольник на два треугольника и используйте алгоритмы, описанные там.
Другой способ - нарисовать линию между двумя точками пересечения, это разделит вашу область на многоугольник (3 или 4 ребра) и круговой сегмент , для обоих вы сможете легче находить библиотеки или выполнять вычисления самостоятельно.
Вот еще одно решение проблемы:
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)));
}
