sorting a list of 3D points in clockwise order

Question

I have a list of 3D points and a center and I want to sort them in (counter)clockwise order around a given normal vector. The points are not coplanar but they and the center are bound to the surface of a sphere and they outline a polygon. the normal vector is the vector from the center of the sphere to the center of the sorting center. I tried this comparison function, but it fails when two points are more than π/2 apart from each other.

How do I get an actual 3D (counter)clockwise ordering for an arbitrary set of points?

This is not a duplicate of Sorting 3D points on the surface of a sphere in clockwise order because this question is specifically about dealing with the lack of transitivity in angular comparisons.

This is not a duplicate of Sorting a List of 3d coplanar points to be clockwise or counterclockwise because that question is more about determining if a point is closer to being clockwise or counterclockwise from another, which while it’s a comparison relation, it does not give a well-defined total ordering.

Answer 1

As you've figured out, a single dot product can't work by itself because it's a scalar cosine, and every value of cosine corresponds to two points of the unit circle.

So one way to a solution is to find two perpendicular reference vectors in the plane given by the normal and take triple products with those. They'll be the sine and cosine of an angle you can use for sorting. So you can either use atan2(y,x) to get an exact angle, or - if speed is important - approximate atan2/(pi/4) with slope and inverse slope.

To get the two vectors you need, first take the longest cross product I x n, J x n and K x n where I, J, K are the unit axis vectors. Call this vector p. It must lie in the plane because it's perpendicular to n. (You're taking the longest to avoid floating point accuracy issues.)

Now compute q = n x p. This also lies in the plane because it's perpendicular to n, but it's also perpendicular to p ... exactly what we need.

To recap, p and q are perpendicular vectors in any plane for which n is a normal.

Now if c is the center, for each point r in the polygon, compute triple products t = n * ((r - c) x p) and u = n * ((r - c) x q). Then atan2(u, t) or its approximation is a sorting metric.

Demo

Just to show this does work, including the atan2 approximation:

public class Sorter3d {

  // Sorting key metric calculator.
  static class Order {
    final Vec n, pp, qp;
    final Pt c;

    Order(Vec n, Pt c) { 
      this.c = c;
      this.n = n;
      pp = n.cross(Vec.I).longer(n.cross(Vec.J)).longer(n.cross(Vec.K));
      qp = n.cross(pp);
    } 

    double getKey(Pt r) {
      Vec rmc = r.minus(c);
      return approxAtan2(n.dot(rmc.cross(pp)), n.dot(rmc.cross(qp)));
    }
  }

  // Affine 3d vectors.
  static class Vec {
    static final Vec I = Vec.of(1, 0, 0);
    static final Vec J = Vec.of(0, 1, 0);
    static final Vec K = Vec.of(0, 0, 1);
    final double x, y, z;
    private Vec(double x, double y, double z) { this.x = x; this.y = y; this.z = z; }
    static Vec of(double x, double y, double z) { return new Vec(x, y, z); }
    Vec cross(Vec o) { return Vec.of(y * o.z - z * o.y, z * o.x - x * o.z, x * o.y - y * o.x); }
    double dot(Vec o) { return x * o.x + y * o.y + z * o.z; }
    double dot(Pt o) { return x * o.x + y * o.y + z * o.z; }
    double len2() { return dot(this); }
    double len() { return Math.sqrt(len2()); }
    Vec scale(double s) { return Vec.of(x * s, y * s, z * s); }
    Vec unit() { return scale(1.0 / len()); }
    Vec longer(Vec o) { return len2() > o.len2() ? this : o; }
    public String toString() { return String.format("[%.3f,%.3f,%.3f]", x, y, z); }
  }

  // Affine 3d points.
  static class Pt {
    static final Pt O = Pt.of(0, 0, 0);
    final double x, y, z;
    private Pt(double x, double y, double z) { this.x = x; this.y = y; this.z = z; }
    static Pt of(double x, double y, double z) { return new Pt(x, y, z); }
    Pt plus(Vec o) { return Pt.of(x + o.x, y + o.y, z + o.z); }
    Vec minus(Pt o) { return Vec.of(x - o.x, y - o.y, z - o.z); }
    public String toString() { return String.format("(%.3f,%.3f,%.3f)", x, y, z); }
  }

  // Return approximation of atan2(y,x) / (PI/2);
  static double approxAtan2(double y, double x) {
    int o = 0;
    if (y < 0) { x = -x; y = -y; o |= 4; }
    if (x <= 0) { double t = x; x = y; y = -t; o |= 2; }
    if (x <= y) { double t = y - x; x += y; y = t; o |= 1; }
    return o + y / x;
  }

  public static void main(String [] args) {
    // Make some random points radially sorted about the Z axis.
    int nPts = 17;
    Pt [] pts = new Pt[nPts];
    for (int i = 0; i < nPts; ++i) {
      double r = 1.0 + 10 * Math.random();
      double theta = i * (2 * Math.PI / nPts);
      pts[i] = Pt.of(r * Math.cos(theta), r * Math.sin(theta), 40.0 * (1 - Math.random()));
    }
    // Pick arbitrary normal vector and center point.
    // Rotate z-axis to normal and translate origin to center.
    Vec normal = Vec.of(-42.0, 17.0, -91.0);
    Vec cx = Vec.J.cross(normal).unit();
    Vec cy = normal.cross(cx).unit();
    Vec cz = normal.unit();
    Vec rx = Vec.of(cx.x, cy.x, cz.x);
    Vec ry = Vec.of(cx.y, cy.y, cz.y);
    Vec rz = Vec.of(cx.z, cy.z, cz.z);
    Pt center = Pt.of(11, 12, 13);
    Vec ofs = center.minus(Pt.O);
    Pt [] xPts = new Pt[nPts];
    for (int i = 0; i < nPts; ++i) {
      xPts[i] = Pt.of(rx.dot(pts[i]), ry.dot(pts[i]), rz.dot(pts[i])).plus(ofs);
    }
    // Check the sort keys returned by the sorter.
    Order order = new Order(normal, center);
    for (int i = 0; i < nPts; ++i) {
      System.out.println(order.getKey(xPts[i]));
    }
  }
}

This prints a valid key order:

4.0
3.9924071330572093
3.982224060033384
3.9612544376696253
3.8080585081381275
0.03457371559793447
0.013026386180392412
0.006090856009723169
0.0018388671161891966
7.99632901621898
7.987892035846782
7.974282237149798
7.93316335979413
4.106158894193932
4.019755500146331
4.008967674404233
4.003810901304664

Answer 2

Okay I’ve managed to find my own solution which uses only dot and cross products and no inverse trig or square roots or anything. You pick the first vertex v in your list and use that as your reference. Then you cross that vector r = v - center with the normal vector to get a half-space partition vector p. If the two inputs are on the same side of p then you can use the triple product without any problems because the cylindrical angle between them will be less than π. There’s some edge cases to look out for though so I figured I’d just share some pseudocode.

let c be the center around which the counterclockwise sort is to be performed
let n be the normal vector 

r := vertices[0] - c // use an arbitrary vector as the twelve o’clock reference
p := cross(r, c)     // get the half-plane partition vector

// returns true if v1 is clockwise from v2 around c
function less(v1, v2):
    u1 := v1 - c
    u2 := v2 - c
    h1 := dot(u1, p)
    h2 := dot(u2, p)

    if      h2 ≤ 0 and h1 > 0:
        return false

    else if h1 ≤ 0 and h2 > 0:
        return true

    else if h1 = 0 and h2 = 0:
        return dot(u1, r) > 0 and dot(u2, r) < 0

    else:
        return dot(cross(u1, u2), c) > 0

    //           h2 > 0     h2 = 0     h2 < 0
    //          ————————   ————————   ————————
    //  h1 > 0 |   *        v1 > v2    v1 > v2
    //  h1 = 0 | v1 < v2       †          *
    //  h1 < 0 | v1 < v2       *          *

    //  *   means we can use the triple product because the (cylindrical)
    //      angle between u1 and u2 is less than π
    //  †   means u1 and u2 are either 0 or π around from the zero reference
    //      in which case u1 < u2 only if dot(u1, r) > 0 and dot(u2, r) < 0

Answer 3

You project the points on a plane perpendicular to the normal (form an orthogonal frame from the normal). Then in this plane, use polar coordinates and sort by angles. Note anyway, that the zero of the angles is arbitrary.