如何计算三角形(指定为三(X,Y)对)和圆(X,Y,R)之间的交叉区域？我做了一些搜索无济于事.这是为了工作,而不是学校.:)
它在C#中看起来像这样:
struct { PointF vert[3]; } Triangle;
struct { PointF center; float radius; } Circle;
// returns the area of intersection, e.g.:
// if the circle contains the triangle, return area of triangle
// if the triangle contains the circle, return area of circle
// if partial intersection, figure that out
// if no intersection, return 0
double AreaOfIntersection(Triangle t, Circle c)
{
...
}
Brian Moths.. 32
首先,我将提醒我们如何找到多边形的区域.完成此操作后,找到多边形和圆形之间交点的算法应该很容易理解.
如何找到多边形的区域让我们看一下三角形的情况,因为那里出现了所有必要的逻辑.假设我们在逆时针绕三角形时有一个带顶点(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的圆与多边形的交点区域,如下图所示:
我们有兴趣找到绿色区域.正如在单个多边形的情况下,我们可以将计算分解为为多边形的每一边找到一个区域,然后向上添加这些区域.
我们的第一个区域将如下所示:
第二个区域看起来像
第三个区域将是
同样,前两个区域在我们的情况下是积极的,而第三个区域将是负面的.希望取消将有效,以便净面积确实是我们感兴趣的区域.让我们看看.
实际上,这些区域的总和将是我们感兴趣的区域.
同样,我们可以更加严格地解释其原因.让我成为交集定义的区域,让P为多边形.然后从前面的讨论中,我们知道我们想要计算出围绕I边界的r ^ 2d？/ 2的积分.但是,这很难做到因为它需要找到交点.
相反,我们在多边形上做了一个整数.我们在多边形的边界上积分了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)的面积.换句话说,它给出了与将dxdy整合在一起的相同答案？(P).
再次使用变量变量,我们知道将dxdy积分超过？(P)与将Jdxdy积分在P上相同,其中J是雅各比的？
现在我们可以将Jdxdy的积分分成两个区域:圆圈中的部分和圆圈外的部分.现在？只留下圆圈中的点,所以J = 1,所以P的这一部分的贡献是位于圆圈中的P部分的面积,即交点的面积.第二个区域是圆圈外的区域.那J = 0？将此部分折叠到圆的边界.
因此,我们计算的确实是交叉区域.
现在我们相对肯定我们在概念上知道如何找到该区域,让我们更具体地讨论如何计算单个区段的贡献.让我们首先看一下我将称之为"标准几何"的片段.如下所示.
In standard geometry, the edge goes horizontally from left to right. It is described by three numbers: xi, the x-coordinate where the edge starts, xf, the x-coordinate where the edge ends, and y, the y coordinate of the edge.
Now we see that if |y| < R, as in the figure, then the edge will intersect the circle at the points (-xint,y) and (xint,y) where xint = (R^2-y^2)^(1/2). Then the area we need to calculate is broken up into three pieces labelled in the figure. To get the areas of regions 1 and 3, we can use arctan to get the angles of the various points and then equate the area to R^2 ??/2. So for example we would set ?i = atan2(y,xi) and ?l = atan2(y,-xint). Then the area of region one is R^2 (?l-?i)/2. We can obtain the area of region 3 similarly.
The area of region 2 is just the area of a triangle. However, we must be careful about sign. We want the area shown to be positive so we will say the area is -(xint - (-xint))y/2.
Another thing to keep in mind is that in general, xi does not have to be less than -xint and xf does not have to be greater than xint.
The other case to consider is |y| > R. This case is simpler, because there is only one piece which is similar to region 1 in the figure.
Now that we know how to compute the area from an edge in standard geometry, the only thing left to do is describe how to transform any edge into standard geometry.
But this just a simple change of coordinates. Given some with initial vertex vi and final vertex vf, the new x unit vector will be the unit vector pointing from vi to vf. Then xi is just the displacement of vi from the center of the circle dotted into x, and xf is just xi plus the distance between vi and vf. Meanwhile y is given by the wedge product of x with the displacement of vi from the center of the circle.
CodeThat completes the description of the algorithm, now it is time to write some code. I will use java.
First off, since we are working with circles, we should have a circle class
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;
}
}
For polygons, I will use java's Shape class. Shapes have a PathIterator that I can use to iterate through the edges of the polygon.
Now for the actual work. I will separate the logic of iterating through the edges, putting the edges in standard geometry etc, from the logic of computing the area once this is done. The reason for this is that you may in the future want to compute something else besides or in addition to the area and you want to be able to reuse the code having to deal with iterating through the edges.
So I have a generic class which computes some property of class T about our polygon circle intersection.
public abstract class CircleShapeIntersectionFinder
It has three static methods that just help compute geometry:
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];
}
There are two instance fields, a Circle which just keeps a copy of the circle, and the currentSquareRadius, which keeps a copy of the square radius. This may seem odd, but the class I am using is actually equipped to find the areas of a whole collection of circle-polygon intersections. That is why I am referring to one of the circles as "current".
private Circle currentCircle; private double currentSquareRadius;
Next comes the method for computing what we want to compute:
public final T computeValue(Circle circle, Shape shape) {
initialize();
processCircleShape(circle, shape);
return getValue();
}
initialize() and getValue() are abstract. initialize() would set the variable that is keeping a total of the area to zero, and getValue() would just return the area. The definition for processCircleShape is
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);
}
Let's take a second to look at initializeForNewCirclePrivate quickly. This method just sets the instance fields and allows the derived class to store any property of the circle. Its definition is
private void initializeForNewCirclePrivate(Circle circle) {
currentCircle = circle;
currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius();
initializeForNewCircle(circle);
}
initializeForNewCircle is abstract and one implementation would be for it to store the circles radius to avoid having to do square roots. Anyway back to processCircleShape. After calling initializeForNewCirclePrivate, we check if the polygon is null (which I am interpreting as an empty polygon), and we return if it is null. In this case, our computed area would be zero. If the polygon is not null then we get the PathIterator of the polygon. The argument to the getPathIterator method I call is an affine transformation that can be applied to the path. I don't want to apply one though, so I just pass null.
Next I declare the double[]s that will keep track of the vertices. I must remember the first vertex because the PathIterator only gives me each vertex once, so I have to go back after it has given me the last vertex, and form an edge with this last vertex and the first vertex.
The currentSegment method on the next line puts the next vertex in its argument. It returns a code that tells you when it is out of vertices. This is why the control expression for my while loop is what it is.
Most of the rest of the code of this method is uninteresting logic related to iterating through the vertices. The important thing is that once per iteration of the while loop I call processSegment and then I call processSegment again at the end of the method to process the edge that connects the last vertex to the first vertex.
Let's look at the code for 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);
}
In this method I implement the steps to transform an edge into the standard geometry as described above. First I calculate segmentDisplacement, the displacement from the initial vertex to the final vertex. This defines the x axis of the standard geometry. I do an early return if this displacement is zero.
Next I calculate the length of the displacement, because this is necessary to get the x unit vector. Once I have this information, I calculate the displacement from the center of the circle to the initial vertex. The dot product of this with segmentDisplacement gives me leftX which I had been calling xi. Then rightX, which I had been calling xf, is just leftX + segmentLength. Finally I do the wedge product to get y as described above.
Now that I have transformed the problem into the standard geometry, it will be easy to deal with. That is what the processSegmentStandardGeometry method does. Let's look at the code
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);
}
}
The first if distinguishes the cases where y is small enough that the edge may intersect the circle. If y is big and there is no possibility of intersection, then I call the method to handle that case. Otherwise I handle the case where intersection is possible.
If intersection is possible, I calculate the x coordinate of intersection, intersectionX, and I divide the edge up into three portions, which correspond to regions 1, 2, and 3 of the standard geometry figure above. First I handle region 1.
To handle region 1, I check if leftX is indeed less than -intersectionX for otherwise there would be no region 1. If there is a region 1, then I need to know when it ends. It ends at the minimum of rightX and -intersectionX. After I have found these x-coordinates, I deal with this non-intersection region.
I do a similar thing to handle region 3.
For region 2, I have to do some logic to check that leftX and rightX do actually bracket some region in between -intersectionX and intersectionX. After finding the region, I only need the length of the region and y, so I pass these two numbers on to an abstract method which handles the region 2.
Now let's look at the code for 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);
}
I simply use atan2 to calculate the difference in angle between leftX and rightX. Then I add code to deal with the discontinuity in atan2, but this is probably unnecessary, because the discontinuity occurs either at 180 degrees or 0 degrees. Then I pass the difference in angle onto an abstract method. Lastly we just have abstract methods and getters:
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;
}
}
Now let's look at the extending class, CircleAreaFinder
public class CircleAreaFinder extends CircleShapeIntersectionFinder{ 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) { }
}
It has a field area to keep track of the area. initialize sets area to zero, as expected. When we process a non intersecting edge, we increment the area by R^2 ??/2 as we concluded we should above. For an intersecting edge, we decrement the area by y*length/2. This was so that negative values for y correspond to positive areas, as we decided they should.
现在好的一点是,如果我们想要跟踪周边,我们就没有必要做更多的工作.我定义了一个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{ 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;
}
无论如何,这是我对算法的描述.我认为这很好,因为它是确切的,并没有那么多的情况要检查.
首先,我将提醒我们如何找到多边形的区域.完成此操作后,找到多边形和圆形之间交点的算法应该很容易理解.
如何找到多边形的区域让我们看一下三角形的情况,因为那里出现了所有必要的逻辑.假设我们在逆时针绕三角形时有一个带顶点(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的圆与多边形的交点区域,如下图所示:
我们有兴趣找到绿色区域.正如在单个多边形的情况下,我们可以将计算分解为为多边形的每一边找到一个区域,然后向上添加这些区域.
我们的第一个区域将如下所示:
第二个区域看起来像
第三个区域将是
同样,前两个区域在我们的情况下是积极的,而第三个区域将是负面的.希望取消将有效,以便净面积确实是我们感兴趣的区域.让我们看看.
实际上,这些区域的总和将是我们感兴趣的区域.
同样,我们可以更加严格地解释其原因.让我成为交集定义的区域,让P为多边形.然后从前面的讨论中,我们知道我们想要计算出围绕I边界的r ^ 2d？/ 2的积分.但是,这很难做到因为它需要找到交点.
相反,我们在多边形上做了一个整数.我们在多边形的边界上积分了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)的面积.换句话说,它给出了与将dxdy整合在一起的相同答案？(P).
再次使用变量变量,我们知道将dxdy积分超过？(P)与将Jdxdy积分在P上相同,其中J是雅各比的？
现在我们可以将Jdxdy的积分分成两个区域:圆圈中的部分和圆圈外的部分.现在？只留下圆圈中的点,所以J = 1,所以P的这一部分的贡献是位于圆圈中的P部分的面积,即交点的面积.第二个区域是圆圈外的区域.那J = 0？将此部分折叠到圆的边界.
因此,我们计算的确实是交叉区域.
现在我们相对肯定我们在概念上知道如何找到该区域,让我们更具体地讨论如何计算单个区段的贡献.让我们首先看一下我将称之为"标准几何"的片段.如下所示.
In standard geometry, the edge goes horizontally from left to right. It is described by three numbers: xi, the x-coordinate where the edge starts, xf, the x-coordinate where the edge ends, and y, the y coordinate of the edge.
Now we see that if |y| < R, as in the figure, then the edge will intersect the circle at the points (-xint,y) and (xint,y) where xint = (R^2-y^2)^(1/2). Then the area we need to calculate is broken up into three pieces labelled in the figure. To get the areas of regions 1 and 3, we can use arctan to get the angles of the various points and then equate the area to R^2 ??/2. So for example we would set ?i = atan2(y,xi) and ?l = atan2(y,-xint). Then the area of region one is R^2 (?l-?i)/2. We can obtain the area of region 3 similarly.
The area of region 2 is just the area of a triangle. However, we must be careful about sign. We want the area shown to be positive so we will say the area is -(xint - (-xint))y/2.
Another thing to keep in mind is that in general, xi does not have to be less than -xint and xf does not have to be greater than xint.
The other case to consider is |y| > R. This case is simpler, because there is only one piece which is similar to region 1 in the figure.
Now that we know how to compute the area from an edge in standard geometry, the only thing left to do is describe how to transform any edge into standard geometry.
But this just a simple change of coordinates. Given some with initial vertex vi and final vertex vf, the new x unit vector will be the unit vector pointing from vi to vf. Then xi is just the displacement of vi from the center of the circle dotted into x, and xf is just xi plus the distance between vi and vf. Meanwhile y is given by the wedge product of x with the displacement of vi from the center of the circle.
CodeThat completes the description of the algorithm, now it is time to write some code. I will use java.
First off, since we are working with circles, we should have a circle class
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;
}
}
For polygons, I will use java's Shape class. Shapes have a PathIterator that I can use to iterate through the edges of the polygon.
Now for the actual work. I will separate the logic of iterating through the edges, putting the edges in standard geometry etc, from the logic of computing the area once this is done. The reason for this is that you may in the future want to compute something else besides or in addition to the area and you want to be able to reuse the code having to deal with iterating through the edges.
So I have a generic class which computes some property of class T about our polygon circle intersection.
public abstract class CircleShapeIntersectionFinder
It has three static methods that just help compute geometry:
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];
}
There are two instance fields, a Circle which just keeps a copy of the circle, and the currentSquareRadius, which keeps a copy of the square radius. This may seem odd, but the class I am using is actually equipped to find the areas of a whole collection of circle-polygon intersections. That is why I am referring to one of the circles as "current".
private Circle currentCircle; private double currentSquareRadius;
Next comes the method for computing what we want to compute:
public final T computeValue(Circle circle, Shape shape) {
initialize();
processCircleShape(circle, shape);
return getValue();
}
initialize() and getValue() are abstract. initialize() would set the variable that is keeping a total of the area to zero, and getValue() would just return the area. The definition for processCircleShape is
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);
}
Let's take a second to look at initializeForNewCirclePrivate quickly. This method just sets the instance fields and allows the derived class to store any property of the circle. Its definition is
private void initializeForNewCirclePrivate(Circle circle) {
currentCircle = circle;
currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius();
initializeForNewCircle(circle);
}
initializeForNewCircle is abstract and one implementation would be for it to store the circles radius to avoid having to do square roots. Anyway back to processCircleShape. After calling initializeForNewCirclePrivate, we check if the polygon is null (which I am interpreting as an empty polygon), and we return if it is null. In this case, our computed area would be zero. If the polygon is not null then we get the PathIterator of the polygon. The argument to the getPathIterator method I call is an affine transformation that can be applied to the path. I don't want to apply one though, so I just pass null.
Next I declare the double[]s that will keep track of the vertices. I must remember the first vertex because the PathIterator only gives me each vertex once, so I have to go back after it has given me the last vertex, and form an edge with this last vertex and the first vertex.
The currentSegment method on the next line puts the next vertex in its argument. It returns a code that tells you when it is out of vertices. This is why the control expression for my while loop is what it is.
Most of the rest of the code of this method is uninteresting logic related to iterating through the vertices. The important thing is that once per iteration of the while loop I call processSegment and then I call processSegment again at the end of the method to process the edge that connects the last vertex to the first vertex.
Let's look at the code for 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);
}
In this method I implement the steps to transform an edge into the standard geometry as described above. First I calculate segmentDisplacement, the displacement from the initial vertex to the final vertex. This defines the x axis of the standard geometry. I do an early return if this displacement is zero.
Next I calculate the length of the displacement, because this is necessary to get the x unit vector. Once I have this information, I calculate the displacement from the center of the circle to the initial vertex. The dot product of this with segmentDisplacement gives me leftX which I had been calling xi. Then rightX, which I had been calling xf, is just leftX + segmentLength. Finally I do the wedge product to get y as described above.
Now that I have transformed the problem into the standard geometry, it will be easy to deal with. That is what the processSegmentStandardGeometry method does. Let's look at the code
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);
}
}
The first if distinguishes the cases where y is small enough that the edge may intersect the circle. If y is big and there is no possibility of intersection, then I call the method to handle that case. Otherwise I handle the case where intersection is possible.
If intersection is possible, I calculate the x coordinate of intersection, intersectionX, and I divide the edge up into three portions, which correspond to regions 1, 2, and 3 of the standard geometry figure above. First I handle region 1.
To handle region 1, I check if leftX is indeed less than -intersectionX for otherwise there would be no region 1. If there is a region 1, then I need to know when it ends. It ends at the minimum of rightX and -intersectionX. After I have found these x-coordinates, I deal with this non-intersection region.
I do a similar thing to handle region 3.
For region 2, I have to do some logic to check that leftX and rightX do actually bracket some region in between -intersectionX and intersectionX. After finding the region, I only need the length of the region and y, so I pass these two numbers on to an abstract method which handles the region 2.
Now let's look at the code for 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);
}
I simply use atan2 to calculate the difference in angle between leftX and rightX. Then I add code to deal with the discontinuity in atan2, but this is probably unnecessary, because the discontinuity occurs either at 180 degrees or 0 degrees. Then I pass the difference in angle onto an abstract method. Lastly we just have abstract methods and getters:
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;
}
}
Now let's look at the extending class, CircleAreaFinder
public class CircleAreaFinder extends CircleShapeIntersectionFinder{ 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) { }
}
It has a field area to keep track of the area. initialize sets area to zero, as expected. When we process a non intersecting edge, we increment the area by R^2 ??/2 as we concluded we should above. For an intersecting edge, we decrement the area by y*length/2. This was so that negative values for y correspond to positive areas, as we decided they should.
现在好的一点是,如果我们想要跟踪周边,我们就没有必要做更多的工作.我定义了一个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{ 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;
}
无论如何,这是我对算法的描述.我认为这很好,因为它是确切的,并没有那么多的情况要检查.
如果你想要一个精确的解决方案(或者至少与使用浮点运算一样精确),那么这将涉及大量的工作,因为有很多情况需要考虑.
我计算了九种不同的情况(在下图中按圆圈内三角形的顶点数量分类,以及圆圈中相交或包含的三角形边缘数量):
(然而,众所周知,这种几何情况的枚举是棘手的,如果我错过了一两个,它就不会让我感到惊讶!)
所以方法是:
确定三角形的每个顶点是否在圆内.我假设你知道怎么做.
确定三角形的每个边缘是否与圆相交.(我在这里写了一个方法,或者看到任何计算几何图书.)你需要计算在步骤4中使用的一个或多个交点(如果有的话).
确定您拥有的九个案例中的哪一个.
计算交叉口的面积.案例1,2和9很容易.在其余六种情况下,我绘制了虚线,以显示如何根据三角形的原始顶点以及在步骤2中计算的交点将交叉区域划分为三角形和圆形线段.
这个算法会相当精细并且容易出错,只会影响其中一个案例,因此请确保您拥有涵盖所有九个案例的测试用例(我建议也可以置换测试三角形的顶点).特别注意三角形的一个顶点位于圆的边缘的情况.
如果你不需要一个精确的解决方案,那么光栅化数字和计算交集中的像素(正如其他几个受访者所建议的那样)似乎是一种更容易的代码方法,相应地也不容易出错.
京公网安备 11010802040832号 | 京ICP备19059560号-6 