球の表面に分布するポイントのボロノイ分割を見つけることを含む問題に取り組んでいます。私が知る限り、視覚的には点の Delaunay 三角形分割を見つけるように見えるので、私の力ずくのアプローチはうまくいきます。ただし、頂点を使用して各面のエッジの順序を定義すると、アルゴリズムが失敗するようです。
私が目指していることの例として、2 つの頂点が複数の形成点を共有しているかどうかを判断することによってエッジを決定するハックを使用して、エッジを正しく決定するバージョンの写真を次に示します。面の立体角を計算し、OpenGL のような 3D レンダリング API 用のジオメトリを生成できるようにするためにテッセレーションを使用したいので、このハックでは十分ではないことに注意してください。
赤い円は、球の表面に分布する点です。黄色の線はこれらの点の Delaunay 三角形分割を示し、緑色の線はボロノイ セル間の頂点を定義するために使用される点を示し、黒色の線は頂点によって形成されるエッジを示します。各セルは、点または線の近くにない各ピクセルを、セルの定義点を色に変換することによって決定される色に設定することによって色付けされます。これは、テッセレーション プロセスとは別に実行されます。面の色の値を比較するにはツールを使用する必要があるかもしれませんが、面が面に正しく囲まれていることを示すことができます。これは、私のコードがドロネー三角形分割とボロノイ分割の頂点を正しく決定していることを示しているようです。
ハックを取り除き、顔のポイントを反時計回りに並べ替えるために書いた関数を使用すると、説明できない結果が得られます。私のプログラムを実行するたびに異なるランダム ポイントのセットが生成されるため、これら 2 つの図は同じポイント分布を表すものではないことに注意してください。
問題を示す顔の周りに赤いボックスを描きました。これらのセルには面を通る黒い線があり、一部のエッジがまったく表示されない可能性があることに注意してください (右下のボックスを参照)。
この StackOverflow questionで説明されているアルゴリズムを使用して、ポイントの反時計回りの順序を決定しています。同じ関数を使用して、セルの周囲の頂点の順序を決定し、3 点の外心を決定します。コードにバグがある場合、コードが 3 点の場合に失敗することが予想され、その結果、Delaunay テッセレーションに問題が発生します (順序にエラーがあると、外心を反対側に配置することになるため)。しかし、何十回もの実行でクラッシュしたり、Delaunay テッセレーションの欠陥を明らかにしたりしたことはありません。何時間もコードと格闘しましたが、問題が見つかりません。誰かがこの問題が発生する理由を理解できますか?
以下は、コードの要約リストであり、すべての重要な点がリストされていることを願っています。これは、何かを機能させるために私が書いた複数のファイルからのコードの融合です。アルゴリズムが機能するまで、コードをクリーンアップしようとしない傾向があります。また、使用されていない場合は、インクルードまたは必要なインターフェイス メソッドの実装も入れませんでした。
public class SphericalVoronoiTessellation {
private Map<Point, List<Point>> faces = new HashMap<>();
private Set<Pair<Point, Point>> edges = new HashSet<>();
private Set<Pair<Point, Point>> neighbors = new HashSet<>();
private Map<Point, Set<Point>> vertices = new HashMap<>();
public SphericalVoronoiTessellation(List<Point> points) {
List<Point> copy = new ArrayList<>(points);
Collections.sort(copy);
for (Point p : copy) {
faces.put(p, new ArrayList<Point>());
}
final int n = points.size();
for (int i = 0; i < n - 2; i++) {
Point p = copy.get(i);
for (int j = i + 1; j < n - 1; j++) {
Point q = copy.get(j);
for (int k = j + 1; k < n; k++) {
Point r = copy.get(k);
Point c = getCircumcenter(p, q, r);
double d = p.getSphericalDistanceTo(c);
if (circleIsEmpty(c, d, i, j, k, copy)) {
faces.get(p).add(c);
faces.get(q).add(c);
faces.get(r).add(c);
neighbors.add(pair(p, q));
neighbors.add(pair(p, r));
neighbors.add(pair(q, r));
Set<Point> formedBy;
if (!vertices.containsKey(c)) {
formedBy = new HashSet<>();
vertices.put(c, formedBy);
} else {
formedBy = vertices.get(c);
}
formedBy.add(p);
formedBy.add(q);
formedBy.add(r);
}
}
}
}
// TODO: Determine why using getCounterClockwiseOrder does not correctly
// order the vertices. It seems to correctly order three vertices
// every time, but that might just be luck...
for (Map.Entry<Point, List<Point>> face : faces.entrySet()) {
List<Point> vertices = getCounterClockwiseOrder(face.getValue());
// Store the vertices in the counter-clockwise order so that they
// can be used to determine the face's surface.
faces.put(face.getKey(), vertices);
// Builds a set of edges for the whole diagram. I use this set for
// duplicate-free testing of the edges on the diagram.
for (int k = 0; k < vertices.size(); k++) {
Point a = vertices.get(k);
Point b = vertices.get(k + 1 == vertices.size() ? 0 : k + 1);
edges.add(pair(a, b));
}
}
}
private static Point getCircumcenter(Point a, Point b, Point c) {
List<Point> ccw = new ArrayList<Point>();
ccw.add(a);
ccw.add(b);
ccw.add(c);
ccw = getCounterClockwiseOrder(ccw);
return
getPlaneNormal(
ccw.get(2),
ccw.get(1),
ccw.get(0)
).times(a.getRadius());
}
// This function is the one that may be broken...
private static List<Point> getCounterClockwiseOrder(List<Point> points) {
List<Point> ordered = new ArrayList<Point>(points);
final Point c = getCentroid(points);
final Point n = c.getNormalized();
final Point s = points.get(0);
final Point toS = s.minusCartesian(c);
Collections.sort(
ordered,
new Comparator<Point>() {
@Override
public int compare(Point o1, Point o2) {
if (o1.equals(o2)) {
return 0;
} else {
return Double.compare(
getDistanceFromS(o1),
getDistanceFromS(o2)
);
}
}
private double getDistanceFromS(Point p) {
if (s.equals(p)) {
return 0;
}
double distance = s.getSphericalDistanceTo(p);
Point toP = p.minusCartesian(c);
Point cross = toS.cross(toP);
if (n.dot(cross) < 0) {
distance = RotationDisplacement.REVOLUTION - distance;
}
return distance;
}
}
);
return ordered;
}
private static Point getCentroid(List<Point> points) {
Point centroid = Point.ORIGIN;
for (Point p : points) {
centroid = centroid.plus(p);
}
return centroid.times(1. / points.size());
}
private static Point getPlaneNormal(Point a, Point b, Point c) {
Point d = a.minusCartesian(b);
Point e = c.minusCartesian(b);
return d.cross(e).getNormalized();
}
private static boolean circleIsEmpty(
Point center,
double distance,
int i,
int j,
int k,
List<Point> points
) {
int m = 0;
for (; m < points.size(); m++) {
if (m == i || m == j || m == k) {
continue;
}
if (center.getSphericalDistanceTo(points.get(m)) < distance) {
break;
}
}
return m == points.size();
}
private static Pair<Point, Point> pair(Point a, Point b) {
if (b.compareTo(a) < 0) {
Point swap = b;
b = a;
a = swap;
}
return new ImmutablePair<Point, Point>(a, b);
}
}
public class Point implements Comparable<Point> {
private double radius;
private RotationDisplacement spherical;
private VectorDisplacement cartesian;
public Point(VectorDisplacement coordinates) {
this.cartesian = coordinates;
this.calculateSpherical();
}
public Point(double radius, RotationDisplacement rotations) {
this.radius = Math.abs(radius);
if (radius < 0) {
rotations = rotations.getNormalizedRepresentation();
rotations = new RotationDisplacement(
Math.PI - rotations.getColatitude(),
rotations.getLongitude() + Math.PI
);
}
this.spherical = rotations.getNormalizedRepresentation();
this.calculateCartesian();
}
private void calculateSpherical() {
this.radius = Math.sqrt(
this.getX() * this.getX() +
this.getY() * this.getY() +
this.getZ() * this.getZ()
);
double c =
this.radius > 0 ?
Math.acos(this.getY() / this.radius) :
0;
double l =
c > 0 && c < Math.PI ?
Math.atan2(-this.getZ(), this.getX()) :
0;
this.spherical =
new RotationDisplacement(
c,
l
).getNormalizedRepresentation();
}
public double getX() {
return this.cartesian.getX();
}
public double getY() {
return this.cartesian.getY();
}
public double getZ() {
return this.cartesian.getZ();
}
private void calculateCartesian() {
this.cartesian = new VectorDisplacement(
this.radius * Math.cos(
this.getLongitude()) * Math.sin(this.getColatitude()
),
this.radius * Math.cos(this.getColatitude()),
this.radius * -Math.sin(
this.getLongitude()) * Math.sin(this.getColatitude()
)
);
}
public double getLongitude() {
return this.spherical.getLongitude();
}
public double getColatitude() {
return this.spherical.getColatitude();
}
public Point plus(Point that) {
return new Point(
(VectorDisplacement) this.cartesian.add(that.cartesian)
);
}
public Point times(double scalar) {
return new Point(this.radius * scalar, this.spherical);
}
public Point getNormalized() {
return new Point(1, this.spherical);
}
public Point minusCartesian(Point that) {
return new Point(
(VectorDisplacement) this.cartesian.subtract(that.cartesian)
);
}
public double getSphericalDistanceTo(Point that) {
if (this.radius == 0 || that.radius == 0) {
return 0;
}
return this.radius * Math.abs(
Math.acos(this.dot(that) / (this.radius * that.radius))
);
}
public double dot(Point that) {
return
this.getX() * that.getX() +
this.getY() * that.getY() +
this.getZ() * that.getZ();
}
@Override
public boolean equals(Object other) {
if (!(other instanceof Point)) {
return false;
}
Point that = (Point) other;
return
this.cartesian.equals(that.cartesian) ||
this.radius == that.radius &&
this.spherical.equals(that.spherical);
}
public Point cross(Point that) {
double ux = this.getX();
double uy = this.getY();
double uz = this.getZ();
double vx = that.getX();
double vy = that.getY();
double vz = that.getZ();
return new Point(
new VectorDisplacement(
uy * vz - uz * vy,
uz * vx - ux * vz,
ux * vy - uy * vx
)
);
}
}
public interface Displacement {
public Displacement add(Displacement that);
public Displacement subtract(Displacement that);
public Displacement scale(double coefficient);
}
public class VectorDisplacement implements Displacement {
private double x;
private double y;
private double z;
public VectorDisplacement(double x, double y, double z) {
this.x = x;
this.y = y;
this.z = z;
}
public double getX() {
return x;
}
public double getY() {
return y;
}
public double getZ() {
return z;
}
@Override
public Displacement add(Displacement that) {
if (!(that instanceof VectorDisplacement)) {
throw new IllegalArgumentException(
"VectorDisplacement.add needs a VectorDisplacement"
);
}
VectorDisplacement other = (VectorDisplacement) that;
return new VectorDisplacement(
this.x + other.x,
this.y + other.y,
this.z + other.z
);
}
@Override
public boolean equals(Object other) {
if (!(other instanceof VectorDisplacement)) {
return false;
}
VectorDisplacement that = (VectorDisplacement) other;
return this.x == that.x && this.y == that.y && this.z == that.z;
}
@Override
public Displacement subtract(Displacement that) {
if (!(that instanceof VectorDisplacement)) {
throw new IllegalArgumentException(
"VectorDisplacement.subtract needs a VectorDisplacement"
);
}
VectorDisplacement other = (VectorDisplacement) that;
return new VectorDisplacement(
this.x - other.x,
this.y - other.y,
this.z - other.z
);
}
}
public class RotationDisplacement implements Displacement {
public static double REVOLUTION = Math.PI * 2;
private double colatitude;
private double longitude;
public RotationDisplacement(double colatitude, double longitude) {
this.colatitude = colatitude;
this.longitude = longitude;
}
public double getColatitude() {
return this.colatitude;
}
public double getLongitude() {
return this.longitude;
}
public RotationDisplacement getNormalizedRepresentation() {
double c = clampAngle(colatitude);
double l = 0;
if (c != 0 && c != Math.PI) {
if (c > Math.PI) {
c = RotationDisplacement.REVOLUTION - c;
l = Math.PI;
}
l = clampAngle(longitude + l);
}
return new RotationDisplacement(c, l);
}
@Override
public boolean equals(Object other) {
if (!(other instanceof RotationDisplacement)) {
return false;
}
RotationDisplacement my = this.getNormalizedRepresentation();
RotationDisplacement his =
((RotationDisplacement) other).getNormalizedRepresentation();
return
my.colatitude == his.colatitude &&
my.longitude == his.longitude;
}
private double clampAngle(double radians) {
radians %= RotationDisplacement.REVOLUTION;
if (radians < 0) {
radians += RotationDisplacement.REVOLUTION;
}
return radians;
}
}
この特定の問題を修正するための洞察をいただければ幸いです。