Spherical Voronoi Tessellation with Java 7: need f

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.

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.

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.

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

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

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.

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