我正在研究一个问题，该问题涉及找到分布在球体表面上的点的 Voronoi 镶嵌。据我所知，我的蛮力方法有效，因为在视觉上它似乎找到了点的 Delaunay 三角剖分。但是，在使用顶点定义每个面的边缘顺序时，我的算法似乎失败了。
作为我想要做的一个例子，这里是一个版本的图片，该版本使用一种通过确定两个顶点是否共享多个形成点来确定边缘的方法来正确确定边缘。请注意，我想使用曲面细分来计算面的立体角并为像 OpenGL 这样的 3D 渲染 API 生成几何图形，所以这个 hack 还不够好。
红色圆圈是分布在球体表面上的点。黄线显示这些点的 Delaunay 三角剖分，绿线显示用于定义 Voronoi 单元之间的顶点的点，黑线显示由顶点形成的边缘。通过将不靠近点或线的每个像素设置为通过将单元格的定义点转换为颜色来确定的颜色，为每个单元格着色；这是与镶嵌过程分开执行的。可能需要使用工具来比较面部颜色值，但可以显示面部正确地被面部包围。这似乎表明我的代码正确地确定了 Delaunay 三角剖分和 Voronoi 细分的顶点。
当我删除黑客并使用我编写的逆时针排序面部点的函数时，我得到了我无法解释的结果。请注意，我的程序的每次运行都会生成一组不同的随机点，因此这两个图表并非旨在表示相同的点分布。
我已经在展示问题的面孔周围绘制了红色框。请注意，这些单元格的表面有黑线，并可能导致根本不显示某些边缘（请参见右下角的框）。
我正在使用这个 StackOverflow 问题中描述的算法来确定点的逆时针顺序。我使用相同的函数来确定单元格周围的顶点顺序和确定三个点的外心。如果代码中存在错误，人们会认为代码在三点情况下会失败，从而引入 Delaunay 镶嵌问题（因为顺序错误会导致将外心放置在sphere)，但数十次运行从未崩溃，也没有发现 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;
}
}
任何有关解决此特定问题的见解将不胜感激。