Current approach to 3D Printing using Trimesh

The default mesh slicing option available in PySLM, Part.getVectorSlice was never particularly aimed at being high throughput at the beginning. It aimed to simply offer the functionality to test hatching approaches, that utilised existing approaches available within Trimesh. The built-in slicing option built-upon using trimesh.intersections.mesh_plane functionality existing already and was used for this purpose. This slicing option has been robust, and the author has considered edge-cases that may occur during the intersection between the slicing plane and triangular meshes. This includes additionally shell and to a limited extent processing non-manifold meshes. Additionally, an arbitrary choice of slicing plane may be used, which for 3D graphics is useful, however, for 3D printing, typically nearly all processes work layer-wise in the Z direction.

In the Trimesh implementation, a few optimisation have been built-in. When slicing multiple planes across a mesh using trimesh.intersections.mesh_multiplane, the dot product is pre-calculated between the mesh vertices and the plane normal and origin may be cached and re-used for subsequent computations across multiple slices.

Existing Accelerated Slicing Approaches

Despite the relative time spent during slicing is likely less than the majority of spent within the final hatching process that fills the part with hatch vectors. Nevertheless, the performance wasn’t adequate especially for larger triangular meshes. Some direct approaches could have been low-level optimisation using specialised Python libraries (e.g. Numba, Pyston) or re-coding in underlying c++ code. A better approach was considered.

Building-upon a recent research paper by Minetto et al. [1] , it is possible to improve the performance by selectively pre-sorting and filter triangles of the mesh the to reduce the overall complexity to O(n+m+k), where n is the number of triangles, m is the number of slices and k is the number of intersections. Their method determines the coverage each triangle and separates these to an individual list of triangles for each slice layer to isolate potential intersections in order to reduce the search space size. Inadvertently, the slicing operation becomes trivially parallel to process.

The second focus on their paper is using hash structure to more efficiently combine the generated line segments/paths from the slicing operation into correctly orientated paths. This process ensures a consistent winding number for clipping operations later. Their approach relies on rounding the coordinates within a tolerance, to identify coincident points. However, is reliant on water-tight meshes. Another author applies a variation, however, uses Graph instead of a Hash table with a modest improvement [2]. Curiously, I noticed a big discrepancy between their reported performance. For PySLM, the polygon construction is internally handled automatically using the ClipperLib or Shapely polygon processing libraries.

The key thing taken from that paper and previous research was that by pre-sorting the triangle mesh based on predefined ‘constant’ layer thickness, potential intersecting triangle candidates within the mesh can be extracting during the computation and interestingly parallelised.

Methodology

Provided mesh, the approach first extracts the minimum and maximum Z positions per triangle, which correspond with the bottom and top most layer that a potential mesh intersection occurs. This corresponds to Algorithm 2 in the previous cited paper. Across the generated list of layer slices, a binary search (numpy.searchsorted) across a pre-sorted list is used for identifying these locations. Across each triangle, their array index is appended to each layer within an inefficient for loop. Overall, the cost of this pre-sorting is only expensive when assembling the indices across the layer. Afterwards, these may then be used during the slicing stage.

# Generate a list of slices for each layer across the entire mesh
zSlices = np.linspace(-5, 5, 200)
k = 1000 # number of slices
zMin, zMax = myMesh.bounds[:,2]
zBox = np.linspace(zMin, zMax, k)


tris = myMesh.triangles

# Obtain the min and maximum Z values across the entire mesh
zVals = tris[:, :, 2]
triMin = np.min(zVals, axis=1)
triMax = np.max(zVals, axis=1)


# Below is a manual approach to sorting and collecting the triangles across each layer 

if False:

    triSortIdx = np.argsort(triMinMax[1,:])
    triMinMaxSort = triMinMax[:,triSortIdx]

    startTime = time.time()

    sortTris = []

    iSects2 = []
    for i in range(len(zBox)):
        minInside = zBox[i].reshape(-1,1) > triMinMax[0, :].reshape(1, -1)
        maxInside = zBox[i].reshape(-1,1) < triMinMax[1, :].reshape(1, -1)
        iSects2.append(np.argwhere((minInside & maxInside).ravel()))

    print('endTime array base', time.time() - startTime)

# An alteratnvie more compact way is to use binary search operator available within
# numpy.searchsorted. This locates the bottom and top layer position 

minIdx = np.searchsorted(zBox, triMin, side='left')
maxIdx = np.searchsorted(zBox, triMax, side='left')

# Attach the corresponding presorted triangles into the 
iSects = [[] for i in range(len(zBox))]

# The iterative part for assigning potential triangles for intersection on a triangle are performed here
# Note: Process is very inefficient in native Python O(n*k) 
for i in range(len(minIdx)):

    startLayer = minIdx[i]
    endLayer = maxIdx[i]
    for layer in iSects[startLayer:endLayer+1]:
        layer.append(i)

Plane-Triangle Intersections

The pre-processing has been completed, now the final slincg may be complete. Performing some micro-optimisations, further refactoring may be done to adapt the code previously present in Trimesh.

The process of slicing or intersection, is simple. Suprisingly, there are a no obvious references for this process. Slicing a triangular mesh, relies firstly computing the potential intersection of each edge of a faceted mesh, in this case it is generalised for an abritrary plane. Firstly, the vector between the mesh vertices (v_i) and the plane origin (v_o) is calculated – red lines in the diagram. The dot product is taken with the slicing plane normal \textbf{n}. The sign of the dot product indicates if the point lies above or below the plane – zero uniquely is the intersection.

plane_origin = np.array([0,0,0])
plane_normal = np.array([0,0,1])

vertex_dots = np.dot(myMesh.vertices - plane_origin, plane_normal)

This is extremely convenient to compute, and it only requires the calculation of the Z component because X,Y components are typically zero for planar or layer-wise 3D printing. This becomes a trivial solution, the sign of the difference between the Z position and plane indicates if it is above or below a plane.

For the simplest case: when the signs are opposite across a triangular edge indicates that one point lies above the plane, and the other below. This reduces the number of further operations required. The actual points for intersection across the edge, is simply the linear interpolation of the vector along the edge based on the Z distance between the mesh and the points on the edge.

This can be achieve from the following

t = plane_origin[2] - endpoints[0,:,2]
b = line_dir[:, 2]
d = np.divide(t, b)
intersection = endpoints[0]
intersection = intersection + np.reshape(d, (-1, 1)) * line_dir

Using the pre-sorted triangles from earlier, the indices can be use to extract triangle information to calculate the potential opportunities each triangle intersection may intersect. In the original Trimesh function intersection.mesh_plane, this can be done by using the local_faces and the cache_dot arguments:

new_dots = vertex_dots[myMesh.faces[intersectTris]] - z

segs = mesh_plane(
    mesh=myMesh,
    plane_origin=new_origin,
    plane_normal=plane_normal,
    local_faces=intersectTris,
    return_faces=False,
    cached_dots=new_dots)

For small meshes, this is likely to have a performance penalty for sorting and grouping of triangles. Nevertheless, for large meshes this should provide a substantial improvement because of the likelihood achieving drastic reduction in the potential number of intersection planes within the mesh. Due to independent separation of triangles within the mesh corresponding to a series of slice planes, this is trivially parallelisable as the journal articles conclude. Further optimisations were performed on the trimesh routines to increase the performance a little further.

Numba Optimisation

For further performance improvements, the intersection code in Trimesh mesh_plane was re-programmed using numba. Numba uses the LLVM compiler to JIT recompilation and convert native python bytecode re-interpretation into an intermediate representation that can be further optimised for the target architecture (CPU/GPU) that includes fine-tuning for optional multi-core processing. Additionally, some internal functions are optimised and Furthermore, type introspection is performed automatically to maxis performance.

Numba is easy to use for transforming existing code with minimal effort for prototyping purposes and has direct equivalents for numpy functions. Existing code can be translated by including the appropriate numba imports (jit, njit, prange) and prepending function (kernels) with the appropriate numba function decorator @njit. Multi-processing can be triggered with an additional argument in the @njit function decorator and fine-grained control with the use of prange. An example used on a binary search function, provided below:

from numba import njit, prange

@njit(parallel=False)
def searchsorted_parallel(a, b):
    res = np.empty(len(b), np.intp)
    for i in prange(len(b)):
        res[i] = np.searchsorted(a, b[i])
    return res

From experience, parallel performance in the particular use case offered little advantage. Another observation is that the use of single precision (float32) compared to the default double precision (float64) type offers negligible performance gains at least on an ARM64 architecture.

Note the implementation described in the post ignores sequential sorting of the generated line segments into correctly orientated connected paths that can construct polygons.

Benchmarking

Conclusions

The approach implemented is relatively simple to introduce and can be built already on existing functionality provided by trimesh. The post does not intend to infer that Trimesh’s slicing functionality is insufficient, rather, by building upon existing functionality it can enhance its performance substantially. This is even without requiring programming optimisation offered by numba. Further work is potentially needed to optimise the creation of ordered connecting paths, this may investigated further.

The functionality will be incorporated into a dedicated slicer class in a future upcoming release of PySLM.