基于LSH的3D点云邻域计数代码异常问题排查

普通哈希代码计数始终为1的Bug分析与修复 #

核心问题 #

• 哈希键未量化：直接用原始浮点坐标作为哈希键，由于浮点数精度差异，每个点的键都是唯一的，哈希桶中只有当前点自己，所以计数始终为1。
• 未遍历邻域网格桶：仅查询当前点所在的哈希桶，没有检查周围3x3x3的相邻网格单元，无法找到真正的邻点。

修复示例代码 #

truncation_dist = 0.5  # 你的截断距离

# 将坐标量化到网格单元，生成整数哈希键
def get_grid_key(point):
    return (
        int(point[0] / truncation_dist),
        int(point[1] / truncation_dist),
        int(point[2] / truncation_dist)
    )

# 构建哈希表
hash_map = {}
for point in points:
    key = get_grid_key(point)
    hash_map.setdefault(key, []).append(point)

# 统计每个点的邻域数量
counts = []
for point in points:
    current_key = get_grid_key(point)
    count = 0
    # 遍历当前网格周围3x3x3的所有邻域网格
    for dx in (-1, 0, 1):
        for dy in (-1, 0, 1):
            for dz in (-1, 0, 1):
                neighbor_key = (current_key[0]+dx, current_key[1]+dy, current_key[2]+dz)
                if neighbor_key not in hash_map:
                    continue
                # 检查邻域桶内的点是否真的在截断距离内，排除自身
                for neighbor in hash_map[neighbor_key]:
                    dist_sq = (point[0]-neighbor[0])**2 + (point[1]-neighbor[1])**2 + (point[2]-neighbor[2])**2
                    if dist_sq <= truncation_dist**2 and dist_sq > 1e-12:
                        count += 1
    counts.append(count)

Morton哈希代码未达预期的Bug分析与修复 #

核心问题 #

• 编码前未量化坐标：直接对原始浮点坐标编码，没有按截断距离量化到整数网格，导致哈希键无法正确分组邻域点。
• 邻域查找逻辑错误：误以为Morton码的相邻数值对应邻域网格，实际上Morton码是二进制位交错生成的，相邻网格的Morton码不是连续值，直接加减会找错桶。
• 编码/解码实现错误：常见位移位数错误，或未处理负坐标（Morton编码仅支持非负整数，需先将坐标平移到非负区间）。

修复示例代码 #

truncation_dist = 0.5  # 你的截断距离

# 正确的Morton编码/解码实现
def morton_encode(x, y, z):
    def spread_bits(val):
        val = (val | (val << 16)) & 0x030000FF
        val = (val | (val << 8)) & 0x0300F00F
        val = (val | (val << 4)) & 0x030C30C3
        val = (val | (val << 2)) & 0x09249249
        return val
    return spread_bits(x) | (spread_bits(y) << 1) | (spread_bits(z) << 2)

def morton_decode(code):
    def compact_bits(val):
        val = val & 0x09249249
        val = (val | (val >> 2)) & 0x030C30C3
        val = (val | (val >> 4)) & 0x0300F00F
        val = (val | (val >> 8)) & 0x030000FF
        val = (val | (val >> 16)) & 0x000003FF
        return val
    return (compact_bits(code), compact_bits(code >> 1), compact_bits(code >> 2))

# 量化坐标并平移到非负区间（处理负坐标）
min_x = min(p[0] for p in points)
min_y = min(p[1] for p in points)
min_z = min(p[2] for p in points)

def get_quantized_point(point):
    return (
        int((point[0] - min_x) / truncation_dist),
        int((point[1] - min_y) / truncation_dist),
        int((point[2] - min_z) / truncation_dist)
    )

# 构建Morton哈希表
hash_map = {}
for point in points:
    q_x, q_y, q_z = get_quantized_point(point)
    key = morton_encode(q_x, q_y, q_z)
    hash_map.setdefault(key, []).append(point)

# 统计邻域数量
counts = []
for point in points:
    q_x, q_y, q_z = get_quantized_point(point)
    count = 0
    # 遍历3x3x3邻域的量化坐标，重新编码后查找哈希桶
    for dx in (-1, 0, 1):
        for dy in (-1, 0, 1):
            for dz in (-1, 0, 1):
                n_qx, n_qy, n_qz = q_x+dx, q_y+dy, q_z+dz
                if n_qx < 0 or n_qy < 0 or n_qz < 0:
                    continue
                neighbor_key = morton_encode(n_qx, n_qy, n_qz)
                if neighbor_key not in hash_map:
                    continue
                # 验证实际距离，排除自身
                for neighbor in hash_map[neighbor_key]:
                    dist_sq = (point[0]-neighbor[0])**2 + (point[1]-neighbor[1])**2 + (point[2]-neighbor[2])**2
                    if dist_sq <= truncation_dist**2 and dist_sq > 1e-12:
                        count += 1
    counts.append(count)

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