按顺时针顺序对3D点列表进行排序

3
我有一个由3D点和中心组成的列表,我想按照给定法向量的(逆)顺时针顺序对它们进行排序。这些点不共面,但它们和中心都绑定在球体表面上,并且它们勾勒出一个多边形。法向量是从球心到排序中心的向量。我尝试了这个比较函数,但当两个点之间的距离超过π/2时会失败。
如何为任意一组点获取实际的3D(逆)顺时针排序? 这不是在球面上按顺时针顺序排序3D点的重复问题,因为这个问题具体涉及处理角度比较的非传递性。 这不是将3D共面点列表按顺时针或逆时针排序的重复问题,因为那个问题更多地涉及确定一个点是否更靠近另一个点的顺时针或逆时针,尽管它是一个比较关系,但它没有给出一个明确定义的总排序。
3个回答

6
正如你已经发现的那样,单个点积不能独立工作,因为它是标量余弦,每个余弦值对应于单位圆的两个点。
因此,解决方案之一是在由法线给出的平面上找到两个垂直的参考向量,并使用它们进行三重乘积。它们将是您可以用于排序的角度的正弦和余弦。因此,您可以使用atan2(y,x)来获取精确角度,或者-如果速度很重要-使用斜率和反斜率近似atan2 /(pi / 4)
要获取所需的两个向量,请先取最长的叉积I x nJ x nK x n,其中IJK是单位轴向量。将此向量称为p。它必须位于平面上,因为它垂直于n。(您正在避免浮点精度问题,因此选择最长的。)
现在计算q = n x p。这也位于平面上,因为它垂直于n,但它也垂直于p...正是我们需要的。
简而言之,pq是任何法线为n的平面中的垂直向量。
现在,如果c是中心,则对于多边形中的每个点r,计算三重积t = n * ((r-c) x p)u = n * ((r-c) x q)。然后,atan2(u,t)或其近似值是排序度量。
演示
只是为了展示这确实起作用,包括atan2的近似值:
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]));
    }
  }
}

这将打印一个有效的键顺序:

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

最后一次调用 atan2(u, t) 可以避免或替换为更简单的函数吗?我不需要角度,只需要总排序。 - taylor swift
就像我说的那样,在atan2的位置使用斜率和反斜率很容易。我会添加更多解释。 - Gene

4

好的,我已经找到了自己的解决方案,它只使用点积和叉积,没有反三角函数或平方根等。您选择列表中的第一个顶点v作为参考点。然后将该向量r = v-center与法向量相交,以获得半空间分割向量p。如果两个输入在p的同一侧,则可以使用三重积问题不大,因为它们之间的圆柱角度将小于π。但是需要注意一些边缘情况,所以我想分享一些伪代码。

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

啊,好的。我在想你需要一个以某个任意零点为极角的法线度量。除非比较函数的结果表示全序,否则许多排序方法都无法可靠地工作。像这样两两进行比较的排序,你必须找到一种打破循环的方法,否则要仔细选择排序算法。 - Gene
我应该补充一下,很久以前我实际上实现了一个基于循环顺序的比较函数,而且并没有多想,但是在我的系统上使用C的qsort()并没有可靠地返回循环顺序。 - Gene
没事。我看到你正在从参考向量中获取总序。请注意,如果参考向量的选择不好,这种方法可能会失败。(还要注意我的方法不需要平方根。我已经去掉了它们。) - Gene
Gene:哪个参考向量会有问题?有一个零分量的向量吗? - stephanmg
@taylor swift:我不知道你实际上在哪里使用普通的“n”。 - stephanmg
我认为 p := cross(r, c) 应该改为 p := cross(r, n) - Dave Durbin

0

将点投影到垂直于法线的平面上(从法线构建正交框架)。然后在该平面上,使用极坐标并按角度排序。请注意,角度的零点是任意的。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接