Сначала я напомню нам, как найти площадь многоугольника. Как только мы это сделаем, алгоритм для нахождения пересечения между многоугольником и окружностью должен быть легок для понимания.
Как найти площадь многоугольника
Давайте посмотрим на случай треугольника, потому что там появляется вся необходимая логика. Давайте предположим, что у нас есть треугольник с вершинами (x1, y1), (x2, y2) и (x3, y3) при обходе треугольника против часовой стрелки, как показано на следующем рисунке:
Тогда вы можете вычислить площадь по формуле
A = (x1 y2 + x2 y3 + x3 y1 - x2y1- x3 y2 - x1y3) /2.
Чтобы понять, почему эта формула работает, давайте изменим ее так, чтобы она была в форме
A = (x1 y2 - x2 y1) / 2 + (x2 y3 - x3 y2) / 2 + (x3 y1 - x1y3) /2.
Теперь первый член - это следующая область, которая в нашем случае является положительной:
Если неясно, что площадь зеленой области действительно (x1 y2 - x2 y1) / 2, тогда читайте this .
Второй член - это область, которая снова положительна:
И третья область показана на следующем рисунке. На этот раз область отрицательная
Сложив эти три, мы получим следующую картинку
Мы видим, что зеленая область, которая была за пределами треугольника, отменяется красной областью, так что чистая область является просто областью треугольника, и это показывает, почему наша формула была верной в этом случае.
То, что я сказал выше, было интуитивным объяснением того, почему формула площади была правильной. Более строгое объяснение будет состоять в том, чтобы заметить, что когда мы вычисляем площадь от ребра, получаемая площадь является той же самой областью, которую мы получили бы из интегрирования r ^ 2dθ / 2, поэтому мы эффективно интегрируем r ^ 2dθ / 2 вокруг границы многоугольника, и по теореме Стокса это дает тот же результат, что и интегрирование rdrdθ по области, ограниченной многоугольником. Поскольку интегрирование rdrdθ по области, ограниченной многоугольником, дает площадь, мы заключаем, что наша процедура должна правильно дать площадь.
Площадь пересечения окружности с многоугольником
Теперь давайте обсудим, как найти область пересечения окружности радиуса R с многоугольником, как показано на следующем рисунке:
Мы заинтересованы в поиске области зеленого региона. Мы можем, как и в случае с единственным полигоном, разбить наш расчет на нахождение области для каждой стороны полигона, а затем сложить эти области.
Наша первая область будет выглядеть так:
Вторая область будет выглядеть
И третья область будет
Опять же, первые две области в нашем случае положительные, а третья - отрицательная. Надеемся, что отмены сработают, так что чистая область - это действительно область, которая нас интересует. Посмотрим.
Действительно, сумма областей будет интересующей нас областью.
Опять же, мы можем дать более строгое объяснение того, почему это работает. Позвольте мне быть областью, определенной пересечением, и позвольте P быть многоугольником. Затем из предыдущего обсуждения мы знаем, что мы хотим вычислить интеграл от r ^ 2dθ / 2 вокруг границы I. Однако это трудно сделать, потому что это требует нахождения пересечения.
Вместо этого мы сделали интеграл по многоугольнику. Мы интегрировали max (r, R) ^ 2 dθ / 2 через границу многоугольника. Чтобы понять, почему это дает правильный ответ, давайте определим функцию π, которая переводит точку в полярных координатах (r, θ) в точку (max (r, R), θ). Не следует путать обращение к координатным функциям π (r) = max (r, R) и π (θ) = θ. Затем мы интегрировали π (r) ^ 2 dθ / 2 по границе многоугольника.
С другой стороны, поскольку π (θ) = θ, это то же самое, что интегрирование π (r) ^ 2 dπ (θ) / 2 через границу многоугольника.
Сейчас занимаюсь чаПри изменении переменной мы находим, что получили бы тот же ответ, если бы интегрировали r ^ 2 dθ / 2 через границу π (P), где π (P) - это образ P под π.
Снова используя теорему Стокса, мы знаем, что интегрирование r ^ 2 dθ / 2 по границе π (P) дает нам площадь π (P). Другими словами, он дает тот же ответ, что и интегрирование dxdy по π (P).
Используя снова изменение переменной, мы знаем, что интегрирование dxdy по π (P) такое же, как интегрирование Jdxdy по P, где J - якобиан π.
Теперь мы можем разбить интеграл Jdxdy на две области: часть в круге и часть вне круга. Теперь π оставляет точки в одной окружности, так что J = 1 там, поэтому вклад этой части P является площадью той части P, которая лежит в окружности, то есть площадью пересечения. Вторая область - область вне круга. Там J = 0, так как π сворачивает эту часть до границы круга.
Таким образом, мы действительно вычисляем площадь пересечения.
Теперь, когда мы относительно уверены, что концептуально знаем, как найти область, давайте поговорим более конкретно о том, как вычислить вклад от одного сегмента. Давайте начнем с рассмотрения сегмента, который я назову «стандартной геометрией». Это показано ниже.
В стандартной геометрии край идет горизонтально слева направо. Он описывается тремя числами: xi, координата x, где начинается ребро, xf, координата x, где заканчивается ребро, и y, координата y ребра.
Теперь мы видим, что если | y |
Площадь области 2 - это площадь треугольника. Тем не менее, мы должны быть осторожны с признаком. Мы хотим, чтобы отображаемая область была положительной, поэтому мы скажем, что эта область - (xint - (-xint)) y / 2.
Еще одна вещь, которую следует иметь в виду, это то, что xi не должно быть меньше -xint, а xf не должно быть больше xint.
Другой случай, который нужно рассмотреть, это | y | > R. Этот случай проще, потому что есть только один кусок, который похож на область 1 на рисунке.
Теперь, когда мы знаем, как вычислить площадь от ребра в стандартной геометрии, остается только описать, как преобразовать любое ребро в стандартную геометрию.
Но это просто простая смена координат. Для некоторых из них с начальной вершиной vi и конечной вершиной vf новый вектор единиц x будет единичным вектором, указывающим от vi к vf. Тогда xi - это просто смещение vi из центра окружности, пунктирной в x, а xf - просто xi плюс расстояние между vi и vf. Между тем, у задается произведением клина на х со смещением vi от центра круга.
код
На этом описание алгоритма завершено, теперь пришло время написать некоторый код. Я буду использовать Java.
Прежде всего, так как мы работаем с кругами, у нас должен быть класс круга
public class Circle {
final Point2D center;
final double radius;
public Circle(double x, double y, double radius) {
center = new Point2D.Double(x, y);
this.radius = radius;
}
public Circle(Point2D.Double center, double radius) {
this(center.getX(), center.getY(), radius);
}
public Point2D getCenter() {
return new Point2D.Double(getCenterX(), getCenterY());
}
public double getCenterX() {
return center.getX();
}
public double getCenterY() {
return center.getY();
}
public double getRadius() {
return radius;
}
}
Для полигонов я буду использовать класс Shape в Java. У Shape есть PathIterator, который я могу использовать для итерации по краям многоугольника.
Теперь для фактической работы. Я отделю логику итерации по краям, помещая ребра в стандартную геометрию и т. Д., От логики вычисления площади, как только это будет сделано. Причина этого в том, что в будущем вы можете захотеть вычислить что-то еще помимо или в дополнение к области, и вы захотите иметь возможность повторно использовать код, имеющий дело с итерацией по краям.
Итак, у меня естьуниверсальный класс, который вычисляет некоторое свойство класса T относительно нашего пересечения многоугольника.
public abstract class CircleShapeIntersectionFinder<T> {
У него есть три статических метода, которые просто помогают вычислить геометрию:
private static double[] displacment2D(final double[] initialPoint, final double[] finalPoint) {
return new double[]{finalPoint[0] - initialPoint[0], finalPoint[1] - initialPoint[1]};
}
private static double wedgeProduct2D(final double[] firstFactor, final double[] secondFactor) {
return firstFactor[0] * secondFactor[1] - firstFactor[1] * secondFactor[0];
}
static private double dotProduct2D(final double[] firstFactor, final double[] secondFactor) {
return firstFactor[0] * secondFactor[0] + firstFactor[1] * secondFactor[1];
}
Есть два поля экземпляра: Circle, который просто хранит копию круга, и currentSquareRadius, который хранит копию квадратного радиуса. Это может показаться странным, но класс, который я использую, на самом деле оборудован для поиска областей целого набора пересечений окружности и многоугольника. Вот почему я называю один из кругов «текущим».
private Circle currentCircle;
private double currentSquareRadius;
Далее следует метод вычисления того, что мы хотим вычислить:
public final T computeValue(Circle circle, Shape shape) {
initialize();
processCircleShape(circle, shape);
return getValue();
}
initialize() и getValue() являются абстрактными. initialize() установит переменную, которая сохраняет общую площадь в ноль, а getValue() просто вернет область. Определение для processCircleShape является
private void processCircleShape(Circle circle, final Shape cellBoundaryPolygon) {
initializeForNewCirclePrivate(circle);
if (cellBoundaryPolygon == null) {
return;
}
PathIterator boundaryPathIterator = cellBoundaryPolygon.getPathIterator(null);
double[] firstVertex = new double[2];
double[] oldVertex = new double[2];
double[] newVertex = new double[2];
int segmentType = boundaryPathIterator.currentSegment(firstVertex);
if (segmentType != PathIterator.SEG_MOVETO) {
throw new AssertionError();
}
System.arraycopy(firstVertex, 0, newVertex, 0, 2);
boundaryPathIterator.next();
System.arraycopy(newVertex, 0, oldVertex, 0, 2);
segmentType = boundaryPathIterator.currentSegment(newVertex);
while (segmentType != PathIterator.SEG_CLOSE) {
processSegment(oldVertex, newVertex);
boundaryPathIterator.next();
System.arraycopy(newVertex, 0, oldVertex, 0, 2);
segmentType = boundaryPathIterator.currentSegment(newVertex);
}
processSegment(newVertex, firstVertex);
}
Давайте на секунду взглянем на initializeForNewCirclePrivate. Этот метод просто устанавливает поля экземпляра и позволяет производному классу хранить любое свойство круга. Его определение
private void initializeForNewCirclePrivate(Circle circle) {
currentCircle = circle;
currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius();
initializeForNewCircle(circle);
}
initializeForNewCircle является абстрактным, и одной из реализаций было бы сохранение радиуса окружностей во избежание появления квадратных корней. В любом случае вернемся к processCircleShape. После вызова initializeForNewCirclePrivate мы проверяем, является ли полигон null (который я интерпретирую как пустой многоугольник), и возвращаем, если это null. В этом случае наша вычисленная площадь будет равна нулю. Если многоугольник не null, то мы получаем PathIterator многоугольника. Аргумент к методу getPathIterator, который я вызываю, является аффинным преобразованием, которое может быть применено к пути. Я не хочу применять один, хотя, поэтому я просто передать null.
Затем я объявляю double[] s, которые будут отслеживать вершины. Я должен помнить первую вершину, потому что PathIterator дает мне каждую вершину только один раз, поэтому я должен вернуться назад после того, как она дала мне последнюю вершину, и сформировать ребро с этой последней и первой вершиной.
Метод currentSegment в следующей строке помещает следующую вершину в свой аргумент. Он возвращает код, который сообщает вам, когда он находится вне вершин. Вот почему контрольное выражение для моего цикла while таково.
Большая часть остальной части кода этого метода представляет собой неинтересную логику, связанную с итерацией по вершинам. Важно то, что один раз за итерацию цикла while я вызываю processSegment, а затем снова вызываю processSegment в конце метода, чтобы обработать ребро, соединяющее последнюю вершину с первой вершиной.
Давайте посмотрим на код для processSegment:
private void processSegment(double[] initialVertex, double[] finalVertex) {
double[] segmentDisplacement = displacment2D(initialVertex, finalVertex);
if (segmentDisplacement[0] == 0 && segmentDisplacement[1] == 0) {
return;
}
double segmentLength = Math.sqrt(dotProduct2D(segmentDisplacement, segmentDisplacement));
double[] centerToInitialDisplacement = new double[]{initialVertex[0] - getCurrentCircle().getCenterX(), initialVertex[1] - getCurrentCircle().getCenterY()};
final double leftX = dotProduct2D(centerToInitialDisplacement, segmentDisplacement) / segmentLength;
final double rightX = leftX + segmentLength;
final double y = wedgeProduct2D(segmentDisplacement, centerToInitialDisplacement) / segmentLength;
processSegmentStandardGeometry(leftX, rightX, y);
}
В этом методе я реализую шаги по преобразованию ребра в стандартную геометрию, как описано выше. Сначала я вычисляю segmentDisplacement, смещение от начальной вершины до конечной вершины. Это определяет ось х стандартной геометрии. Я делаю досрочный возврат, если это смещение равно нулю.
Далее я вычисляю длину смещения, потому что это необходимо для получения вектора единицы x. Получив эту информацию, я вычисляю смещение от центра круга к начальной вершине. Точечный продукт этого с segmentDisplacement дает мне leftX, который я называл xi. Тогда rightX, который я называл xf, это просто leftX + segmentLength. Наконец я делаю клиновое произведение, чтобы получить y, как описано выше.
Теперь, когда я преобразовал задачу в стандартную геометрию, с ней будет легко иметь дело. Это то, что делает метод processSegmentStandardGeometry. Давайте посмотрим на код
private void processSegmentStandardGeometry(double leftX, double rightX, double y) {
if (y * y > getCurrentSquareRadius()) {
processNonIntersectingRegion(leftX, rightX, y);
} else {
final double intersectionX = Math.sqrt(getCurrentSquareRadius() - y * y);
if (leftX < -intersectionX) {
final double leftRegionRightEndpoint = Math.min(-intersectionX, rightX);
processNonIntersectingRegion(leftX, leftRegionRightEndpoint, y);
}
if (intersectionX < rightX) {
final double rightRegionLeftEndpoint = Math.max(intersectionX, leftX);
processNonIntersectingRegion(rightRegionLeftEndpoint, rightX, y);
}
final double middleRegionLeftEndpoint = Math.max(-intersectionX, leftX);
final double middleRegionRightEndpoint = Math.min(intersectionX, rightX);
final double middleRegionLength = Math.max(middleRegionRightEndpoint - middleRegionLeftEndpoint, 0);
processIntersectingRegion(middleRegionLength, y);
}
}
Первый if различает случаи, когда y достаточно мал, чтобы ребро могло пересекать окружность. Если y большое и нет возможности пересечения, то я вызываю метод для обработки этого случая. В противном случае я рассматриваю случай, когда возможно пересечение.
Если пересечение возможно, я вычисляю координату x пересечения, intersectionX, и делю ребро на три части, которые соответствуют областям 1, 2 и 3 стандартного геометрического рисунка выше. Сначала я работаю с регионом 1.
Чтобы обработать область 1, я проверяю, действительно ли leftX меньше -intersectionX, иначе в противном случае не было бы области 1. Если существует область 1, то мне нужно знать, когда она заканчивается. Это заканчивается как минимум rightX и -intersectionX. После того как я нашел эти x-координаты, я имею дело с этой непересекающейся областью
Я делаю то же самое для обработки области 3.
Для региона 2 я должен сделать некоторую логику, чтобы проверить, что leftX и rightX действительно заключают в скобки некоторую область между -intersectionX и intersectionX. После нахождения региона мне нужна только длина региона и y, поэтому я передаю эти два числа абстрактному методу, который обрабатывает регион 2.
Теперь давайте посмотрим на код для processNonIntersectingRegion
private void processNonIntersectingRegion(double leftX, double rightX, double y) {
final double initialTheta = Math.atan2(y, leftX);
final double finalTheta = Math.atan2(y, rightX);
double deltaTheta = finalTheta - initialTheta;
if (deltaTheta < -Math.PI) {
deltaTheta += 2 * Math.PI;
} else if (deltaTheta > Math.PI) {
deltaTheta -= 2 * Math.PI;
}
processNonIntersectingRegion(deltaTheta);
}
Я просто использую atan2, чтобы вычислить разницу в углах между leftX и rightX. Затем я добавляю код для разрешения разрыва в atan2, но это, вероятно, не нужно, потому что разрыв происходит либо при 180 градусах, либо при 0 градусах. Затем я передаю разницу в углах на абстрактный метод. Наконец, у нас есть только абстрактные методы и методы получения:
protected abstract void initialize();
protected abstract void initializeForNewCircle(Circle circle);
protected abstract void processNonIntersectingRegion(double deltaTheta);
protected abstract void processIntersectingRegion(double length, double y);
protected abstract T getValue();
protected final Circle getCurrentCircle() {
return currentCircle;
}
protected final double getCurrentSquareRadius() {
return currentSquareRadius;
}
}
Теперь давайте посмотрим на класс расширения, CircleAreaFinder
public class CircleAreaFinder extends CircleShapeIntersectionFinder<Double> {
public static double findAreaOfCircle(Circle circle, Shape shape) {
CircleAreaFinder circleAreaFinder = new CircleAreaFinder();
return circleAreaFinder.computeValue(circle, shape);
}
double area;
@Override
protected void initialize() {
area = 0;
}
@Override
protected void processNonIntersectingRegion(double deltaTheta) {
area += getCurrentSquareRadius() * deltaTheta / 2;
}
@Override
protected void processIntersectingRegion(double length, double y) {
area -= length * y / 2;
}
@Override
protected Double getValue() {
return area;
}
@Override
protected void initializeForNewCircle(Circle circle) {
}
}
Имеется поле area для отслеживания области. initialize устанавливает область на ноль, как и ожидалось. Когда мы обрабатываем непересекающийся край, мы увеличиваем площадь на R ^ 2 Δθ / 2, как мы и пришли к выводу, что мы должны выше. Для пересекающегося края мы уменьшаем площадь на y*length/2. Это было так, что отрицательные значения для y соответствуют положительным областям, как мы решили, что они должны.
Теперь, если мы хотим отслеживать периметр, нам не нужно делать намного больше работы. Я определил AreaPerimeter класс:
public class AreaPerimeter {
final double area;
final double perimeter;
public AreaPerimeter(double area, double perimeter) {
this.area = area;
this.perimeter = perimeter;
}
public double getArea() {
return area;
}
public double getPerimeter() {
return perimeter;
}
}
и теперь нам просто нужно снова расширить наш абстрактный класс, используя AreaPerimeter в качестве типа.
public class CircleAreaPerimeterFinder extends CircleShapeIntersectionFinder<AreaPerimeter> {
public static AreaPerimeter findAreaPerimeterOfCircle(Circle circle, Shape shape) {
CircleAreaPerimeterFinder circleAreaPerimeterFinder = new CircleAreaPerimeterFinder();
return circleAreaPerimeterFinder.computeValue(circle, shape);
}
double perimeter;
double radius;
CircleAreaFinder circleAreaFinder;
@Override
protected void initialize() {
perimeter = 0;
circleAreaFinder = new CircleAreaFinder();
}
@Override
protected void initializeForNewCircle(Circle circle) {
radius = Math.sqrt(getCurrentSquareRadius());
}
@Override
protected void processNonIntersectingRegion(double deltaTheta) {
perimeter += deltaTheta * radius;
circleAreaFinder.processNonIntersectingRegion(deltaTheta);
}
@Override
protected void processIntersectingRegion(double length, double y) {
perimeter += Math.abs(length);
circleAreaFinder.processIntersectingRegion(length, y);
}
@Override
protected AreaPerimeter getValue() {
return new AreaPerimeter(circleAreaFinder.getValue(), perimeter);
}
}
У нас есть переменная perimeter для отслеживания периметра, мы помним значение radius, чтобы избежать необходимости многократно вызывать Math.sqrt, и мы делегируем вычисление площади нашему CircleAreaFinder , Мы видим, что формулы для периметра просты.
Для справки вот полный код CircleShapeIntersectionFinder
private static double[] displacment2D(final double[] initialPoint, final double[] finalPoint) {
return new double[]{finalPoint[0] - initialPoint[0], finalPoint[1] - initialPoint[1]};
}
private static double wedgeProduct2D(final double[] firstFactor, final double[] secondFactor) {
return firstFactor[0] * secondFactor[1] - firstFactor[1] * secondFactor[0];
}
static private double dotProduct2D(final double[] firstFactor, final double[] secondFactor) {
return firstFactor[0] * secondFactor[0] + firstFactor[1] * secondFactor[1];
}
private Circle currentCircle;
private double currentSquareRadius;
public final T computeValue(Circle circle, Shape shape) {
initialize();
processCircleShape(circle, shape);
return getValue();
}
private void processCircleShape(Circle circle, final Shape cellBoundaryPolygon) {
initializeForNewCirclePrivate(circle);
if (cellBoundaryPolygon == null) {
return;
}
PathIterator boundaryPathIterator = cellBoundaryPolygon.getPathIterator(null);
double[] firstVertex = new double[2];
double[] oldVertex = new double[2];
double[] newVertex = new double[2];
int segmentType = boundaryPathIterator.currentSegment(firstVertex);
if (segmentType != PathIterator.SEG_MOVETO) {
throw new AssertionError();
}
System.arraycopy(firstVertex, 0, newVertex, 0, 2);
boundaryPathIterator.next();
System.arraycopy(newVertex, 0, oldVertex, 0, 2);
segmentType = boundaryPathIterator.currentSegment(newVertex);
while (segmentType != PathIterator.SEG_CLOSE) {
processSegment(oldVertex, newVertex);
boundaryPathIterator.next();
System.arraycopy(newVertex, 0, oldVertex, 0, 2);
segmentType = boundaryPathIterator.currentSegment(newVertex);
}
processSegment(newVertex, firstVertex);
}
private void initializeForNewCirclePrivate(Circle circle) {
currentCircle = circle;
currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius();
initializeForNewCircle(circle);
}
private void processSegment(double[] initialVertex, double[] finalVertex) {
double[] segmentDisplacement = displacment2D(initialVertex, finalVertex);
if (segmentDisplacement[0] == 0 && segmentDisplacement[1] == 0) {
return;
}
double segmentLength = Math.sqrt(dotProduct2D(segmentDisplacement, segmentDisplacement));
double[] centerToInitialDisplacement = new double[]{initialVertex[0] - getCurrentCircle().getCenterX(), initialVertex[1] - getCurrentCircle().getCenterY()};
final double leftX = dotProduct2D(centerToInitialDisplacement, segmentDisplacement) / segmentLength;
final double rightX = leftX + segmentLength;
final double y = wedgeProduct2D(segmentDisplacement, centerToInitialDisplacement) / segmentLength;
processSegmentStandardGeometry(leftX, rightX, y);
}
private void processSegmentStandardGeometry(double leftX, double rightX, double y) {
if (y * y > getCurrentSquareRadius()) {
processNonIntersectingRegion(leftX, rightX, y);
} else {
final double intersectionX = Math.sqrt(getCurrentSquareRadius() - y * y);
if (leftX < -intersectionX) {
final double leftRegionRightEndpoint = Math.min(-intersectionX, rightX);
processNonIntersectingRegion(leftX, leftRegionRightEndpoint, y);
}
if (intersectionX < rightX) {
final double rightRegionLeftEndpoint = Math.max(intersectionX, leftX);
processNonIntersectingRegion(rightRegionLeftEndpoint, rightX, y);
}
final double middleRegionLeftEndpoint = Math.max(-intersectionX, leftX);
final double middleRegionRightEndpoint = Math.min(intersectionX, rightX);
final double middleRegionLength = Math.max(middleRegionRightEndpoint - middleRegionLeftEndpoint, 0);
processIntersectingRegion(middleRegionLength, y);
}
}
private void processNonIntersectingRegion(double leftX, double rightX, double y) {
final double initialTheta = Math.atan2(y, leftX);
final double finalTheta = Math.atan2(y, rightX);
double deltaTheta = finalTheta - initialTheta;
if (deltaTheta < -Math.PI) {
deltaTheta += 2 * Math.PI;
} else if (deltaTheta > Math.PI) {
deltaTheta -= 2 * Math.PI;
}
processNonIntersectingRegion(deltaTheta);
}
protected abstract void initialize();
protected abstract void initializeForNewCircle(Circle circle);
protected abstract void processNonIntersectingRegion(double deltaTheta);
protected abstract void processIntersectingRegion(double length, double y);
protected abstract T getValue();
protected final Circle getCurrentCircle() {
return currentCircle;
}
protected final double getCurrentSquareRadius() {
return currentSquareRadius;
}
Во всяком случае, это мое описание алгоритма. Я думаю, что это хорошо, потому что это точно, и не так много случаев, чтобы проверить.