计算在一组x，y，z坐标之间的值的数量

Question

我正试图写一个方法，让我能够计算物体的数量，在三维空间中，落在另一个物体的三维坐标。你可以说这个有值的物体也有半径，所以我试着计算一个球体内物体的数量。你知道吗

我不会发布我当前的脚本，但我将用一个例子来尝试：我有一个具有三维坐标gal_pos和半径gal_rad的星系。你知道吗

import numpy as np
gal_pos = np.array(
  [[ 528.1373291 ,  432.18615723,  443.8348999 ],
   [ 540.12231445,  450.08154297,  442.07891846],
   [ 590.73675537,  234.6769104 ,  296.02798462],
   [ 529.98809814,  161.75544739,  567.58203125],
   [ 552.45446777,  312.1973877 ,  375.42492676],
   [ 700.94335938,   65.46828461,  172.71842957],
   [ 463.43258667,   73.57706451,  285.4147644 ],
   [ 547.74414062,  330.9855957 ,  401.49771118],
   [ 591.89801025,  196.19670105,  274.60073853],
   [ 581.28320312,  376.70013428,  359.81851196],
   [ 520.09820557,  302.17849731,  371.68771362],
   [ 812.84539795,   97.41672516,  150.87428284],
   [ 541.6552124 ,   17.40070724,  373.07562256],
   [ 523.34509277,  302.18151855,  503.6333313 ]])

gal_rad = np.array(
  [ 1.14752779,  1.02471195,  0.79648002,  0.6085083 ,  0.78725676,
    1.07809084,  0.57744866,  0.93733404,  0.76053329,  0.68979678,
    0.61188519,  1.07989271,  0.83872035,  0.59899661])

然后我还得到了具有3D位置的星星star_pos。你知道吗

star_pos = np.array(
  [[ 517.0300293 ,  264.54165649,  547.87835693],
   [ 530.37280273,  358.40835571,  455.68734741],
   [ 530.42211914,  358.20803833,  455.80908203],
   [ 530.86737061,  324.91717529,  407.96405029],
   [ 547.05175781,  333.9262085 ,  403.82403564],
   [ 530.61053467,  325.91259766,  407.04153442],
   [ 533.9979248 ,  331.18804932,  451.3795166 ],
   [ 531.20678711,  326.75308228,  406.44711304],
   [ 550.81237793,  340.88101196,  408.75830078],
   [ 519.52880859,  299.91259766,  516.25140381],
   [ 525.82739258,  301.46209717,  501.66738892],
   [ 524.87988281,  268.88357544,  510.0123291 ],
   [ 524.43371582,  299.99725342,  512.36077881],
   [ 524.40429688,  299.8979187 ,  512.57452393],
   [ 524.40765381,  299.89120483,  512.5032959 ],
   [ 545.57440186,  331.59066772,  401.20291138],
   [ 532.29016113,  306.27557373,  491.26434326],
   [ 530.77410889,  326.18057251,  407.06216431],
   [ 524.14819336,  306.60586548,  509.55993652]])

以上只是我拥有的价值观的一个样本。你知道吗

xmax_rad = gal_pos[:,0]+gal_rad
xmin_rad = gal_pos[:,0]-gal_rad

ymax_rad = gal_pos[:,1]+gal_rad
ymin_rad = gal_pos[:,1]-gal_rad

zmax_rad = gal_pos[:,2]+gal_rad
zmin_rad = gal_pos[:,2]-gal_rad

tot_pop = [] # Total population found each galaxy

Nind = [(x,y,z) for x,y,z in enumerate(star_pos) 
        if any(xmin_rad <=x<= xmax_rad) and 
        any(ymin_rad<=y<=ymax_rad) 
        and any(zmin_rad<=x<=zmax_rad)]
tot_pop.append(Nind)

print tot_pop

我正在尝试的这个方法对我来说是最有意义的，它分解了所有的内容，但是它用于大小为~300的数组，但是返回ValueError: need more than 2 values to unpack用于Nind。可能是因为我的迭代无法解包3个对象？你知道吗

我也尝试过其他方法，在这些方法中，我获取每个位置的大小，但是它返回不正确的结果，以及通过直方图计算值，但同样返回不正确的结果（我通过在2d直方图中投影所有内容进行检查）。我为每个星系建立索引的方法会为每个星系返回空数组：

tot_pop = []
for k in np.arange(len(gal_pos)):
    Nind = [(x,y) for x,y in enumerate(star_pos) 
        if xmin_rad[k] <=x<= xmax_rad[k]) and 
        ymin_rad[k]<=y<=ymax_rad[k]]

    tot_pop.append(Nind)

Answer 1

这里有一种几乎矢量化的方法，利用有效的^{}和slicing来帮助解决问题-

# Define low and high limits
l = gal_pos - gal_rad[:,None]
h = gal_pos + gal_rad[:,None]

# Get mask of valid ones for each row of star_pos
mask = np.ones(star_pos.shape[0], dtype=bool)
for i in range(star_pos.shape[1]):
    mask &= ((l[:,i,None] <= star_pos[:,i]) & (h[:,i,None] >= star_pos[:,i])).any(0)

# Finally use the mask to select valid rows off star_pos
out = star_pos[mask]

称之为几乎矢量化，因为我们仍在迭代star_pos中的列数。但是，既然我们处理的是X,Y,Z数据，那就是3。因此，可以安全地称之为几乎矢量化。你知道吗

对于给定的样本，我得到的是-

In [302]: out
Out[302]: array([], shape=(0, 3), dtype=float64)

因此，没有来自star_pos的点满足极限。你知道吗

Answer 2

您可以使用zip遍历星系+半径，然后使用广播和布尔索引查找匹配项：

result = []
for galaxy, galaxy_radius in zip(gal_pos, gal_rad):
    # With broadcasting you can simply subtract the positions from the galaxy center
    # and using abs avoids checking lower and upper bound.
    rel_star_pos = abs(star_pos - galaxy)
    # Check which distances are below the radius and keep these which are
    # within the radius for x, y and z
    matches = (rel_star_pos <= galaxy_radius).all(axis=1)
    # use boolean indexing to append the stars which satisfy the above condition
    result.append(star_pos[matches])
print(result)

如果要附加索引（而不是实际的星坐标），可以将append行更改为：

result.append(np.where(matches)[0])

或者如果你只想知道匹配的数量：

result.append(np.sum(matches))

但是，我找不到与给定数据的任何匹配项：

[array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64),
 array([], shape=(0, 3), dtype=float64)]