Площадь пересечения круга и прямоугольника

Учитывая 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

Чтобы рассчитать эти площади:

• Большинство точек, которые вам понадобятся, можно найти, найдя пересечение линии и круга.

• Области, которые вам нужно вычислить, можно найти, вычислив площадь кругового сегмента и вычисление площади треугольника.

• Вы можете определить, находится ли вершина внутри круга, вычислив, меньше ли ее расстояние от центра радиуса.

Я понимаю, что на это был дан ответ некоторое время назад, но я решаю ту же проблему и не могу найти готовое работоспособное решение, которое я мог бы использовать. Обратите внимание, что мои блоки выровнены по оси, это не совсем указано OP. Приведенное ниже решение является полностью общим и будет работать для любого количества перекрестков (а не только для двух). Обратите внимание, что если ваши прямоугольники не выровнены по оси (но все же прямоугольники с прямыми углами, а не общие четырехугольники), вы можете воспользоваться округлыми кругами, повернуть все координаты так, чтобы прямоугольник оказался выровненным по оси, и затем используйте этот код.

Я хочу использовать интеграцию - это хорошая идея. Начнем с написания очевидной формулы построения круга:

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, и мы уже сказали unit, что означает 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 add) для блоков, которые пересекают. Обратите внимание, что деления делаются по радиусу и могут быть сокращены до двух делений и множества умножений. Если рассматриваемый прямоугольник пересекает ось 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 #, реализующий описанную выше концепцию. Обратите внимание, что есть особый случай, когда указанный x лежит за пределами круга. Я просто использую простой обходной путь (который не во всех случаях дает правильные ответы)

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 проблема: Мухобойка. Вы также можете щелкнуть ссылки оценки, чтобы загрузить исходный код решения.

Спасибо за ответы,

Я забыл упомянуть, что оценок площадей было достаточно. Тот; s почему, в конце концов, просмотрев все варианты, я пошел с оценкой Монте-Карло, где я генерирую случайные точки в круге и проверяю, находятся ли они в коробке.

В моем случае это, вероятно, более эффективно. (У меня есть сетка, размещенная над кругом, и я должен измерить соотношение круга, принадлежащего каждой из ячеек сетки.)

Спасибо

Возможно, вы можете использовать ответ на этот вопрос, где задается площадь пересечения круга и треугольника. Разделите прямоугольник на два треугольника и используйте описанные там алгоритмы.

Другой способ - провести линию между двумя точками пересечения, при этом ваша область разбивается на многоугольник (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)));
}