Efficient method for creation of large point-clouds and meshes from numpy arrays

[kleinsimon]

Hi there,
first of all, I have to thank you for providing this package, it really provides a very nice interface to occ.
In my project, I have to display point clouds and meshes with about 10M vertices. While the display of these entities works nice, the creation is very slow. I'm using the following approach:

occ_points = Graphic3d_ArrayOfPoints(self.points.shape[0])
for p in self.points:
    occ_points.AddVertex(*p)
pc = AIS_PointCloud()
pc.SetPoints(occ_points)

and for the mesh:

n_verts = len(vertices)
c_pts = TColgp_Array1OfPnt(1, n_verts)

for i in range(n_verts):
    c_pts.SetValue(i+1, gp_Pnt(*vertices[i,:]))

n_faces = len(faces)
c_faces = Poly_Array1OfTriangle(1, n_faces)

for i in range(n_faces):
    c_faces.SetValue(i+1, Poly_Triangle(int(faces[i, 0])+1, int(faces[i, 1])+1, int(faces[i, 2])+1))

occ_tris = Poly_Triangulation(c_pts, c_faces)

While this works, the loops are pretty slow. I imagine, there are better ways to send the arrays through swig. I think numpy provides swig extensions for such cases, but I think it's not used here.
On the other side, since the arrays on C++-side are simply structs, it would also be possible to simply point it to the numpy data in memory. As far as I can see, the NCollection_Array1 provides three options for that:

• Initialize the length and provide a pointer to the first element
• Use the "Move" method of an already initialized object to point it to some memory
• Use the "Assign" method to copy the data from a pointer

I can get a PyCapsule object for the pointer from a numpy-array using the array_struct property, but this is not accepted as a type by the swig interface.

Are there other approaches for reducing data copy and loops like that? Has anybody successfully transferred data using a numpy pointer or other techniques?

Thanks in advance

[tpaviot]

I agree that the creation of large meshes using python loops and pythonocc is quite long. Actually, the current python wrapper is not optimized for performances: every call to an occ is put into a try statement to catch occ exceptions. If you have one million calls to the same method, then it might slows down the process.

The idea to use numpy and to pass/get numpy arrays to/from occ is interesting. We talked about that many years ago, but no one, as far as I know, never did it. I'm not familiar with PyCapsule.

Note that I recently pushed an addition (see #992) that seems to be related, you could use this approach. See https://github.com/tpaviot/pythonocc-demos/blob/master/examples/core_mesh_data_source.py and tpaviot/pythonocc-demos@2061eed You may have a look at it, perhaps it does the job without any further effort

[kleinsimon]

Thanks for the hint. I stumbled across this demo, but by now I only used the vanilla builds from conda. But I think here the data loading is not done in python but in the stl-reader.

I managed to get a decent 10x speedup using cython and interacting with the direct interface calls:

import numpy as np
cimport numpy as np
DTYPEF = np.float
DTYPEI = np.int

ctypedef np.float_t DTYPE_F_t
ctypedef np.int_t DTYPE_I_t

from cython.operator cimport dereference as deref

from OCC.Core._gp import new_gp_Pnt
from OCC.Core._TColgp import new_TColgp_Array1OfPnt, TColgp_Array1OfPnt_SetValue
from OCC.Core._Poly import \
     new_Poly_Array1OfTriangle, \
     Poly_Array1OfTriangle_SetValue, \
     new_Poly_Triangle, \
     new_Poly_Triangulation
     
from OCC.Core.Poly import Poly_Triangulation


def create_triangulation(np.ndarray[DTYPE_F_t, ndim=2] vertices, np.ndarray[DTYPE_I_t, ndim=2] faces):
    cdef int n_pts = vertices.shape[0]
    p_array = new_TColgp_Array1OfPnt(0, n_pts-1)
    
    for i in range(n_pts):
        TColgp_Array1OfPnt_SetValue(p_array, i, new_gp_Pnt(vertices[i,0],vertices[i,1],vertices[i,2]))
        
    cdef int n_faces = faces.shape[0]
    f_array = new_Poly_Array1OfTriangle(0, n_faces-1)
    
    for i in range(n_faces):
        Poly_Array1OfTriangle_SetValue(f_array, i, new_Poly_Triangle(faces[i,0], faces[i,1], faces[i,2]))
        
        
    return Poly_Triangulation(p_array, f_array)
        

This is nice, because I have to do nothing else but create the pyx file. I'm rather new to cython and think with better typing speed may increase further. Found something related here:
https://cython-users.narkive.com/8rtYS8bI/using-cython-with-swig-arguments

update: the max-bound was wrong, now fixed to be length-1

[tpaviot]

Yes, in the demo data is loaded through the STL file, but you can also add your own data, I will upload another demo.

The cython/SWIG approach is very interesting, thank you, I did not know it was possible to do that. Does it require some compilation stage?

[tpaviot]

It might also be possible to use numpy together with occt with the same approach

[kleinsimon]

Yes, it must be compiled. I had to install cython, as well as some c++ build tools from VS (for windows). Then the setup.py has to be called like this python .\setup.py build_ext --inplace. That way, the compiled pyd file is created (other types on Linux I think), which is more or less a dll.
In python, its very easy to use:

from poly import create_triangulation

occ_tris = create_triangulation(vertices, faces)

which gives the Poly_Triangulation object back.

I attached the sources as well as my compiled file for python 3.8 and windows 64.

occ_wrapper.zip
(EDIT: Uploaded a defect file... sorry)

[kleinsimon]

I think there is a lot to optimize, I could not manage to type the swig objects correctly, so python objects will still be generated. Also there seems to be a newer method for accessing the numpy arrays. I will look further into this.

[kleinsimon]

Maybe this is not possible without diving into the swig internals... it seems that swig only exposes calls to method, which create or await SwigPyObject objects, which are python objects that wrap the pointer to the c object. There seems to be no way to circumvent the creation of these objects in cython.

So maybe the easier method would be to allow your MeshDataSource to be generated with numpy arrays for the faces and vertices directly. From what I can see, this seems not to be to complicated. I can not tell, if new dependecies are necessary. It seems, a .i file from numpy is required, as well as some tweaks to the typemap.

https://numpy.org/doc/stable/reference/swig.interface-file.html#using-numpy-i

I will see if I can manage to get the build working and implement this.

[kleinsimon]

So I got this working using numpy.i and some modifications to MeshDataSouce with a huge speedup. But by now, I can not get OCC to draw the mesh in a shaded view, I always get a wireframe with marks at the vertices, regardless what flags I set in the builder or what display mode I set in the MeshVS object. This is contrary to AIS_Triangulation, where I could only get a shaded view, but could not set color or other aspects.
For compilation to work, numpy has to be installed in the python env and its include directory has to be in path. I think a dependency to numpy is generally no problem, since its installed nearly by default, but maybe there is a way to make this optionally. As soon as I get the mesh to display without errors and I'm conviced theres nothing wrong in the data, I could make me pull request

The numpy.i file is to be assumed in src\SWIG_files\common

MeshDataSource Numpy.zip

[tpaviot]

Great, please open a PR. After I checked it works properly, I will certainly add this feature as a compilation option, since I do not want to force an additional dependency to be installed by pthonocc users. BTW, did you try numba (http://numba.pydata.org/)?

[kleinsimon]

jeah, I'm using Numba regulary, but it seems it cannot work with the swig libraries.

[kleinsimon]

I just noticed, that the example core_mesh_data_source produces a wireframe output without shaded faces, too. Is this inteded? Or is something still missing in the datasource? Must Get3DGeom be implemented for this? Trying to figure it out using the docs...

[tpaviot]

Fixed