Voxelize non-uniformly gridded points

[laserman781]

I would like to voxelize non-uniformly gridded points. This is the code I have written so far:

vtkpoints = PVGeo.points_to_poly_data(unregular_blocks)
voxelizer = PVGeo.filters.VoxelizePoints()
p = pv.Plotter(notebook = False)
grid = voxelizer.apply(vtkpoints)
p.add_mesh(grid)
p.show_axes()
p.plot()

How would I go about doing this?

This is the error I am receiving:

ERROR: In C:\VPP\standalone-build\VTK-source\Filters\Python\vtkPythonAlgorithm.cxx, line 191
vtkPythonAlgorithm (0000026171FAEB30): Failure when calling method: "ProcessRequest": 
Traceback (most recent call last):
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\vtk\util\vtkAlgorithm.py", line 151, in ProcessRequest
    return vtkself.ProcessRequest(request, inInfo, outInfo)
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\vtk\util\vtkAlgorithm.py", line 197, in ProcessRequest
    return self.RequestData(request, inInfo, outInfo)
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\PVGeo\filters\voxelize.py", line 204, in RequestData
    self.__dx, self.__dy, self.__dz, grid=pdo)
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\PVGeo\filters\voxelize.py", line 110, in points_to_grid
    x,y,z = self.estimate_uniform_spacing(xo, yo, zo)
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\PVGeo\filters\voxelize.py", line 88, in estimate_uniform_spacing
    xr, yr, zr, dx, dy, angle = r.estimate_and_rotate(x, y, z)
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\PVGeo\filters\xyz.py", line 572, in estimate_and_rotate
    angle, dx, dy = self._estimate_angle_and_spacing(pts)
  File "C:\Users\User\AppData\Local\Continuum\anaconda3\lib\site-packages\PVGeo\filters\xyz.py", line 533, in _estimate_angle_and_spacing
    distall, ptsiall = tree.query(pts, k=2)
  File "sklearn\neighbors\binary_tree.pxi", line 1309, in sklearn.neighbors.kd_tree.BinaryTree.query
ValueError: k must be less than or equal to the number of training points

[banesullivan]

Can you also include the point set that you are working with? (Attach as a file here)

[laserman781]

I would prefer not to share the files, is this a problem @banesullivan ?

Thank you.

[banesullivan]

I'm not sure I can track down this issue without seeing your data - perhaps you could email them to me and I'll delete the file when finished testing

[laserman781]

@banesullivan email sent.

[banesullivan]

Thanks, @laserman781 - The VoxelizePoints filter from PVGeo is built for regularly spaced grids however it is possible to use it to create voxels of uniform size throughout a non-uniform/unstructured point cloud like yours.

I'm not sure if this is what you want, but here's a go at it:

import pyvista as pv
import PVGeo
import pandas as pd

df = pd.read_csv(filename)
point_cloud = PVGeo.points_to_poly_data(df)

voxelizer = PVGeo.filters.VoxelizePoints()
voxelizer.set_deltas(5, 5, 2) # Your block sizes in dx, dy, dz
voxelizer.set_estimate_grid(False) # This is crucial for this point cloud
grid = voxelizer.apply(point_cloud)
grid.plot(notebook=False)

[banesullivan]

From looking at your dataset, it appears to be a slice/threshold through a structured grid. If possible, I'd recommend trying to figure out where that data came from and if there is an original grid that you could use instead. If you do acquire that additional info, you could follow appoaches similar to those in the examples in #16 and #28

[laserman781]

@banesullivan In this case, each point represents a block centroid, so since the grid is non-uniform, the blocks will not be all the same sizes, some will be larger or smaller than others. But, they should all be touching... How can we accomplish this?

Thank you so much!

[banesullivan]

Do you know the block sizes? Because the VoxelizePoints filter allows you to manually set the block sizes for every block

[banesullivan]

All you would need are three arrays: the X, Y, and Z sizes for each voxel

[laserman781]

@banesullivan that is an unknown and is based on the distances between the points, the goal is to be able to plot the blocks for any grid, whether the blocks all be the same size or have several different sizes. I attempted this with a regular grid and it was able to determine automatically the recovered cell size and recovered angle but only if they are in a regular grid.

Is it possible to do the same thing with several different block sizes depending on the different block centroid distance?

[banesullivan]

I attempted this with a regular grid and it was able to determine automatically the recovered cell size and recovered angle but only if they are in a regular grid.

That's good to hear as the VoxelizePoints filter was designed for exactly that.

Is it possible to do the same thing with several different block sizes depending on the different block centroid distance?

Perhaps - but the VoxelizePoints filter is likely not the algorithm of choice as it is intedended for regular/uniform grids.

How would you compute the block centroid distance?

[banesullivan]

So I've been experimenting with this and I don't have an immediate solution - I did, however, fix a few bugs with the VoxelizePoints filter which is updated on PVGeo's master branch.

In the meantime, I think PyntCloud might have some features that could help recover the grid from the point cloud. I'd recommend checking out that library. FYI, I built some interoperability between PyntCloud and PyVista a while back, but we might need to add in a few additional features if PyntCloud can properly recreate the grid to be able to put that grid into a PyVista object

@daavoo, would you be able to chime in here?

Using code above:

from pyntcloud import PyntCloud

cloud = PyntCloud.from_pyvista(point_cloud)

# Now leverage PyntCloud on the `cloud` object

[banesullivan]

If you can find a way to compute the size sizes then all you would need to do is this (leverage version 2.0.3 of PVGeo which I just released a few minutes ago):

import pyvista as pv
import PVGeo
import pandas as pd
import numpy as np

df = pd.read_csv('sub_blocks.csv')
point_cloud = PVGeo.points_to_poly_data(df)

# Replace this with some sort of way to compute size of each voxel
x_cell_sizes = np.random.randint(3, 6, point_cloud.n_points)
y_cell_sizes = np.random.randint(3, 6, point_cloud.n_points)
z_cell_sizes = np.random.randint(3, 6, point_cloud.n_points)

voxelizer = PVGeo.filters.VoxelizePoints(unique=True, estimate=False)
voxelizer.set_deltas(x_cell_sizes, y_cell_sizes, z_cell_sizes)
grid = voxelizer.apply(point_cloud)

p = pv.Plotter(notebook=False)
p.add_mesh(grid)
p.show_grid()
p.add_axes()
p.show()

[laserman781]

Upon revising the thread, when you mentioned "A recovered connected grid with proper geometry". This is what is needed. The block point data is based on their centroids and the goal is to recover their proper geometry. They are all touching in reality and should, in general, be the same size. Although, that is not always the case and in some areas, the blocks can become smaller and in more quantity (denser) as we have seen in the point when they were displayed.
Does this make a difference in the solution to the problem? I leveraged the PVGeo voxelizer filter in order to accomplish this for a uniform grid but all I am trying to do is recover the original geometry.

Does the solution change if we want to recover the original geometry of the blocks instead of simply Voxelize them?

[banesullivan]

Does this make a difference in the solution to the problem?

Well, sort of. What changes is how we estimate the cell sizes which has to be done in the X, Y, and Z directions separately for every cell/voxel. Once you have a way to calculate what the cell sizes should be for every point, then we can use the VoxelizePoints filter in just the same way.

Computing those cell sizes is a lot easier said than done. I've tried to develop an algorithm to do this in the past with no luck - even for just rectilinear grids like this:

import pyvista as pv
import numpy as np

x = np.array([0, 2, 2.5, 2.75, 3.0, 4, 6])
y = np.array([0, 1, 1.5, 1.75, 2, 2.5, 4.5])
z = np.array([0])
grid = pv.RectilinearGrid(x, y, z)

p = pv.Plotter()
p.add_mesh(grid, show_edges=True, )
p.add_mesh(grid.cell_centers(), color='black')
p.show_grid()
p.enable_parallel_projection()
p.show(cpos='xy')

Trying to recover the original cell sizes given just the cell centers is incredibly challenging because the cell edges are not directly in the middle of two nodes. The solver would have to account for almost all the nodes on a mesh in order to make a decision about the size of a single cell.

If someone has an idea of how to implement an algorithm that can solve that problem, I will very happily code it up!

[banesullivan]

The mesh type you emailed me, @laserman781, is more like a TreeMesh or UnstructuredGrid so solving that kind of geometry is a whole other beast that I'm not sure how to approach

[laserman781]

I am currently taking a stab at an algorithm that can solve this. But if I already have the incrementation between blocks (the block sizes), will we be able to produce the result?
(See email)

[banesullivan]

Yes, just pass those cell sizes into arrays as I do in #41 (comment)

[laserman781]

After trying to pass the x,y,z cell sizes as an array to the voxelizer I keep receiving this error:
PVGeoError: X-Cell spacings are not properly defined for all points.

Code:

x_cell_sizes = df[['XINC']].values
y_cell_sizes = df[['YINC']].values
z_cell_sizes = df[['ZINC']].values

voxelizer = PVGeo.filters.VoxelizePoints()
voxelizer.set_deltas(x_cell_sizes.values , y_cell_sizes.values, z_cell_sizes.values)
voxelizer.set_estimate_grid(False) # This is crucial for this point cloud
grid = voxelizer.apply(vtkpoints)
grid.plot(notebook=False)

How can I fix this?

[banesullivan]

Make sure you have the latest version of PVGeo and don't call .values twice:

x_cell_sizes = df[['XINC']].values
y_cell_sizes = df[['YINC']].values
z_cell_sizes = df[['ZINC']].values

voxelizer = PVGeo.filters.VoxelizePoints()
voxelizer.set_deltas(x_cell_sizes , y_cell_sizes, z_cell_sizes)
voxelizer.set_estimate_grid(False) # This is crucial for this point cloud
grid = voxelizer.apply(vtkpoints)
grid.plot(notebook=False)

[laserman781]

Same error occurs for me

[banesullivan]

Is vtkpoints the proper object? - that is, does this condition hold?

vtkpoint.GetNumberOfPoints() == len(x_cell_sizes) == len(y_cell_sizes) == len(z_cell_sizes)

[laserman781]

This is the code I use to create vtkpoints:

iregular_blocks = pd.read_csv('blocks.csv')
vtkpoints = PVGeo.points_to_poly_data(iregular_blocks)

They all seem to have the same number of points although it returns False...

[banesullivan]

Something is wrong with your code - Try restarting your kernel and running in sequential order. It looks like something modified the cell sizes arrays between when you print the lengths and when you ran the conditional

[laserman781]

I figured out the problem, turns out the points data frame had had fewer points than the size data frame. It now plots exactly the way I want it to.

Thank you so much for your help!

[ammaryasirnaich]

Hi there, I have a question related to the voxelization of the 3D scene. I have a 3D scene that I want to voxelize and want to analyze the point clouds inside those voxels. The original 3D lidar point clouds look something as

and when I run the code

cloud = PyntCloud.from_file(fault_file)
vtkpoints = PVGeo.points_to_poly_data(cloud.points)
voxelizer = PVGeo.filters.VoxelizePoints()
voxelizer.set_deltas(0.4, 0.2, 0.2) # Your block sizes in dx, dy, dz
voxelizer.set_estimate_grid(False) 
grid = voxelizer.apply(vtkpoints)
grid.plot(notebook=False)

Its plots the voxels that are solid gray

Is there any way that I can plot voxel with boundaries as being transparent to see point clouds inside them?

[RichardScottOZ]

Use opacity=0.5 or whatever suits in plots

[ammaryasirnaich]

HI @RichardScottOZ. thanks for the tip. I tried it with opaciy=05 and 0.1 but the plot makes the solid face of the voxel transparent but the point clouds occupied inside the voxel are still invisible

[RichardScottOZ]

You are saying each voxel has something inside it?

[ammaryasirnaich]

Yes, you can see from the first image that i have plotted. As these are the point clouds with intensity values. After that, I have applied a grid structure for generating voxels in the 3D scene.

[RichardScottOZ]

#482 might be of interest?

[RichardScottOZ]

Still not sure I see, what does one of your voxels and cloud inside look like?

[ammaryasirnaich]

Sorry for not being clear, actually i am trying to be something like this

so far i have just managed to plot the point clouds which look like

Not i want to apply voxel upon the above point cloud similarly as the first image i have shown.

[RichardScottOZ]

Thanks - so the plot above is just the point clouds?

[RichardScottOZ]

Now possibly - voxelise them as above

Add that to a plot (with transparent opacity)
them add the original point clouds mesh to the same plot?

[ammaryasirnaich]

Thanks - so the plot above is just the point clouds?

Yes absolutely

[ammaryasirnaich]

Now possibly - voxelise them as above

Add that to a plot (with transparent opacity) them add the original point clouds mesh to the same plot?

Thanks for the tip really appreciate that. I will give it a try as you suggested

[banesullivan]

Now possibly - voxelise them as above

Add that to a plot (with transparent opacity)
them add the original point clouds mesh to the same plot?

Yep - specifically, you may want to use style='wireframe' when plotting the voxels: https://docs.pyvista.org/api/plotting/_autosummary/pyvista.Plotter.add_mesh.html#add-mesh