在Python中将点云分成象限

Question

假设我有一个由30个三维点组成的数组，用Python中的numpy数组表示：

import numpy as np
b = np.round(np.random.random((30,3))*20)
b = b - b.mean(axis=0)

我希望将这些点以点云的“重心”为参考分成八个象限。我可以这样做：

for a in (b, -b):
    q_list = a[a[:,0]>0.][a[a[:,0]>0.][:,1]>0][a[a[:,0]>0.][a[a[:,0]>0.][:,1]>0][:,2]>0]
    quad_list.append(q_list * sign)
    q_list = a[a[:,0]>0.][a[a[:,0]>0.][:,1]>0][a[a[:,0]>0.][a[a[:,0]>0.][:,1]>0][:,2]<=0]
    quad_list.append(q_list * sign)             
    q_list = a[a[:,0]>0.][a[a[:,0]>0.][:,1]<=0][a[a[:,0]>0.][a[a[:,0]>0.][:,1]<=0][:,2]>0]
    quad_list.append(q_list * sign)              
    q_list = a[a[:,0]>0.][a[a[:,0]>0.][:,1]<=0][a[a[:,0]>0.][a[a[:,0]>0.][:,1]<=0][:,2]<=0]
    quad_list.append(q_list * sign)
    sign *= -1

当然，这个函数可以正常工作。它返回一个包含八个位置坐标数组的列表，每个数组只包含落在特定象限中的点。但我觉得应该有更清晰、更简洁的方法来处理这个问题。您有什么建议吗？

- carsonc

---

类似这样的代码可能很不错：np.packbits(b > 0.0, axis=1) >> 5。它会将象限标记为0-7，但可能不是所需的顺序。 - user2379410

如果它没有按照所需的顺序标记它们，那么将期望的顺序映射到实际编号上的查找表添加起来将相对容易。 - holdenweb

Answer 1

这里有一个递归解决方案。它适用于任意维度的情况。

import numpy as np

def split_into_quadrants(points, idx=0):
    if idx < points.shape[-1]:
        positive = points[points[:, idx] >= 0]
        negative = points[points[:, idx] < 0]
        return (split_into_quadrants(positive, idx+1) +
                split_into_quadrants(negative, idx+1))
    else:
        return [points]

b = np.round(np.random.random((30,3))*20)
b = b - b.mean(axis=0)
print(split_into_quadrants(b))

- YXD

---

谢谢。我刚刚对照了我的已更正的代码进行了测试，结果是一致的。这正是我希望学到的东西。我甚至从未想过递归方法。 - carsonc

Answer 2

递归对于大型数组来说可能计算量过大。

更好的选择可能是一个函数，它可以将点本身分成四个部分返回，或者仅返回四分之一的索引而不使用递归（仍然应该适用于任意维数）。

import numpy as np

def split_into_quadrants(points, get_quadrants_index=False):
   ndims = points.shape[1]
   existing_quadrants = list(itertools.product([-1, 1], repeat=ndims))
   quadrants_dict = dict(zip(existing_quadrants,range(len(existing_quadrants))))
   sign_points = np.sign(points)
   quadrants_idx = np.apply_along_axis(lambda x: quadrants_dict[(tuple(x.astype(int)))], 1, sign_points)
   if get_quadrants_index:
       return quadrants_idx
   points_split_into_quadrants = []
   for i in set(quadrants_idx):
       points_split_into_quadrants.append(points[quadrants_idx==i])
   return points_split_into_quadrants

b = np.round(np.random.random((30,3))*20)
b = b - b.mean(axis=0)
print(split_into_quadrants(b,get_quadrants_index=False))

- Baragorn