0

我有一个表面上有点的 3D 球体。这些点以球坐标表示,如方位角、仰角和 r。

例如,我的数据集是一个矩阵,其中包含给定球体上的所有可用点:

  azimuth       elevation     r
[[  0.          90.           1.47      ]
 [  0.          85.2073787    1.47      ]
 [  0.          78.16966379   1.47      ]
 [  0.          70.30452954   1.47      ]
 [  0.          62.0367492    1.47      ]
 [  0.          53.56289304   1.47      ]
 [  0.          45.           1.47      ]
 [  0.          36.43710696   1.47      ]
 [  0.          27.9632508    1.47      ]
 [  0.          19.69547046   1.47      ]
 [  0.          11.83033621   1.47      ]
 [  0.           4.7926213    1.47      ]
 [  0.           0.           1.47      ]
 [  0.          -4.7926213    1.47      ]
 [  0.         -11.83033621   1.47      ]
 [  0.         -19.69547046   1.47      ]
 [  0.         -27.9632508    1.47      ]
 [  0.         -36.43710696   1.47      ]
 [  0.         -45.           1.47      ]
 [  0.         -53.56289304   1.47      ]
 [  0.         -62.0367492    1.47      ]
 [  0.         -70.30452954   1.47      ]
 [  0.         -78.16966379   1.47      ]
 [  0.         -85.2073787    1.47      ]
 [  0.         -90.           1.47      ]
 [  1.64008341  -1.6394119    1.47      ]
 [  1.64008341   1.6394119    1.47      ]
 [  2.37160039   8.01881788   1.47      ]
 [  2.37160039  -8.01881788   1.47      ]
 [  2.80356493 -15.58649429   1.47      ]
 [  2.80356493  15.58649429   1.47      ]
 [  3.16999007  23.70233802   1.47      ]
 [  3.16999007 -23.70233802   1.47      ]
 [  3.56208248 -32.09871039   1.47      ]
 [  3.56208248  32.09871039   1.47      ]
 [  4.04606896  40.63141594   1.47      ]
 [  4.04606896 -40.63141594   1.47      ]
 [  4.1063771   -4.09587122   1.47      ]
ecc...

注意:我故意省略了完整的数据矩阵,因为它包含大量数据。如果需要/要求使问题完全可重现,我将提供完整数据。

该矩阵表示如下图像:

球体点

给定一个任意点,我想计算数据集中“包含”输入点的 3 个最近点。

到目前为止,我的代码如下:

def compute_three_closest_positions(self, azimuth_angle, elevation):
    requested_position = np.array([azimuth_angle, elevation, 0])

    # computing the absolute difference between the requested angles and the available one in the dataset
    result = abs(self.sourcePositions - requested_position) #subtracting between dataset and requested point
    result = np.delete(result, 2, 1) # removing the r data 
    result = result.sum(axis=1) #summing the overall difference for each point

    # returning index of the closest points
    indexes = result.argsort()[:3] 

    closest_points = self.sourcePositions[indexes]
    
    return closest_points 

基本上,我从矩阵数据集中的所有点中减去所请求的方位角和仰角self.sourcePositions

代码工作正常,问题是有时我得到 3 个不包含请求点的最近点。

例子:

错误一:

  Requested point: azimut, elevation, distance
    [200   0   1.47]
    # As you might notice, the requested point is not inside the triangle created by the 3 closest points
    Three closest points: azimut, elevation, distance
    [[199.69547046   0.           1.47      ]
     [199.40203659   5.61508214   1.47      ]
     [199.40203659  -5.61508214   1.47      ]]

好一个:

Requested position:
[190   0   1.47]
# As you can notice, in this case the requested point is inside the triangle generated by the closest 3 points
Three closest points:
[[191.83033621   0.           1.47      ]
 [188.02560265   2.34839855   1.47      ]
 [188.02560265  -2.34839855   1.47      ]]

我该如何解决这个问题?我想获得“三角形”(我在球面上,所以它不是真正的三角形)包含我请求的点的 3 个最近点。

4

1 回答 1

1

从一开始azimut+elevation听起来不适合这个,我更喜欢latitude+longitudetherms,因为你使用的是其他东西。

现在如评论中所述,如果您的点形成规则网格,您可以制作一个方程式,将拓扑映射到点并返回,其中拓扑由数组中点的整数索引(或 2 个索引)描述。但是,如果您无法推断这一点或网格不规则,那么您可以这样做:

  1. 重新排序您的数组,使其按纬度和经度排序

    所以经度<0,2*PI>和纬度在<-PI/2,+PI/2>所以我会按两者对你的数组进行排序,而纬度具有更大的优先级(权重)。让我们调用这个新数组pnt[]

  2. 映射点p到最近的球体顶点

    简单地进行二分搜索pnt[],直到找到ix具有最大小于或等于纬度的点索引p

    然后从线性搜索ix(如果您重新排序pnt[]切片或记住每个纬度有多少点,则可以使用二进制搜索)直到找到小于或等于的最大纬度p

    现在pnt[ix]是球体上“最接近”的点p

  3. 列出最近的邻居pnt[ix]

    所以只是pnt[ix]从经度的左侧开始,所以pnt[ix+1]应该是下一个点(如果你越过数组大小你或极需要用蛮力检查点,但只是数组中的最后几个)

    现在我们只需要找到低于或高于这两个点的对应点(取决于你在哪一边p)。p因此,以相同的方式找到 2 个最接近的点,#2但纬度越来越小(前一个切片,后一个切片)。这将为您3*2 points提供 4 个潜在三角形的形式(始终使用首先找到的 2 个点)。

  4. 测试可能的三角形

    所以你有潜在的三角形p0,p1,p2,它与球面“平行”。因此,只需将您的点投影到其平面上:

    u = p1-p0
    v = p2-p0
    v = cross(u,v)
    v = cross(u,v)
    p' = p0 + dot(p-p0,u)*u + dot(p-p0,v)*v
    

    u,v基向量也是如此,并且p0是平面的原点......现在只需测试p'在三角形内,因此要么使用 2D 和重心坐标,要么利用叉积并检查 CW/CCW,如:

     if (dot(cross(p1-p0,p'-p0),p')>=0)
      if (dot(cross(p2-p1,p'-p1),p')>=0)
       if (dot(cross(p0-p2,p'-p2),p')>=0)
        point_is_inside;
     if (dot(cross(p1-p0,p'-p0),p')<=0)
      if (dot(cross(p2-p1,p'-p1),p')<=0)
       if (dot(cross(p0-p2,p'-p2),p')<=0)
        point_is_inside;
    

    因此,如果三角形的所有 3 个边都具有相同的 CW/CCW ness,p'那么您就找到了三角形。

于 2021-01-14T08:35:39.053 回答