首先,我将提醒我们如何找到多边形的面积。完成此操作后,查找多边形和圆的交点的算法应该很容易理解。
如何求多边形的面积
让我们看一下三角形的情况,因为所有基本逻辑都出现在那里。假设我们有一个顶点为 (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,请阅读本文。
第二项是这个领域,这也是积极的:
第三个区域如下图所示。这次面积为负
将这三个加起来我们得到下面的图片
我们看到三角形之外的绿色区域被红色区域抵消了,因此净面积就是三角形的面积,这说明了为什么我们的公式在这种情况下是正确的。
我上面说的是对面积公式为什么正确的直观解释。更严格的解释是观察当我们从边缘计算面积时,我们得到的面积与我们从积分 r^2dθ/2 得到的面积相同,因此我们有效地在边界周围积分 r^2dθ/2根据斯托克斯定理,这给出了与在多边形边界区域上积分 rdrdθ 相同的结果。由于在由多边形包围的区域上积分 rdrdθ 给出了面积,我们得出结论,我们的程序必须正确给出面积。
圆与多边形相交的面积
现在让我们讨论如何找到半径为 R 的圆与多边形的交点面积,如下图所示:
我们有兴趣找到绿色区域的面积。就像在单个多边形的情况下一样,我们可以将计算分解为为多边形的每一边找到一个区域,然后将这些区域相加。
我们的第一个区域将如下所示:
第二个区域看起来像
第三个领域将是
同样,在我们的案例中,前两个领域是积极的,而第三个领域是消极的。希望取消能够解决,以便净面积确实是我们感兴趣的区域。让我们看看。
事实上,这些区域的总和将是我们感兴趣的区域。
同样,我们可以更严格地解释为什么会这样。令 I 为交点定义的区域,令 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 相同。
现在做一个变量的改变,我们发现如果我们在 π(P) 的边界上积分 r^2 dθ/2 会得到相同的答案,其中 π(P) 是 P 在 π 下的图像。
再次使用斯托克斯定理,我们知道在 π(P) 的边界上积分 r^2 dθ/2 得到了 π(P) 的面积。换句话说,它给出了与在 π(P) 上积分 dxdy 相同的答案。
再次使用变量的变化,我们知道 dxdy 在 π(P) 上的积分与 Jdxdy 在 P 上的积分相同,其中 J 是 π 的雅可比。
现在我们可以将 Jdxdy 的积分分成两个区域:圆内部分和圆外部分。现在 π 只在圆中留下点,所以在那里 J=1,所以 P 的这一部分的贡献是 P 位于圆中的部分的面积,即交点的面积。第二个区域是圆外的区域。J=0,因为 π 将这部分折叠到圆的边界。
因此,我们计算的确实是交叉点的面积。
现在我们相对确定我们在概念上知道如何找到该区域,让我们更具体地讨论如何计算单个段的贡献。让我们从我称之为“标准几何”的部分开始。如下所示。
在标准几何中,边缘从左到右水平移动。它由三个数字描述:xi,边缘开始的 x 坐标,xf,边缘结束的 x 坐标,以及 y,边缘的 y 坐标。
现在我们看到如果 |y| < R,如图所示,那么边将在点 (-xint,y) 和 (xint,y) 处与圆相交,其中 xint = (R^2-y^2)^(1/2)。然后我们需要计算的区域被分成图中标记的三个部分。要获得区域 1 和 3 的面积,我们可以使用 arctan 来获得各个点的角度,然后将面积等同于 R^2 Δθ/2。例如,我们将设置 θi = atan2(y,xi) 和 θl = atan2(y,-xint)。那么区域一的面积是R^2 (θl-θi)/2。类似地,我们可以得到区域 3 的面积。
区域 2 的面积就是三角形的面积。但是,我们必须小心标志。我们希望显示的面积是正的,所以我们会说面积是 -(xint - (-xint))y/2。
要记住的另一件事是,通常,xi 不必小于 -xint,xf 不必大于 xint。
另一个要考虑的情况是 |y| > R. 这种情况比较简单,因为只有一块与图中区域1相似。
既然我们知道如何从标准几何中的边计算面积,剩下要做的就是描述如何将任何边转换为标准几何。
但这只是坐标的简单变化。给定一些初始顶点 vi 和最终顶点 vf,新的 x 单位向量将是从 vi 指向 vf 的单位向量。那么 xi 就是 vi 从点入 x 的圆心的位移,而 xf 就是 xi 加上 vi 和 vf 之间的距离。同时 y 由 x 与 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;
}
}
对于多边形,我将使用 java 的Shape
类。Shape
s 有一个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];
}
有两个实例字段,aCircle
只保留圆的副本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 轴。如果这个位移为零,我会提前返回。
接下来我计算位移的长度,因为这是获得 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 度。然后我将角度差异传递给抽象方法。最后,我们只有抽象方法和 getter:
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;
}
无论如何,这就是我对算法的描述。我认为这很好,因为它是准确的,并且没有太多需要检查的案例。