voxel_down_sample_and_trace not return the full index of points

[faultaddr]

Describe the bug I want to get the voxelized point cloud with its corresponding index in original point cloud. but i found that the total index count does not equals the original point cloud number. in the case below,it expected to be 40000,but actually i got 6519.which means i cannot get each voxel's original coordinate.maybe the trace function become useless.

To Reproduce Steps to reproduce the behavior: voxel_down_sample_and_trace(pcd,voxel_size,min_bound,max_bound)

Expected behavior get the downsampled pointcloud and its corresponding index in original point cloud

Screenshots

Environment (please complete the following information):

• OS: Manjaro 18.04
• Python version:3.6.9
• Open3D version:0.8.0
• Is this remote workstation?: no
• How did you install Open3D?: build_from_source ( gcc)

[os-gabe]

I am also confused by this

[HM102]

This is not a bug, you will not get len(new_index) = len(pcd.points) for sure. The new_index represent the original point id of the higher resolution point cloud, but for the points in the new_pcd only, which is definitely less points then pcd, since it is downsampled!

Now the question is, why do you get more then one index, I guess these are the index of points that contributed to the downsampled point in the new_pcd! @syncle this looks like a reoccurring question, can you clarify, if that's the case.

[faultaddr]

okay,I know what u mean,maybe the function just not the one that could solve my problem. well,I expecte to get the all point index in a small grid after downsample. so that i could find the one to many mapping:from a downsampled point to original points.(n-->1 then trace 1-->n) thank u any way! best wishes!

[os-gabe]

Lets say I have an original pointcloud with N points in it. When I run voxel_down_sample_and_trace on the original pointcloud I get back a new pointcloud with M points in it where M <= N. I also get back the trace matrix which is a mapping from the downsampled cloud to the original cloud. This all makes sense so far.

In my case the trace matrix is always Mx8. The fact that it has M rows makes sense becaust it is a mapping from the downsampled pointcloud but why does it have 8 columns? Is it assuming at most 8 points from the original cloud contribute to each point in the new cloud? If so, that seems like a bad assumption.

Here is an example to produce some data we can consider concretely:

N = 1000
vox_size = 0.5

pts = np.zeros((N, 3))
pts[:, 0] = np.linspace(0, 1, N)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(pts)
pcd_ds, trace = pcd.voxel_down_sample_and_trace(
    vox_size, pcd.get_min_bound(), pcd.get_max_bound(), False
)
# Offset the new points a little in y so we can see them in the visualization
pcd_ds.points = o3d.utility.Vector3dVector(np.asarray(pcd_ds.points) + [0.0, 0.1, 0.0])
print("N = %d" % (len(pcd.points)))
print("M = %d" % (len(pcd_ds.points)))
print("Trace Shape = ", trace.shape)
print(trace)

o3d.visualization.draw_geometries([pcd, pcd_ds])

If you run the example you can see that you get the following trace matrix:

[[999  -1  -1  -1  -1  -1  -1  -1]
 [249 499  -1  -1  -1  -1  -1  -1]
 [749 998  -1  -1  -1  -1  -1  -1]]

Am I to interpret this as saying point 0 in the new cloud came from point 999 in the old cloud, point 1 in the new cloud came from points 249 and 499 in the old cloud, and point 2 in the new cloud came from points 749 and 998 in the old cloud?

This interpretation appears to be incorrect to me so what am I missing?

[eddienewton]

@os-gabe looking at the code in function VoxelDownSampleAndTrace, the trace matrix always has 8 rows because the trace represents a cube (8 vertices).

From what I can discern, the column chosen is dependent on the closest cubic vertex to that point.

Regarding getting the mapped voxel index, you'll probably just need to implement a python equivalent to (found in DownSample.cpp):

auto ref_coord = (points_[i] - voxel_min_bound) / voxel_size;
auto voxel_index = Eigen::Vector3i(int(floor(ref_coord(0))),
    int(floor(ref_coord(1))),
    int(floor(ref_coord(2))));

(Edited): I also just noticed that VoxelDownSampleAndTrace and VoxelDownSample have different index implementations. You'll need to adjust your index calculation accordingly.

Hope this helps

[os-gabe]

Thanks @eddienewton for your help. For what it's worth here is a python implementation I came up with to do what I needed:

def voxel_filter(pts, scores, grid_size, use_3d=False, min_num_pts=4):
    mins = pts.min(axis=0) - grid_size
    maxs = pts.max(axis=0) + grid_size
    bins = [np.arange(mins[i], maxs[i], grid_size) for i in range(len(mins))]

    si = 2
    if use_3d:
        si = 3

    counts, edges, binnumbers = stats.binned_statistic_dd(
        pts[:, :si],
        values=None,
        statistic="count",
        bins=bins[:si],
        range=None,
        expand_binnumbers=False
    )

    ub = np.unique(binnumbers)
    pts_ds = []
    scores_ds = []
    for b in ub:
        if len(np.where(binnumbers == b)[0]) >= min_num_pts:
            pts_ds.append(pts[np.where(binnumbers == b)[0]].mean(axis=0))
            scores_ds.append(scores[np.where(binnumbers == b)[0]].mean())
    pts_ds = np.vstack(pts_ds)
    scores_ds = np.vstack(scores_ds).reshape(-1)
    return pts_ds, scores_ds

I wanted a voxel filter that found the mean position of all the points within each voxel and also the mean score (where pts.shape = (n, 3) and scores.shape = (n, )). I'm assuming my implementation will not be as fast as VoxelDownSampleAndTrace but it got me what I needed and was relatively easy to understand.

[eddienewton]

@os-gabe thanks for this! Looks pretty useful actually. I might have to borrow it :)

[CuberrChen]

感谢@eddienewton的帮助。对于这里值得的是我想出的Python实现来做我需要做的事情：

def voxel_filter(pts, scores, grid_size, use_3d=False, min_num_pts=4):
    mins = pts.min(axis=0) - grid_size
    maxs = pts.max(axis=0) + grid_size
    bins = [np.arange(mins[i], maxs[i], grid_size) for i in range(len(mins))]

    si = 2
    if use_3d:
        si = 3

    counts, edges, binnumbers = stats.binned_statistic_dd(
        pts[:, :si],
        values=None,
        statistic="count",
        bins=bins[:si],
        range=None,
        expand_binnumbers=False
    )

    ub = np.unique(binnumbers)
    pts_ds = []
    scores_ds = []
    for b in ub:
        if len(np.where(binnumbers == b)[0]) >= min_num_pts:
            pts_ds.append(pts[np.where(binnumbers == b)[0]].mean(axis=0))
            scores_ds.append(scores[np.where(binnumbers == b)[0]].mean())
    pts_ds = np.vstack(pts_ds)
    scores_ds = np.vstack(scores_ds).reshape(-1)
    return pts_ds, scores_ds

我想要一个体素滤波器，该滤波器可以找到每个体素内所有点的平均位置以及均值score（where pts.shape = (n, 3)和scores.shape = (n, )）。我以为我的实现将不会很快，VoxelDownSampleAndTrace但是它满足了我的需求，并且相对容易理解。

can tou tell me the 'stats' is which package? and how to use this function?

[os-gabe]

Sorry about that: from scipy import stats

Do you have an example of what you would like to use it for?

[CuberrChen]

Sorry about that: from scipy import stats

Do you have an example of what you would like to use it for?

Thank you for your quick reply！yep..I want to input (x,y,z) pointcloud which is numpy array to get downsampled result. but i can't use your function directly.the example of my usage is as follows: points = np.random.rand(10000, 3) down_points = voxel_filter(points,0.05,use_3d=True)

Can you direct me？I didn't need scores, so I deleted that part。the result is as follow：

Traceback (most recent call last): File "/home/chenxiangbo/fast/hahhahahha.py", line 44, in <module> down_points = voxel_filter(points,0.05,use_3d=True) File "/home/chenxiangbo/fast/hahhahahha.py", line 23, in voxel_filter expand_binnumbers=False File "/home/chenxiangbo/.local/lib/python3.7/site-packages/scipy/stats/_binned_statistic.py", line 518, in binned_statistic_dd if not np.isfinite(values).all() or not np.isfinite(sample).all: TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

[CuberrChen]

Sorry about that: from scipy import stats

Do you have an example of what you would like to use it for? Thank you very much for your help. I have solved this problem! Although i did't use your way，your work is great. my solve：

points = np.random.rand(10000, 3)
point_cloud = open3d.PointCloud()
point_cloud.points = open3d.Vector3dVector(points)
open3d.visualization.draw_geometries([point_cloud])
down_sample_point_cloud = open3d.geometry.voxel_down_sample(point_cloud,voxel_size = 0.07)
open3d.visualization.draw_geometries([down_sample_point_cloud])
down_sample_point_cloud = np.asarray(down_sample_point_cloud.points)

[germanros1987]

@griegler do you think this is something you could look into?

[Abdurrasith]

how about open3d.geometry.voxel_down_sample_and_trace with the input open3d.io.read.point_cloud? it seems always getting AttributeError: module 'open3d.cpu.pybind.geometry' has no attribute 'voxel_down_sample_and_trace'

how to solve it? thanks

[aminullah6264]

@Abdurrasith

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(coord)
pcd.colors = o3d.utility.Vector3dVector(color)   
pcd, iddx = pcd.voxel_down_sample_and_trace(g_s, pcd.get_min_bound(), pcd.get_max_bound(), False)  
coord_down = np.array(pcd.points)
color_down = np.array(pcd.colors)

[matekosztricz]

Hi,

I was wondering If there has been any updates regarding this issue. I'm trying to count points inside a VoxelGrid's Voxels. with little success. After a long debugging session and reading this thread my issue is most likely with the different indexing in the voxel_down_sample_and_trace() method and the VoxelGrid.create_from_point_cloud_within_bounds() method. My general idea was that If I use these 2 methods with the same given bounding box, it would result in the same VoxelGrid being inspected, but that is incorrect. Does anyone have a solution for this? Is there ever going to be an implemented way to solve this in an easier way?

Here's reproduction of what I'm trying to do, with a funky visualization that aims to show that there is a problem (this is the only way I could easily show Voxels one-by-one, since adding voxels to an empty grid did not show anything upon displaying)

import numpy as np
import open3d as o3d
import copy

# Define the parameters
num_points = 200  # Number of random points
cube_size = 5  # Size of the cube
voxel_size = 0.5

voxel_points_index_array = []
voxel_point_count_array = []

# Generate random points for the first point cloud
points_cloud1 = cube_size * np.random.rand(num_points, 3)
pcd1 = o3d.geometry.PointCloud()
pcd1.points = o3d.utility.Vector3dVector(points_cloud1)

min_bound = pcd1.get_min_bound()
max_bound = pcd1.get_max_bound()

voxel_grid1 = o3d.geometry.VoxelGrid.create_from_point_cloud_within_bounds(pcd1, voxel_size=voxel_size,
                                                                               min_bound=min_bound,
                                                                               max_bound=max_bound)
pcd_copy1 = copy.deepcopy(pcd1)
_, _, points_index_list = pcd_copy1.voxel_down_sample_and_trace(voxel_size=voxel_size, min_bound=min_bound,
                                                               max_bound=max_bound)

for i, idx in enumerate(points_index_list):
    voxel_points_index_array.append(np.asarray(idx))
    voxel_point_count_array.append(len(idx))

voxs = voxel_grid1.get_voxels()
coord_frame = o3d.geometry.TriangleMesh.create_coordinate_frame(size=1.0, origin=[0, 0, 0])

for n in range(10):
    vg = copy.deepcopy(voxel_grid1)
    for i in range(n, len(voxs), 1):
        vg.remove_voxel(voxs[i].grid_index)
    print(points_index_list[n])
    print(voxel_point_count_array[n])
    pcd2 = pcd1.select_by_index(points_index_list[n])
    o3d.visualization.draw_geometries([vg, pcd2])