Как можно посчитать площадь сложной фигуры?

Question

[kate]

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

Увеличила левый верхний угол и закрасила непокрытые участки для наглядности:


Из окружностей известны только их радиусы, но можно по графику точки пересечения найти.

Не понимаю, как можно посчитать через интегралы, если не задана f(x).
На всякий случай прикрепляю построение окружности:

R = 2;
X = 2;
Y = 2;
phi = linspace(0,2*pi);
x = X + R*cos(phi);
y = Y+ R*sin(phi);    
plot(x, y,'r', X,Y, 'r*');

Comments

[d'Ivan]

Катечка, снимай маску.

[Wataru]

Ох и сомневаюсь я, если матлаб так умеет.

[kate]

Wataru, Да матлаб, наверное, не столь важен, как сам подсчет таких фигур(

[d'Ivan]

Нужно подсчитать именно аналитически?

[kate]

Роман Мирр, да хоть как-то понять, как это можно посчитать, если фигура состоит из дуг.

[Dmtm]

не уверен про матлаб, но если достаточно точности на уровне пикселей, то завести массив соответствующий прямоугольнику и заносить туда пиксели каждого круга (через принадлежность точки), потом найти пустые позиции в массиве

[Lander]

Можно использовать метод Монте-Карло. :)

Answer 1

[Илья]

Код на Python

import numpy as np
from skimage import draw

height, width = 700, 1300
image = np.zeros((height, width), dtype=np.uint8)

# генерируем рандомные окружности
N = 250
np.random.seed(10)
xy = np.random.randint((0, 0), (height, width), (N, 2))
radii = np.random.randint(40, 60, N)


# закрашиваем окружности (цветом 1, фон 0)
for x, y, r in np.column_stack((xy, radii)):
    rr, cc = draw.disk((x, y), r, shape=(height, width))
    image[rr, cc] = 1


# подсчитываем количество ненулевых пикелей
# или просто суммируем, если единицы
# и получаем долю от всех пикселей
area_percent = np.count_nonzero(image)/(height * width)
print(f'площадь занимаемая окружностями {area_percent*100}%')

Вывод: площадь занимаемая окружностями 88.35087912087913%

Можно увеличить разрешение изображения, отмасштабировав окружности и получить высокую точность площади

Answer 2

[Wataru]

Как написал Илья, можно растеризовать картинку и подсчитать незакрашенные пиксели. Это простой в реализации метод, но, если нужна высокая точность, он будет очень медленным и требовательным к памяти (экспоненциальное время и память от количества нужных знаков). Однако, есть алгоритм, который позволяет получать площадь в символьном виде (типа 3/2+5113/11*sqrt(31)). Ну, или любое нужное вам количество знаков без огромного потребления памяти или вычислений (метод полиномиален от количества вычисляемых знаков).

Это некоторое обобщение метода выделения областей или нахождения граней планарного графа.

Сначала вам надо пересечь все пары окружностей и окружности с прямоугольникм, назначить всем уникальным точкам уникальные номера. Потом надо на каждой окружности составить список всех точек пересечений на ней и отсортировать его вдоль окружности (допустим, против часовой стрелки). Тут же надо добавить самую левую и самую правую точку на окружности (X+-R, Y) - это поможет считать площадь позже. также надо составить список всех точек на прямоугольнике, отсортированных также против часовой стрелки. В этот список надо включить вершины прямоугольника.

Теперь надо построить граф: вершины в нем - точки пересечения, вершины прямоугольника, горизональные точки на окружностях, а ребра - дуги окружностей и отрезки прямоугольника. Для этого мы и сортировали точки против часовой стрелки. Тупо для каждой окружности/прямоугольника добавляем в каждой вершине ребра до следующуй и предыдущей вершины в списке.

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

После этой сортировки вы в каждой вершине имеете все ребра отсортированные по углу. Теперь можно для каждого исходяжего из вершины ребра найти следующее по часовой стрелке. Это то ребро, которое получится, если стоя в точке пересечения пойти не по этому ребру, а по тому, что идет чуть-чуть левее:

вот картинка mad paint skillz

Также, если не очевидно, граф-ориентированный. И по время построения вам надо для каждого ребра запомнить обратное.

Теперь можно выделить все области на рисунке: надо, как в алгоритме по ссылке выше, взять любое пока не обработанное ребро и в цикле, пока не вернетесь в него же, переходить к обратному ребру и от него к следущему по углу из вершины. Таким образом симулируется обход области с поворотом как можно левее. Или как-будто вы находитесь на плоскости внутри области и обходите ее приложив правую руку к стене.

В итоге каждая замкнутая область будет обойдена против часовой стрелки. Кроме одного исключения - вся картинка будет обойдена по границе по часовой стрелке, но это можно легко проверить через векторные произведения соседних дуг (касательных) и выкинуть эту область из рассмотрения. У вас будет список дуг/отрезков прямоугольника. Теперь осталось понять, какие области не покрыты окружностями. Во-первых, если область имеет выпуклую дугу (дуга против часовой стрелки вдоль окружности), то область 100% внутри окружности. Во-вторых, надо взять любую вершину области и проверить, не лежит ли она строго внутри какой-либо окружности.

Если область оказалась не покрыта, то можно взять ее площадь и прибавить к ответу. Тут нам помогает, что мы нанесли горизонтальные точки на окружности - все ребра области или сторого вертикальные отрезки или идут горизонтально, как график функции. А значит можно взять интеграл по x (вы знаете, какой окружности принадлежит ребро, знаете границы по x и верхняя или нижняя грань это).

Вот так и считается ответ. Все координаты можно считать в символьном виде - они будут вида a/b+sqrt(c)/d. Вектора касательных будут такого же типа. Такие числа можно перемножать, складывать, вычитать и сравнивать. Правда, после операций вы получите множество радикалов, но тут нет проблем, потому что сравнивать результаты вам уже не придется. Ну, или тут можно использовать числа с плавающей запятой и спокойно получить 12-14 знаков точности.