まず、多角形の面積を求める方法を思い出してください。これが完了すると、多角形と円の交点を見つけるアルゴリズムは理解しやすくなります。
多角形の面積を見つける方法
三角形の場合を見てみましょう。基本的なロジックはすべてそこにあるからです。次の図に示すように、頂点が (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 であることが明らかでない場合は、これを読んでください。
2 番目の項はこの領域で、これもまた正です。
そして、3 番目の領域を次の図に示します。今回は面積がマイナス
これら3つを追加すると、次の図が得られます
三角形の外側にあった緑色の領域が赤色の領域によってキャンセルされるため、正味の領域は三角形の領域にすぎないことがわかります。これは、この場合に式が真であった理由を示しています.
上で述べたことは、なぜ面積式が正しいのかについての直感的な説明でした。より厳密な説明は、エッジから面積を計算するときに得られる面積は、積分 r^2dθ/2 から得られる面積と同じであるため、境界の周りで r^2dθ/2 を効果的に積分していることを観察することです。であり、ストークスの定理により、これは、多角形に囲まれた領域で rdrdθ を積分した場合と同じ結果をもたらします。多角形で囲まれた領域で rdrdθ を積分すると面積が得られるので、この手順で正しく面積が得られるはずであると結論付けます。
円と多角形の交点の面積
次の図に示すように、半径 R の円と多角形との交点の面積を求める方法について説明します。
緑の領域の面積を見つけることに興味があります。単一の多角形の場合と同様に、計算を分割して多角形の各辺の面積を見つけ、それらの面積を合計します。
最初の領域は次のようになります。
2番目の領域は次のようになります
そして3つ目のエリアは
この場合も、最初の 2 つの領域はプラスですが、3 番目の領域はマイナスになります。キャンセルがうまくいき、正味の面積が実際に私たちが関心のある面積になることを願っています。見てみましょう。
実際、面積の合計が、関心のある面積になります。
繰り返しますが、なぜこれが機能するのかについて、より厳密な説明を行うことができます。交差によって定義される領域を I とし、ポリゴンを 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 のイメージです。
再びストークスの定理を使用すると、r^2 dθ/2 を π(P) の境界で積分すると π(P) の面積が得られることがわかります。つまり、dxdy を π(P) で積分した場合と同じ答えが得られます。
再び変数の変更を使用すると、dxdy を π(P) で積分することは、Jdxdy を P で積分することと同じであることがわかります。ここで、J は π のヤコビアンです。
これで、Jdxdy の積分を 2 つの領域 (円内の部分と円外の部分) に分割できます。ここで、π は円内の点のみを残すため、そこで J=1 となるため、P のこの部分からの寄与は、円内にある P の部分の面積、つまり交点の面積になります。2 番目の領域は、円の外側の領域です。π がこの部分を円の境界まで折りたたむため、J=0 になります。
したがって、私たちが計算するのは実際に交差点の面積です。
領域を見つける方法を概念的に理解したので、単一のセグメントからの寄与を計算する方法についてより具体的に説明しましょう。私が「標準ジオメトリ」と呼ぶもののセグメントを見ることから始めましょう。以下に示します。
標準的なジオメトリでは、エッジは左から右に水平になります。これは 3 つの数値で表されます。xi はエッジの開始位置の x 座標、xf はエッジの終了位置の x 座標、y はエッジの y 座標です。
|y| の場合 < R、図のように、エッジは点 (-xint,y) と (xint,y) で円と交差します。ここで、xint = (R^2-y^2)^(1/2) です。次に、計算する必要がある面積は、図でラベル付けされた 3 つの部分に分割されます。領域 1 と 3 の面積を取得するには、arctan を使用してさまざまな点の角度を取得し、その面積を R^2 Δθ/2 と等しくします。たとえば、θi = atan2(y,xi) と θl = atan2(y,-xint) を設定します。この場合、領域 1 の面積は R^2 (θl-θi)/2 です。領域 3 の面積も同様に求めることができます。
領域 2 の面積はちょうど三角形の面積です。ただし、符号には注意が必要です。表示される面積を正にしたいので、面積は -(xint - (-xint))y/2 とします。
もう 1 つ注意すべき点は、一般に、xi が -xint より小さくなくても、xf が xint より大きくなくてもよいということです。
考慮すべきもう 1 つのケースは |y| です。> R. この場合は、図の領域 1 に類似するピースが 1 つしかないため、より単純です。
標準ジオメトリのエッジから面積を計算する方法がわかったので、あとはエッジを標準ジオメトリに変換する方法を説明するだけです。
しかし、これは単純な座標の変更にすぎません。最初の頂点 vi と最後の頂点 vf を持つものがあると、新しい x 単位ベクトルは vi から vf を指す単位ベクトルになります。この場合、xi は x に点線で囲まれた円の中心からの vi の変位であり、xf は xi に vi と vf の間の距離を加えたものです。一方、y は x と円の中心からの変位 vi のくさび積によって与えられます。
コード
これでアルゴリズムの説明は完了です。次はコードを記述します。ジャバを使用します。
まず、円を扱っているので、円クラスが必要です。
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にPathIteratorは、ポリゴンのエッジを反復処理するために使用できる があります。
では実際の作業です。これが完了したら、面積を計算するロジックから、エッジを反復してエッジを標準ジオメトリに配置するなどのロジックを分離します。この理由は、将来、領域のほかに、または領域に加えて何かを計算したい場合があり、エッジを反復処理する必要があるコードを再利用できるようにしたいからです。
Tしたがって、ポリゴン円の交差に関するクラスのプロパティを計算するジェネリック クラスがあります。
public abstract class CircleShapeIntersectionFinder<T> {
ジオメトリの計算に役立つ 3 つの静的メソッドがあります。
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を保持する と、正方形の半径のコピーを保持する の2 つのインスタンス フィールドがあります。これは奇妙に思えるかもしれませんが、私が使用しているクラスは実際には、円とポリゴンの交差のコレクション全体の領域を見つける機能を備えています。そのため、サークルの 1 つを「現在」と呼んでいます。
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は抽象的であり、1 つの実装は、平方根を実行する必要がないように円の半径を格納することです。とにかくに戻るprocessCircleShape。を呼び出した後initializeForNewCirclePrivate、多角形がnull(空の多角形と解釈しています) かどうかを確認し、そうであれば戻りnullます。この場合、計算された面積はゼロになります。ポリゴンがそうでない場合、ポリゴンnullの を取得しPathIteratorます。呼び出すメソッドの引数getPathIteratorは、パスに適用できるアフィン変換です。ただし、適用したくないので、渡すだけnullです。
double[]次に、頂点を追跡するを宣言します。は各頂点を 1 回しか与えないため、最初の頂点を覚えておく必要があります。そのためPathIterator、最後の頂点が与えられた後に戻って、この最後の頂点と最初の頂点でエッジを形成する必要があります。
次の行のcurrentSegmentメソッドは、次の頂点を引数に入れます。頂点がなくなったことを知らせるコードを返します。これが、while ループの制御式がそのままになっている理由です。
このメソッドの残りのコードのほとんどは、頂点の反復処理に関連する興味深いロジックではありません。重要なことは、while ループの反復ごとに 1 回呼び出し、メソッドの最後でもう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 単位ベクトルを取得するために必要だからです。この情報を取得したら、円の中心から最初の頂点までの変位を計算します。これと の内積は、私が xi と呼んでいたものをsegmentDisplacement与えてくれます。leftX次に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);
}
}
1 つ目は、エッジが円と交差するほど が小さいif場合を区別します。が大きく、交差する可能性がないy場合は、そのケースを処理するメソッドを呼び出します。yそれ以外は交差が可能な場合を扱います。
交差が可能である場合、交差の x 座標を計算しintersectionX、エッジを 3 つの部分に分割します。これらの部分は、上の標準ジオメトリ図の領域 1、2、および 3 に対応します。まず、領域 1 を処理します。
リージョン 1 を処理するには、leftXが実際に より小さいかどうかを確認-intersectionXします。そうでない場合、リージョン 1 はありません。リージョン 1 がある場合は、それがいつ終了するかを知る必要があります。rightXとの最小値で終了し-intersectionXます。これらの x 座標を見つけた後、この非交差領域を処理します。
領域 3 を処理するために同様のことを行います。
領域 2 については、それをチェックするロジックを実行する必要がleftXありrightX、実際に と の間の領域を囲み-intersectionXますintersectionX。領域を見つけた後は、領域の長さと のみが必要なyので、領域 2 を処理する抽象メソッドにこれら 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の間の角度の差を計算するために使用します。次に、 の不連続性を処理するコードを追加しますが、不連続性は 180 度または 0 度のいずれかで発生するため、これはおそらく不要です。次に、角度の違いを抽象メソッドに渡します。最後に、抽象メソッドとゲッターのみを使用します。leftXrightXatan2
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;
}
とにかく、それがアルゴリズムの私の説明です。正確で、チェックするケースがそれほど多くないので、いいと思います。
