球的Voronoi镶嵌与Java 7:需要修复的缠绕面顶点顶点、Voronoi、Java

2023-09-08 01:14:13 作者:我請客你買單

我的工作,涉及到寻找沃罗努瓦​​剖分为分布在球体的表面点的一个问题。据我所知,我的蛮力方法有效,因为在视觉上似乎找到点的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.

由于我什么都为一个例子,这里是一个版本,可以正确地判断使用黑客工具,通过确定是否两个顶点共享多个形成点确定的边缘边缘的图片。请注意,我想用镶嵌才能够计算出脸部的立体角,并产生几何形状像OpenGL的一个3D渲染API,所以这个黑客是不够的。

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单元之间的顶点,和黑线显示由顶点所形成的边缘。每个小区通过将每个像素不靠近点或线来转化细胞的限定指向一个颜色确定的颜色着色;这是从分割过程分开进行。这可能需要使用工具来比较脸色值,但是可以示出,该面​​由面正确地包围。这似乎表明,我的code正确决定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.

当我删除了黑客和使用我写的订购一张脸的逆时针点的功能,我得到的结果我无法解释。请注意,我的节目的每一次运行产生一组不同的随机点,因此这两个图不旨在重新present相同点分布

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.

我身边的面孔,展示问题绘制红色方框。需要注意的是,这些细胞具有通过其面运行的黑线,并且可能会导致一些边不被psented在所有重$ P $(见右下框)。

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).

我使用的this StackOverflow的问题确定点的逆时针顺序。我用的是同样的功能,以确定细胞周围的顶点顺序并确定三点外心。如果在code中的错误,人们可能会期望code失败了三点情况,因此引进与德劳内镶嵌有问题(因为在订货会导致将外心错误在球体的相反侧),但数十奔跑从未应声也不透露与德洛奈镶嵌任何瑕疵。我斟酌着我的code几个小时,我不能发现问题。 有人可以明白为什么这个问题会出现?

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?

下面是在code什么希望都列出的要点汇总清单。它是从多个文件我写信试图让一些工作code的混合体;我倾向于不要试图清理code,直到我的算法工作。我还没有把在包括或所需的接口方法实现的,如果他们不使用。

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;
    }
}

任何洞察到解决这个特定的问题将是AP preciated。

Any insight into fixing this specific problem would be appreciated.

推荐答案

当然,它需要寻求帮助,看解决自己和LT;叹息&GT;

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;
}