Extract and resample point cloud to a given size from a binary mask

[sohmandal]

Hi,

From a given binary mask of a 3D volume containing one object, does vedo support any pipeline, which can used to extract the surface point cloud (of length M) of this object at then resample this point cloud to a desired size (of length N), where N>M?

I have attached one such binary mask. mask.zip

[marcomusy]

hi, what about this

from vedo import *

arr = np.load("Downloads/mask.npy")
vol = Volume(arr)

iso = vol.isosurface().lw(0.1)
iso.smooth().subdivide(2)
iso.decimate(N=2000)
print(iso.points().shape)

show(iso, axes=1).close()

[sohmandal]

Hi @marcomusy, thank you very much for your reply.

The mask that is being considered, is of shape [64,128,128]. Now from running the following snippet

arr = np.load("mask.npy")
vol = Volume(arr)
iso = vol.isosurface().lw(0.1)

print('surface point cloud coordinates on axis 0-', np.sort(iso.points()[:,0]))

the output is surface point cloud coordinates on axis 0- [70.333336 70.333336 70.333336 ... 91.666664 91.666664 91.666664]

Since the length of the mask at axis 0 is 64, shouldn't the above output be bounded between [0,63] ? Or am I missing something obvious?

[marcomusy]

Hi, no, because the isosurface is a polygonal object and it's bounds do not correspond to the volume bounds. When creating the Volume object you can specify the spacing keyword, i.e. the 3 sizes of the voxel, which are 1 by default.

[marcomusy]

PS, sorry maybe this is what you meant?

from vedo import *

arr = np.load("/home/musy/Downloads/mask.npy")
# arr = np.transpose(arr, axes=[2, 1, 0])
print(arr.shape)
vol = Volume(arr)
print(vol.bounds())

iso = vol.isosurface().lw(0.1)
iso.smooth().subdivide(2)
iso.decimate(N=2000)
print(iso.points().shape)

show(iso, vol.box(), axes=1).close()

numpy and vtk have a different convention so you may need to transpose the array.

[sohmandal]

Hi, that is it. Perfect, thank you very much! I missed the difference in ordering convention between numpy and vtk.

I have couple more questions-

1. Why do you set the width of mesh edges at a particular value (lw(0.1)) after you get the isosurface?
2. Is it necessary to smooth the surface before doing the subdivision? Because I am assuming the smoothing could eliminate certain fine features of the surface.
3. Is the decimate(N = k) supposed to reduce the surface to k+1, not to k points as given as the parameter to the decimate() method?
4. I am intending to perform this pipeline on (instance) masks containing multiple objects on-the-fly for training a neural network. Since this is an expensive pipeline in terms of the time that it takes, would you suggest any performance speedup in any way over this pipeline?

[marcomusy]

Please find a new version where I fixed the transpose problem: pip install -U git+https://github.com/marcomusy/vedo.git now you can remove the transpose bit.

1. that was just to visualize the mesh edges
2. yes, you need to play around with smoothing (if necessary at all) and decimate until you get the desirable result (there is no objective "right or wrong" here..)
3. I never tested it, I think the decimate may only generate about the desired number of vertices depending on the specific problem
4. Uhm if you think the bottleneck is in this part of the pipeline, you may strip down the vedo wrapping to the native vtk and you may get some marginal speedup. Smoothing and subdividing can be the most expensive passes.

[sohmandal]

Fantastic, thank you for the update!

3. For all the instances that I tried, it seems like decimate(N = k) reduces the surface to (k + 1) points
4. Thank you for the suggestions.

[sohmandal]

Hi @marcomusy, is there a way to make decimate output a specific value? I am noticing for decimate(N = k) reduces the surface to either (k + 1) points or k points from case to case basis.

[marcomusy]

Hi Soham, sorry I'm afraid you can only control the behavior of decimate(N=...) only up to +-1 precision. You may give it a try with decimate(fraction=...) where you compute the fraction based on the nr of points.

[sohmandal]

Thanks @marcomusy! Seems like it would still be uncertain. However for my use case, I could accommodate +-1 points in the point cloud.

On this note, I was wondering - how does decimate remove points from the point cloud? I assume this happens in some uniform manner?

[marcomusy]

how does decimate remove points from the point cloud? I assume this happens in some uniform manner?

from https://vtk.org/doc/nightly/html/classvtkQuadricDecimation.html

The algorithm is based on repeated edge collapses until the requested mesh reduction is achieved. Edges are placed in a priority queue based on the "cost" to delete the edge. The cost is an approximate measure of error (distance to the original surface)–described by the so-called quadric error measure. The quadric error measure is associated with each vertex of the mesh and represents a matrix of planes incident on that vertex. The distance of the planes to the vertex is the error in the position of the vertex (originally the vertex error iz zero). As edges are deleted, the quadric error measure associated with the two end points of the edge are summed (this combines the plane equations) and an optimal collapse point can be computed. Edges connected to the collapse point are then reinserted into the queue after computing the new cost to delete them. The process continues until the desired reduction level is reached or topological constraints prevent further reduction. Note that this basic algorithm can be extended to higher dimensions by taking into account variation in attributes (i.e., scalars, vectors, and so on).

This paper is based on the work of Garland and Heckbert who first presented the quadric error measure at Siggraph '97 "Surface Simplification Using Quadric Error Metrics". For details of the algorithm Michael Garland's Ph.D. thesis is also recommended. Hughues Hoppe's Vis '99 paper, "New Quadric Metric for Simplifying Meshes with Appearance Attributes" is also a good take on the subject especially as it pertains to the error metric applied to attributes.

---

Note that you can also use .subsample() directly on your pointcloud to make it very spatially uniform (for large clouds this might become slow).

[sohmandal]

Thanks @marcomusy! This is very helpful.

Would the extracted point cloud from the 3D volume have any particular consistent orientation in terms of ordering in the point cloud (such as top-to-bottom-in-clockwise, or anything similar)? If not, is there any way to enforce some sort of consistent ordering?

Further, would doing all these operations (subdivide, decimate) on the extracted point cloud hamper the ordering of the points in the resultant point cloud, if there exists any ordering?

[marcomusy]

hi

extracted point cloud from the 3D volume

i'm sure I understand this... you can isosurface a volume and obtain a polygonal mesh, this mesh has a defined orientation

if not, is there any way to enforce some sort of consistent ordering?

you can compute the normals from a pure pointcloud (no faces) with computeNormalsWithPCA

would doing all these operations (subdivide, decimate) on the extracted point cloud hamper the ordering of the points in the resultant point cloud, if there exists any ordering?

in principle such operations should retain the original information, or interpolate it.