Python统计3D网格单元格内点数：np.histogramdd索引问题及替代方案咨询

我完全懂你的困扰——np.histogramdd确实能算出每个单元格的粒子数，但要把返回数组的索引和实际的三维单元格区间对应起来，实在太不直观了。下面几个方法能帮你清晰地看到每个单元格对应的空间范围和里面的粒子数量：

方法1：手动计算单元格标签，用Counter统计 #

这个方法直接给每个粒子打上所属单元格的区间标签，统计结果一目了然，完全不需要猜索引对应的位置。

import numpy as np
from collections import Counter

# 你的原始数据
testdata = np.array([[0.5,0.5,0.5],[0.6,0.6,0.6],[0.7,0.7,0.7],[1.5,0.5,0.5],[1.5,0.6,0.6],[0.5,1.5,0.5],[0.5,1.5,1.5]])

# 定义各维度的网格边界（可以根据需求调整）
x_edges = [0, 1, 2]
y_edges = [0, 1, 2]
z_edges = [0, 1, 2]

# 计算每个粒子在x/y/z维度上的单元格索引（减1让索引从0开始）
x_idx = np.digitize(testdata[:, 0], x_edges) - 1
y_idx = np.digitize(testdata[:, 1], y_edges) - 1
z_idx = np.digitize(testdata[:, 2], z_edges) - 1

# 生成每个粒子所属单元格的区间描述（直观易懂）
cell_labels = [
    f"x∈[{x_edges[i]},{x_edges[i+1]}), y∈[{y_edges[j]},{y_edges[j+1]}), z∈[{z_edges[k]},{z_edges[k+1]})"
    for i, j, k in zip(x_idx, y_idx, z_idx)
]

# 统计每个单元格的粒子数
cell_counts = Counter(cell_labels)

# 打印结果
for cell, count in cell_counts.items():
    print(f"{cell}: {count} 个粒子")

运行后会直接输出类似这样的结果：

x∈[0,1), y∈[0,1), z∈[0,1): 3 个粒子
x∈[1,2), y∈[0,1), z∈[0,1): 2 个粒子
x∈[0,1), y∈[1,2), z∈[0,1): 1 个粒子
x∈[0,1), y∈[1,2), z∈[1,2): 1 个粒子

方法2：用Pandas的分箱+分组统计 #

如果你习惯用DataFrame处理数据，这个方法会非常顺手，结果是结构化的表格，可读性拉满。

import numpy as np
import pandas as pd

testdata = np.array([[0.5,0.5,0.5],[0.6,0.6,0.6],[0.7,0.7,0.7],[1.5,0.5,0.5],[1.5,0.6,0.6],[0.5,1.5,0.5],[0.5,1.5,1.5]])

# 转换为DataFrame，方便后续处理
df = pd.DataFrame(testdata, columns=['x', 'y', 'z'])

# 定义网格边界
bins = [0, 1, 2]

# 对每个维度进行分箱，生成区间标签
df['x_bin'] = pd.cut(df['x'], bins=bins, include_lowest=True)
df['y_bin'] = pd.cut(df['y'], bins=bins, include_lowest=True)
df['z_bin'] = pd.cut(df['z'], bins=bins, include_lowest=True)

# 按三个维度的区间分组，统计每个单元格的粒子数
cell_counts = df.groupby(['x_bin', 'y_bin', 'z_bin']).size().reset_index(name='particle_count')

# 打印结构化结果
print(cell_counts)

输出的结果是一个清晰的表格，能直接看到每个三维区间对应的粒子数量。

方法3：改进np.histogramdd的使用，建立索引与区间的映射 #

如果你已经习惯用np.histogramdd，其实可以利用它返回的边界信息，把索引和实际区间对应起来，不用再瞎猜：

import numpy as np

testdata = np.array([[0.5,0.5,0.5],[0.6,0.6,0.6],[0.7,0.7,0.7],[1.5,0.5,0.5],[1.5,0.6,0.6],[0.5,1.5,0.5],[0.5,1.5,1.5]])

xcoord = testdata[:,0]
ycoord = testdata[:,1]
zcoord = testdata[:,2]
xedg = [0,1,2]
yedg = [0,1,2]
zedg = [0,1,2]

# 注意这里接收两个返回值：计数数组和各维度的边界
histo, edges = np.histogramdd([xcoord,ycoord,zcoord], bins=(xedg,yedg,zedg), range=[[0,2],[0,2],[0,2]])

# 遍历所有单元格索引，映射到对应的区间
for i in range(histo.shape[0]):
    for j in range(histo.shape[1]):
        for k in range(histo.shape[2]):
            x_interval = f"[{edges[0][i]}, {edges[0][i+1]})"
            y_interval = f"[{edges[1][j]}, {edges[1][j+1]})"
            z_interval = f"[{edges[2][k]}, {edges[2][k+1]})"
            count = histo[i,j,k]
            if count > 0:
                print(f"x{x_interval}, y{y_interval}, z{z_interval}: {int(count)} 个粒子")

这里的edges是np.histogramdd返回的边界数组，通过索引i,j,k就能直接找到对应的x/y/z区间，完美解决了索引无意义的问题。

内容的提问来源于stack exchange，提问作者Smtl