我正在研究一个问题，涉及找到分布在球体表面上的点的 Voronoi 曲面分割。据我所知，我的蛮力方法是有效的，因为从视觉上看，它似乎找到了点的 Delaunay 三角剖分。然而，当使用顶点定义每个面的边缘顺序时，我的算法似乎失败了。
作为我想要的示例，这里是一个版本的图片，该版本使用一种黑客方法正确确定边缘，该黑客通过确定两个顶点是否共享多个形成点来确定边缘。请注意，我希望使用曲面分割来计算面的立体角并为 OpenGL 等 3D 渲染 API 生成几何图形，因此这个 hack 还不够好。
红色圆圈是分布在球体表面的点。黄线显示这些点的 Delaunay 三角剖分，绿线显示哪些点用于定义 Voronoi 单元之间的顶点，黑线显示由顶点形成的边。通过将不靠近点或线的每个像素设置为通过将单元格的定义点转换为颜色而确定的颜色，对每个单元格进行着色；这是与曲面分割过程分开执行的。可能需要使用工具来比较面部颜色值，但可以显示面部被面部正确包围。这似乎表明我的代码正确确定了 Delaunay 三角剖分和 Voronoi 曲面分割的顶点。
当我删除黑客并使用我编写的用于逆时针排序脸部点的函数时，我得到了无法解释的结果。请注意，我的程序的每次运行都会生成一组不同的随机点，因此这两个图并不代表相同的点分布。
我在脸部周围画了红色框来说明问题。请注意，这些单元格的表面有黑线，可能会导致某些边缘根本无法显示(请参见右下框)。
我正在使用 this StackOverflow question 中描述的算法确定点的逆时针顺序。我使用相同的函数来确定单元格周围的顶点顺序并确定三个点的外心。如果代码中存在错误，人们会期望代码在三点情况下失败，从而引入 Delaunay 曲面分割的问题(因为顺序错误会导致将外心放置在曲面分割的另一侧)球体)，但数十次运行从未崩溃，也未暴露出 Delaunay 曲面分割的任何缺陷。我已经研究了我的代码几个小时，但找不到问题。 有人能明白为什么会出现这个问题吗？
以下是代码的摘要列表，我希望列出所有要点。它是我为了让某些东西正常工作而编写的多个文件中的代码的混合体；在我的算法起作用之前，我倾向于不会尝试清理代码。如果不使用它们，我也没有放入包含或所需的接口(interface)方法实现。
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;
}
}
任何有关解决此特定问题的见解将不胜感激。
最佳答案
当然，需要寻求帮助才能自己找到解决方案<叹息>。
问题是我使用从球体中心到表面(顶点坐标)的向量来确定顶点之间的角度，而不是使用从面的质心到顶点的向量。后一种方法将给出 [0, 2 * PI) 范围内的结果，因为点围绕质心旋转，而前一种方法仅检索顶点之间的大圆距离。
我已经修复了 getCounterClockwiseOrder 方法，如下所示，现在它可以工作了。如果其他人正在寻找如何使用 Java 确定球形 Voronoi 镶嵌，我将保留这个问题。
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).getNormalized();
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;
}
Point toP = p.minusCartesian(c).getNormalized();
double distance = toS.getSphericalDistanceTo(toP);
Point cross = toS.cross(toP).getNormalized();
if (n.dot(cross) < 0) {
distance = RotationDisplacement.REVOLUTION - distance;
}
return distance;
}
}
);
return ordered;
}
关于3d - 使用 Java 7 进行球形 Voronoi 曲面分割 : need fix for winding vertices around faces，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/19011166/