使用 Java 7 进行球形 Voronoi 镶嵌:需要修复围绕面缠绕顶点的问题 [英] Spherical Voronoi Tessellation with Java 7: need fix for winding vertices around faces
问题描述
我正在解决一个问题，该问题涉及为分布在球体表面上的点找到 Voronoi 细分.据我所知，我的蛮力方法有效，因为在视觉上它似乎找到了点的 Delaunay 三角剖分.但是，在使用顶点定义每个面的边缘顺序时，我的算法似乎失败了.
I am working on a problem that involves finding the Voronoi tessellation for points distributed over the surface of a sphere. As far as I can tell, my brute-force approach works, since visually it seems to find the Delaunay triangulation of the points. However, my algorithm seems to fail when defining each face's edge order using the vertices.
作为我要使用的示例，这里是一个版本的图片，该版本使用 hack 正确确定边缘，该 hack 通过确定两个顶点是否共享多个形成点来确定边缘.请注意，我想使用曲面细分来计算面的立体角并为 OpenGL 等 3D 渲染 API 生成几何图形，因此这种 hack 还不够好.
As an example of what I'm going for, here is a picture of a version that correctly determines the edges using a hack that determines edges by determining whether two vertices share more than one forming point. Note that I want to use the tessellation to be able to calculate the solid angle of a face and to generate geometry for a 3D-rendering API like OpenGL, so this hack is not good enough.
红色圆圈是分布在球体表面上的点.黄线显示这些点的 Delaunay 三角剖分，绿线显示哪些点用于定义 Voronoi 单元之间的顶点，黑线显示顶点形成的边.通过将不靠近点或线的每个像素设置为通过将单元格的定义点转换为颜色来确定的颜色来为每个单元格着色；这是与曲面细分过程分开执行的.可能需要使用工具来比较面部颜色值，但可以显示面部被面部正确包围.这似乎表明我的代码正确确定了 Delaunay 三角剖分和 Voronoi 曲面细分的顶点.
The red circles are points distributed over the surface of the sphere. The yellow lines show the Delaunay triangulation for those points, the green lines show which points are used to define the vertices between the Voronoi cells, and the black lines show the edges formed by the vertices. Each cell is colored by setting each pixel not near a point or a line to a color determined by transforming the cell's defining point to a color; this is performed separately from the tessellation process. It may take using a tool to compare the face color values, but it can be shown that the faces are correctly enclosed by the faces. This seems to indicate that my code correctly determines the Delaunay triangulation and the vertices for the Voronoi tessellation.
当我删除 hack 并使用我编写的逆时针排序面部点的函数时，我得到了无法解释的结果.请注意，我的程序每次运行都会生成一组不同的随机点，因此这两个图并不打算表示相同的点分布.
When I remove the hack and use the function I wrote for ordering a face's points counter-clockwise, I get results I can't explain. Please note that every run of my program generates a different set of random points, so these two diagrams are not intended to represent the same point distribution.
我在说明问题的面孔周围画了红框.请注意，这些单元格的表面有黑色线条，可能会导致某些边缘根本无法显示(请参见右下框).
I've drawn red boxes around faces that demonstrate the problem. Note that these cells have black lines running through their faces, and can result in some edges not being represented at all (see the lower-right box).
我正在使用 这个 StackOverflow 问题 确定点的逆时针顺序.我使用相同的函数来确定单元格周围的顶点顺序和确定三点的外心.如果代码中存在错误，人们会期望代码在三点情况下失败，从而引入 Delaunay 曲面细分的问题(因为顺序中的错误会导致将外心放置在三点的另一侧)球体)，但数十次运行从未崩溃，也没有发现 Delaunay 曲面细分的任何缺陷.我已经与我的代码搏斗了几个小时，但我找不到问题所在.有人能明白为什么会出现这个问题吗?
I'm using an algorithm described at this StackOverflow question to determine the counter-clockwise order of points. I use the same function for determining the vertex order around cells and for determining the circumcenter of three points. If there is a bug in the code, one would expect the code to fail for the three-point case and thus introduce a problem with the Delaunay tessellation (since an error in the order would result in placing the circumcenter on the opposite side of the sphere), yet dozens of runs have never crashed nor revealed any flaws with the Delaunay tessellation. I've wrestled with my code for hours, and I cannot find the problem. Can somebody see why this problem occurs?
以下是代码的汇总列表，我希望列出所有要点.它是我编写的多个文件的混合代码，试图让某些东西正常工作；在我的算法起作用之前，我倾向于不尝试清理代码.如果不使用，我也没有放入包含或必需的接口方法实现.
The following is a summarized listing of the code with what I hope are all the salient points listed. It is an amalgam of code from multiple files I wrote trying to get something working; I tend not to try to clean up the code until my algorithm works. I also did not put in the includes or required interface method implementations if they're not used.
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;
}
}
对解决此特定问题的任何见解将不胜感激.
Any insight into fixing this specific problem would be appreciated.
推荐答案
自然需要求助自己才能看到解决方案
Naturally, it takes asking for help to see the solution oneself <sighs>.
问题是我使用从球体中心到表面(顶点坐标)的向量来确定顶点之间的角度，而不是从面的质心到顶点的向量.后一种方法将给出 [0, 2 * PI) 范围内的结果，因为点都围绕质心旋转，而前一种方法只是检索顶点之间的大圆距离.
The problem is that I am using the vectors from the center of the sphere to the surface (the vertices' coordinates) for determining the angle between the vertices as opposed to the vector from a face's centroid to the vertices. The latter approach will give a result in the range of [0, 2 * PI), since the points rotate all around the centroid, whereas the prior approach just retrieves the great circle distance between the vertices.
我已经修正了 getCounterClockwiseOrder 方法如下，现在它可以工作了.如果其他人正在寻找如何使用 Java 确定球形 Voronoi 细分，我会留下这个问题.
I've fixed the getCounterClockwiseOrder method as follows, and now it works. I'll leave this question up in case somebody else is looking for how to determine a spherical Voronoi tessellation with Java.
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;
}
这篇关于使用 Java 7 进行球形 Voronoi 镶嵌:需要修复围绕面缠绕顶点的问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
