Tag: Slicing

Super Slicing Performance in 3D Printing with PySLM

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.


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]:

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.

3D Printing Slicer - Slicing Process for a Triangular Mesh in Additive Manufacturing

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.

3D Printing Slicer - Slicing Process for a Triangular Mesh in Additive Manufacturing. Process interpolates the Z-Position of the edge intercepts to form a line segment.

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(

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

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.


The optimised slicing approach was applied on a double gyroid lattice structure that was subdivided consisting of 6.3 million triangles, with triangle intersections evenly distributed across the volume (k=10 800 ± 1090). The slicing was performed across 1000 individual slices, with the reported slicing time (extracting un-ordered segments) that was measured with the various approaches tried, that are currently used in PySLM:

The above approaches were based on a single processor/core using Python 3.8, and the performance was evaluated across 7 runs using %timeit in the ipython console and the results are presented on the side. Unlike the default implementation in PySLM, these do not generated ordered polygons via the shapely library.

It can be observed that there is a improvement in the slicing performance by ordering and grouping of intersected triangles.

Benchmarking the slicing performance using three approaches on a gyroid matrix mesh (6.3 million triangles)


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.


1 Minetto, R., Volpato, N., Stolfi, J., Gregori, R. M. M. H., & da Silva, M. V. G. (2017). An optimal algorithm for 3D triangle mesh slicing. Computer-Aided Design92, 1–10. https://doi.org/10.1016/j.cad.2017.07.001
2 Bhandari, S. (2022). A graph-based algorithm for slicing unstructured mesh files. Additive Manufacturing Letters3, 100056. https://doi.org/10.1016/j.addlet.2022.100056

GPU 3D Printing Slicer for DLP/Jetting using PySLM

Anycubic Digital Light Projection (DLP) 3D Printer System used for Slicing
Anycubic DLP 3D Printer System

Digital Light Projector (DLP) 3D Printers are an exceptionally productive technique for producing highly detailed (30<𝝁m) parts at high speeds at minimal costs.The CLIP process is a further enhancement in build speeds.

Briefly, the DLP process is similar to Stereolithography (SLA). It cures a vat of UV curable polymer material above a flexible transparent PTFE membrane. Instead of a single exposure (UV laser) into the resin, a monochrome LCD screen is used to mask the UV exposure source underneath. A greyscale bitmap image is used for each layer. Typically for most systems, after exposing the layer (1.5-3 s), the upper build-platform retracts, and mechanically pulls the cured layer away from the flexible membrane and the process is repeated. Surprisingly simple, but effective in cost and the production speed.

Additionally, bitmap based approaches are used amongst Material Jetting (MJ) technologies predominantly used within our research group CfAM, at Nottingham. Both DLP and Material Jetting offer high resolution between 30-100 𝝁m both in the XY slice plane dependent on the printer, and for Inkjet downwards of 1-10 𝝁m layer thickness depending on the choice of ink loading. Accordingly, these high resolutions are demanding to print. I came accustomed to using these printers in our CfAM lab at Nottingham on a recent project. The affordability of these printers is genuinely remarkable, owing to their mechanical simplicity.

Based on a previous post back in 2016 by Dr Matt Keeter, this is an excellent reference to an approach using WebGL implementation. Their post introduced the method, but the approach was obscured by its WebGL based implementation. Frutstratingly, I never came across an implementation for use in a research environment. These appraoches are most likely used in the free slicer software provided for desktop DLP 3D Printers.

DLP 3D Printer - Anycubic Mono 4k - Nottingham Lab
Anycubic Mono 4k DLP Printer at the University of Nottingham’s CfAM Lab loaded with a composite ink. 2022.

Interestingly, this approach can also be extended for generating 3D voxel models, by applying the project across multiple directions. However, the reliability of such method for non-manifold meshes would likely be limiting.

Method for Bitmap Slicing

The method is similar and use the same infrastructure to that used in the previous post for performing height map ray projection. Likewise, to provide a cross-platform compatibility, the use of Vispy and OpenGL 2.0 GLSL shaders are utilised within a single script. As such, the resolution of the output is limited to the maximum framebuffer size supported by the GPU driver on the system.

The approach for generating slices relies on having a connected watertight with surface triangles normals correctly orientated (fixable natively using Trimesh). The approach uses a combination of Stencil buffers integrated natively in GPU hardware.

By choosing an appropriate Z-clipping plane for the camera, the Stencil buffer is used to keep and discard rasterised triangles with the z-clipping range based on the Z-order. In order to determine if the fragments rendered are inside or outside the mesh. The render pipeline uses three passes:

  • Pass 1: stencil buffer increments on front facing fragments
  • Pass 2: stencil buffer decrements on back facing fragments
  • Pass 3: discard fragments where the stencil buffer is zero

During all the render passes, GL Depth tests are turned off. Typically in 3D Programs, triangles that are obscured from view of the 3D camera, or hidden behind other triangles are culled and the fragment is discarded prior to rendering . In this method, depth testing is turned off. The full approach is detailed further in the excerpt below inside the on_draw call.

def on_draw(self, event):

    with self._fbo:
        # Set the GL state
        gloo.set_state(blend=False, stencil_test=True, depth_test=False, polygon_offset_fill=False, cull_face=False)

        # Set the size of the framebuffer to fit the geometry with the correct aspect ratio
        gloo.set_viewport(0, 0, self._visSize[0], self._visSize[1])
        gloo.set_clear_color((0.0, 0.0, 0.0, 0.0))
        # Clear the framebuffer

        self.program['bounds'] = self.bbox[0,2], self.bbox[1,2]
        self.program['aspect'] = self.physical_size[1] /  self.physical_size[0]
        # The position of the slice position passed to the GLSL shader
        self.program['frac'] = self._z * 2.0 

        # Draw twice, adding and subtracting values in the stencil buffer

        # Render Pass 1 (Increment Stencil Buffers)
        gloo.set_stencil_func('always', 0, 0xff)
        gloo.set_stencil_op('keep', 'keep', 'incr', 'back')
        gloo.set_stencil_op('keep', 'keep', 'keep', 'front')
        self.program.draw('triangles', self.filled_buf)

        # Render Pass 2 (Decrease Stencil Buffers)
        gloo.set_stencil_op('keep', 'keep', 'decr', 'front')
        gloo.set_stencil_op('keep', 'keep', 'keep', 'back')
        self.program.draw('triangles', self.filled_buf)

        # Clear only the color buffer
        gloo.clear(color=True, depth=False, stencil=False)

        # Render Pass 3
        gloo.set_stencil_func('notequal', 0, 0xff)
        gloo.set_stencil_op('keep', 'keep', 'keep')
        self.program.draw('triangles', self.filled_buf)

        # Store the final framebuffer 
        self.rgb = _screenshot((0, 0, self._visSize[0], self._visSize[1]))

The GLSL shaders are not particularly interesting. Focus should be given to the Vertex shader, rather than the Fragment shader. This Vertex shader processes vertices of the mesh and applies the Model View Projection (MVP) transformation matrix onto the input mesh. The MVP matrix is chosen to scale the entire geometry so that it fits within the Z-clipping range of Z = -1 to +1, and is within the scope of rendering into Stencil buffer whilst using the 3D Orthographic Camera. Finally, the model is transformed based on a fractional range (0-1) to obtain the required Z-slicing plane. An epsilon value is provided for round-off purposes.

uniform   mat4 u_model; // Model transform matrix
uniform   mat4 u_view;
uniform   mat4 u_projection;

uniform  vec2 bounds; // Z bounds
uniform  float frac;  // Z fraction (0 to 1)
uniform  float aspect; // Aspect ratio

attribute vec3 a_position;

#define EPSILON 0.001

void main() {

    vec3 pos = a_position;
    // Ensure the bottom of the part is positioned to z=0 using the bottom bounding box
    pos.z -= bounds[0];
    // Scale the so that it fits within the clipping range (-1.0 < z < 1.0)
    pos.z *= -2.0/(bounds[1]-bounds[0]);
    // Adjust the position of  the verticies 
    pos.z -= frac;  
    gl_Position = u_projection * u_view * u_model * vec4(pos, 1.0);
    gl_Position.z += 1.0 - EPSILON;

The remainder of the script sets up the infrastructure for Vispy. This is performed within the initialisation call for the script. This methods sets up the correct OpenGL state, viewport size including the use of an off-screen render and specific selection of a separate Stencil framebuffer used to render onto. Both the vertex and fragment shaders are compiled and the transformation matrix is generated based on an Orthographic projection sized to the bounding box of the geometry.

    # Window Size
    shape = int(self._visSize[1]), int(self._visSize[0])

    # Create the render texture used by default in the pipeline
    self._rendertex = gloo.Texture2D((shape + (4,)), format='rgba', internalformat='rgba32f')
    # These are not used but are for reference
    #self._colorBuffer = gloo.RenderBuffer(self.shape, format='color')
    #self._depthRenderBuffer = gloo.RenderBuffer(shape, format='depth')

    # Create the stencil buffer (8 bit component)
    self._stencilRenderBuffer = gloo.RenderBuffer(shape, format='stencil')
    self._stencilRenderBuffer.resize(shape, format=gloo.gl.GL_STENCIL_INDEX8)

    # Create FBO, attach the color buffer and depth buffer
    self._fbo = gloo.FrameBuffer(self._rendertex)

    # Set the size of the view port based on the size of the window (the bounding box)
    gloo.set_viewport(0, 0, self.physical_size[0], self.physical_size[1])
    gloo.set_viewport(0, 0, self._visSize[0], self._visSize[1])

    # Create the initial orthographic view projection transformation based on the bounding box of the geometry
    self.projection = ortho(self.bbox[1, 0], self.bbox[0, 0], self.bbox[1, 1], self.bbox[0, 1], 2, 40)
    # Identity matrix
    self.model = np.eye(4, dtype=np.float32)

     # Set MVP variables for shaders
    self.program['u_projection'] = self.projection
    self.program['u_model'] = self.model
    self.program['u_view'] = self.view

Other operations are processing and the Trimesh and correctly transformed into the correct position:

The script was applied to a porous aerofoil structure with an XY resolution of 20 µm that was used previously on an Anycubic DLP system. Below is an example cross-section taken using this approach. Notice the high resolution

GPU 3D Printing Slicer used on an aerofoil structure


The overall approach may have a limited use by itself. Generally, the need to bespoke high resolution slices are limited at this stage. For reference, the full excerpt of the script is temporarily located here. In future, I will consider including this as another option within PySLM.

The source code can be found below or on GitHub:

PySLM Scan Path Iterator

Historical Background

An upcoming key feature in PySLM is the Iterators primarily useful for simulation studies, such as predicting thermo-mechanical behavior of scan strategies. Much of this builds upon ideas in former work that was done during my PhD for investigating the generation of residual stress in selective laser melting. In that study, MSC Marc, a commercial Finite Element analysis package was used to predict residual stresses generated during the process. The discretised position and laser parameters of the exposure from the laser was controlled by combination of Fortran User Subroutines and libSLM, the former c++ library.

Prior to running, a configuration file was passed to the simulation, specifying the SLM machine build file used for simulating the scan paths. libSLM parsed the compatible build-file and then based on the current time and increment would interpolate the position of the exposure point and laser parameters during firing. Beyond this main functionality, there were additional house keeping required for running the simulation, including passing information between the programs, and also additional tools to efficiently seek at an arbitrarily point in time, the state and position of the laser. This was necessary for restarting simulations on a HPC and for adaptive time-stepping required to keep numerical stability. For efficient seeking across the layers and each layer geometry structure was cached within a tree, that could be parsed on demand if necessary. Much of the infrastructure was excessive, although the implementation had to be written in c++ or Fortran to be used by integrated with the commercial solver.

Although it is difficult to perceive the full benefit of having a Python version of the same functionality, there are some instances and some analysis codes where this could be of benefit for modelling this and other processes as well.

Implementation of PySLM Iterators

The implementation builds upon the existing design from the original libSLM library. For all the Iterator classes, similar to most of the other pyslm.analysis module’s tools, the list of Layers and Models with the Laser Parameters should be passed:

# Iterates across layer geometries
layerGeomIter = Iterator(models, layerList) 

# Iterates across individual scan vectors - currently only ContourGeometry/HatchGeometry
scanVectorIter = ScanVectorIterator(models, layerList) 

# Generates an scan exposure point iterator
scanIter = ScanIterator(models, layerList) 

The first stage is building a time cache tree across each LayerGeometry. In practice, the cache tree structure is not necessary if the scan iterator iteratively increments along in time. Having a cache structure enables non-linear movement of the iterator across the entire build . It also provides a fast random-access lookup to seek to a specific Layer or LayerGeometry for use in simulations or analyses.

This structure is formed by iteratively measuring the total time taken to scan an individual LayerGeometry group which is stored in a tree node (TimeNode). The cumulative time taken to scan across each LayerGeometry TimeNode to provide the total scan time across the Layer. The TimeNode can be assigned child and parent nodes using the attributes (TimeNode.parent and TimeNode.children) in order to navigate across the entire tree. Each TimeNode provides a key-value pair (id, value) to store the reference LayerGeometry or Layer for simplified access.

The Cache Tree is generated and stored in the Base Class, Iterator and is generated in the private method (Iterator._generateCache) and stored in the attribute Iterator.tree.

def _generateCache(self) -> None:
 self._tree = TimeNode() for layerId, layer in enumerate(self.layers):     # Create the layer     layerNode = TimeNode(self._tree, id=layerId, value=layer)     self._tree.children.append(layerNode)     for layerGeomId, layerGeom in enumerate(layer.geometry):         geomNode = TimeNode(layerNode, id=layerGeomId, value=layerGeom)         geomNode.time = getLayerGeometryTime(layerGeom, self._models)         layerNode.children.append(geomNode) self._cacheValid = True

The Iterator class has many useful facilities, such as build-time estimation, seeking access to the Layer or LayerGeometry at an arbitrary point in time. The class stores additional info such as the layer dwellTime – this can be re-implemented in a derived class. For implementing the iterator behavior used across all dependent classes it also stores the current time and reference pointers to the current Layer and LayerGeometry. Essentially the Iterator class can be used to iterate across each LayerGeometry within a build as a foundation to the other class. Each of these Iterator classes builds upon the magic methods available in Python: __iter__ and __next__ . The __iter__ method simply sets up the object and re-initialises the Iterators attributes. Once the cache tree is generated internally, it offers no penalty to generate a new iterator . Below is an excerpt taken from the ScanVectorIterator:

def __iter__(self):
    self._time = 0.0
    self._layerGeomTime = 0.0
    self._layerInc = 0
    self._layerGeomInc = 0

    return self

def __next__(self):
     if self._layerScanVecIt < len( self._layerScanVectors):
         scanVector = self._layerScanVectors[self._layerScanVecIt]
         self._layerScanVecIt += 1
         return scanVector
         # New layer
         if self._layerIt < len(self._layers):
             layerVectors = self.getLayerVectors(self._layers[self._layerIt])
             self._layerScanVectors = self.reshapeVectors(layerVectors)
             self._layerScanVecIt = 0
             self._layerIt += 1
             return self.next()
             raise StopIteration

The Iterator class and ScanVectorIterator class do not require much further attention, as the pointer to the geometry is incremented only. The ScanIterator class, however, is more useful for simulation and will be discussed further.

Scan Iterator Class

The ScanIterator class is used for incrementally advancing the exposure source across each scan vector. This is particularly important for visualising or simulating the AM process. The time increment is based on a chosen but adjustable timestep, and the laser parameters across each scan vector (i.e. the effective scan velocity) obtained from the assigned BuildStyle.

The exposure point is linearly interpolated across each scan vector based on the current time within the LayerGeometry depending on the type. For identifying the position, the cumulative distance is captured and the current timeOffset for the layer geometry is used to estimate the distance covered by the exposure source across the entire LayerGeometry section. For simplicity this assumes no acceleration terms and uses a constant velocity profile. Based on the timeOffset, the scan vector is obtained and then the final position is interpolated across the scan vector.

laserVelocity = getEffectiveLaserSpeed(buildStyle)

# Find the cumulative distance across the scan vectors in the LayerGeometry (Contour)
delta = np.diff(layerGeom.coords, axis=0)
dist = np.hypot(delta[:,0], delta[:,1])
cumDist = np.cumsum(dist)
cumDist2 = np.insert(cumDist, 0,0)

# If the offsetDist calculated is outside of the cumulative distance then some error has occured
if offsetDist > cumDist2[-1]:
    raise Exception('Error offset distance > cumDist {:.3f}, {:.3f}'.format(offsetDist, cumDist2[-1]))

id = 0

# Find the id of the current scan vector given the time offset
for i, vec in enumerate(cumDist2):
    if offsetDist < vec:
        id = i

# interpolate the position based on offset distance in the scan vector
linearOffset = (offsetDist - cumDist2[id-1]) / dist[id-1]
point = layerGeom.coords[id-1] + delta[id-1] * linearOffset

The above example is specifically for the contour geometry. Note the for loop is not particularly efficient but serves its purpose for identifying the Iterator’s current scan vector.

Iterator Use:

Each iterator can be subsequently called after using the iter method in a variety of pythonic ways:

#Create a scan vector iterator
ScanVectorIterator(models, layerList)

# Create a python iter object from a ScanVectorIterator
scanIter = iter(scanVectorIter)

# Get a single scan vector
firstScanVec = more(scanIter)

# Collect all the remaining scan vectors
scanVectors = np.array([point for point in scanIter])

Current Limitations:

Note the current implementation of the iterators currently only consider ContourGeometry and HatchGeometry and does not include PointGeometry groups. The jump vectors are ignored, which will have a small but in most situations a negligible effect on the the overall accuracy of the timing used for the iterators.

Another obvious limitation is that this only accounts for single exposure source systems. It is not known to myself, how multiple-exposure systems scan (i.e. are they truly in parallel based on the laser number) or is there is some built-in machine heuristic which balances the scanning across all laser sources and spatially – e.g. to prevent overheating. This depends on the SLM system such as if multiple exposure sources are limited by zones or have full areal access to the bed. Anyone’s comments or experiences on this aspect would be sincerely welcomed.


An example showing the basic usage and functions available with the Iterator classes are available in the Github Repo examples/example_laser_iterator.py

Multi-threading Slicing & Hatching in PySLM

In PySLM, the slicing and hatching problem is inherently parallelisable simply because we can discretise areas the geometry into disrcete layers that for most situations behave independent. However, the actual underlying algorithms for slicing, offsetting the boundaries, clipping the hatch vectors is serial (single threaded). In order to significantly reduce the process time, multi-threaded processing is very desirable

Multi-threading in Python

Unfortunately, Python like most scripting or interpreter languages of the past are not inherently designed or destined to be multi-threaded. Perhaps, this may change in the future, but other scripting languages may fill this computational void (Rust, Julia, Go). Python, by intentions limits any multi-threaded use in scripts by using a construct known as the GIL – Global Interpreter Lock. This is a shared situation in other common scripting languages Matlab (ParPool), Javascript (Worker) where the parallel computing capability of multi-core CPUs cannot be exploited in a straightforward manner.

To some extent special distributions such as via the Anaconda distribution, core processing libraries such as numpy, scipy, scikit libraries internally are generally multi-threaded and support vectorisation via native CPU extensions. More computational mathematical operations and algorithms can to some extent be optimised to run in parallel automatically using numba, numexpr, and however, this cannot cover more broad multi-functional algorithms, such as those used in PySLM.

Python has the thread module for multi-threaded processing, however, for CPU bound processing it has very limited use This is because Python uses the global interpreter lock – GIL and this only allows one programming thread (i.e. one line) to be executed at any instance. It is mainly used for asynchronous IO (network or filesystem) operations which can be processed in the background.

Use of Multiprocessing Library in PySLM

The other option is to use the multiprocessing library built into the core Python distribution. Without going into too much of the formalities, multi-processing spawns multiple python processes and assign batches of work. The following programming article I found as a useful reference to the pitfalls of using the library.

In this implementation, the Pool and Manager modules are used to more optimally process the geometry. The most important section is to initialise the multiprocessing library with the ‘spawn‘ method, which stops random interruptions during the operation as discussed in the previous article.

from multiprocessing import Manager
from multiprocessing.pool import Pool
from multiprocessing import set_start_method


The Manager.dict acts as a ‘proxy‘ object used to more efficiently store data which is shared between each process that is launched. Without using manager, for each process launch, a copy of the objects passed are made. This is important for the geometry or Part object, which if it were to contain a lattice of a set ofs complex surface would become expensive to copy.

d = Manager().dict()
d['part'] = solidPart
d['layerThickness'] = layerThickness # [mm]

A Pool object is used to create a set number of processes via setting the parameter processes=8 (typically one per CPU core). This is a fixed number re-used across a batch through the entire computation which removes the cost of additional copying and initialising many process instances. A series of z slice levels are created representing the layer z-id. These are then merged into a list of tuple pairs with the Manager dictionary and is stored in processList.

Pool.map is used to perform the slice function (calculateLayer) and collect all computed layers following the computation.

p = Pool(processes=8)

numLayers = solidPart.boundingBox[5] / layerThickness
z = np.arange(0, numLayers).tolist()

processList = list(zip([d] * len(z), z))

# Run the pro
layers = p.map(calculateLayer, processList)

The slicing function is fairly straightforward and just unpacks the arguments and performs the slicing and hatching operation. Note: each layer needs to currently initialise its own instance of a Hatcher class because this is not shared across all the processes. This carries a small cost, but means each layer can process entirely independently; in this example the change is the hatchAngle across layers. The layer position is calculated using the layer position (zid) and layerThickness.

def calculateLayer(input):
    # Typically the hatch angle is globally rotated per layer by usually 66.7 degrees per layer
    d = input[0]
    zid= input[1]

    layerThickness = d['layerThickness']
    solidPart = d['part']

    # Create a StripeHatcher object for performing any hatching operations
    myHatcher = hatching.Hatcher()

    # Set the base hatching parameters which are generated within Hatcher
    layerAngleOffset = 66.7
    myHatcher.hatchAngle = 10 + zid * 66.7
    myHatcher.volumeOffsetHatch = 0.08
    myHatcher.spotCompensation = 0.06
    myHatcher.numInnerContours = 2
    myHatcher.numOuterContours = 1
    myHatcher.hatchSortMethod = hatching.AlternateSort()

    #myHatcher.hatchAngle += 10

    # Slice the boundary
    geomSlice = solidPart.getVectorSlice(zid*layerThickness)

    # Hatch the boundary using myHatcher
    layer = myHatcher.hatch(geomSlice)

    # The layer height is set in integer increment of microns to ensure no rounding error during manufacturing
    layer.z = int(zid*layerThickness * 1000)
    layer.layerId = int(zid)

    return zid

The final step to use multiprocessing in Python is the inclusion of the python __name__ guard i.e:

if __name__ == '__main__':

The above is unfortunate because it makes debugging slightly more tedious in editors, but is the price for extra performance.

Performance Improvement

The performance improvement using the multiprocesssing library is shown in the table below for a modest 4 core laptop (my budget doesn’t stretch that far).

PySLM: A matplotlib showing hatching and slicing across multiple layers for the a cubic geometry using the python multi-processing library.
Matplotlib Figure showing every subsequent 10 layers of hatching for the geometry but is shown reduced scale.

This was performed on the examples/inversePyramid.stl geometry with an overall bounding box size [90 \times 90 \times 60] mm, hatch distance h_d=0.08mm and the layer thickness set at 40 μm.

Number of ProcessesRun Time [s]
Base-line (Simple For loop)121
Approximate timings for 4 core CPU i7 Processors Using the Multi-Processing Library.

Given these are approximate timings, it is nearly linear performance improvement for the simple example. However, it can be seen choosing more processes beyond cores, does squeeze some extra performance out – perhaps due to Intel’s hyperthreading. Below shows that the CPU is fully utilised.

PySLM: Multi-threading options


This post shows how one can adapt existing routines to generate multi-processing slicing and hatching with PySLM. In the future, it is desirable to explore a more integrated class structure for hooking functions onto. Other areas that are of interest to explore are potentially the use of GPU computing to parallelise some of the fundamental algorithms.


The example showing this that can be run is example_3d_multithread.py in the github repo.

Slicing and Hatching for Selective Laser Melting (L-PBF)

Much of slicing and hatching process is already taken for granted in commercial software mostly offered by the OEMs of these systems rarely discussed amongst academic research. Already we observe practically the implications direct control over laser parameters and scan strategy on the quality of the bulk material – reduction in defects, minimising distortion due to residual stress, and the surface quality of parts manufactured using these process. Additionally, it can have a profound impact the the metallic phase generation, micro-structural texture driven via physics-informed models [1], grading of the bulk properties and offer precise control over manufacturing intricate features such as thin-wall or lattice structures [2].

This post hopefully highlights to those unfamiliar some of the basis process encountered in the generation of machine build files used in AM systems and get a better understanding to the operation behind PySLM. I have tried my best to generalise this as much as possible, but I imagine there are subtleties I have not come across.

This post is to provide some reference into the generation of hatches or scan vectors are created for use in AM processes such as selective laser melting (SLM), which uses a point energy source to raster across a medium. Some people prefer to more generally to classify the family of processes using the technical ASTM F42 committee standards 52900 and 52911 – Powder Bed Fusion (PBF). I won’t go into the basic process of the manufacturing processes such as EBM, SLM, SLA, BJF, as there are many excellent articles already that explain these in far greater detail.

Machine Build Files

AM processes require a digital representation to manufacture an object. These tend to be computed offline – separate from the 3D Printer, using specialist or dedicated pre-processing software. I expect this will become a closed-loop system in the future, such that the manufacturing integrated directly into the machine.

For some AM process families, the control operations may be exceedingly granular – i.e. G-code. G-code formats state specific instructions or functional commands for the 3D printer to sequentially or linearly execute. These tend to fit with deposition methods such as Filament Extrusion, Direct-Ink-Writing (robo-casting) and direct energy deposition (DED) methods. Typically, these tends to be for deposition with a machine systems, which requires coordination of physical motion in-conjunction with some mechanised actuation to deposit/fuse material.

Machine Build File Formats for L-PBF

For exposure (laser, electron-beam) based AM processes, commercial systems use a compact notation solely for representing the scan path the exposure source will traverse . The formats are often binary to aid their compactness.

To summarise, within these build files, an intermediate representation consists of index-based referenceable parameters for the build. The remainder consists of a series of layers, that contain geometric entities (points, vectors) that are used to to control the exposure for the border or contour or raster or infill the interior region. For L-PBF processes, the digital files, commonly referred as “machine build file” comes in various flavours dependent on the machine manufacture:

  • Renshaw .mtt,
  • SLM Solution .slm,
  • DMG Mori Realizer .rea
  • EOS .sli
  • Aconity .cli+ or .ilt wrapper

Some file formats, such as Open Beam Path format can specify bezier curves [3]. Another recently proposed open source format created by RWTH Aachen in 2022 called OpenVector Format based on Google’s Protobuf schema. The format aims to offer a specification universally compatible across a swathe of PBF processes and supplement existing commercial formats with additional build-process meta-data (e.g. build, platform temperature, dosing) and detailed definition with further advancements in the process, such as multi-beam builds.

Build-File Formats

Higher level representations that describe the distribution of material(s) defining geometry – this could be bitmap slices or even a 3D model. Processes such as Jetting, BJF, High Speed Sintering, DLP Vat-polymerisation currently available offer this a reality. With time, polymer and metal processes will evolve to become 2D:, diode aerial melting [4] or more aerial based scanning based on holographic additive manufacturing methods, such as those proposed by Seurat AM [5] based off research at LLNL, and recently at University of Cambridge [6] . In the future, we can already observe the exciting prospect of new processes such as computed axial lithograph [7] that will provide us near instantaneous volumetric additive manufacturing.

For now, single and multi point exposure systems for the imminent future will remain with us as the currently available process. PySLM uses an intermediate representation – specifying a set of points and lines to control the exposure of energy into a layer.

The Slicing and Hatching Process in L-PBF

With nearly most conventional 3D printing process, it begins with a 3D representation of a solid volume or geometry. 2D planar slices or layers are extracted from a 3D mesh or B-Rep surface in CAD by taking cross-sections from a geometry. Each slice layer consist of a set of boundaries and holes describing the cross-section of an object. Note: non-planar deposition does exist for DED/Filament processes, such as this Curved Layer Fused Deposition Modeling [ref] and a spherical slicing technique [8].

For consolidating material, an exposure beam must raster across the surface medium (metal or polymer powder, or a photo-polymer resin) depending on the process. Currently this is a single or multiple point which moves at a velocity vwith a power P across the surface. The designated exposure or energy deposited into the medium is principally a function of these two parameters, depending on the type of laser:

  • (Quasi)-Continious Wave: The laser remains active switched on (typically modulated using a form of PWM) across the entire length of the scan vector
  • Pulsed Mode (Q-Switched): Laser is pulsed at set distances and exposure times across the scan vector

Numerous experiments often tend to result in parametric power/speed maps to the achieved part bulk density, that result in usually optimal processing windows that produce stable and consistent melt-tracks [9][10]. Recently, process maps are based on a non-dimensional parameter such as the normalised enthalpy approach, that more reliably assist selecting a suitable process windows [11].

Illustration of a scan vector commonly used in Laser Powder-Bed Fusion (SLM)

However, the complexity of the process extends further and is related to a many additional variables dependent on the process such as layer thickness, absorption coefficient (powder and material), exposure beam profile etc.. Additionally, the cumulative energy deposited spatially over a period of time must consider overlap of scan vectors within an area.

Scan Vector Generation

Each boundary polygon is offset initially to account for the the radius of the beam exposure, which is termed a ‘spot compensation factor‘. Some processes such as SLS or BJF account for global part shrinkage volumetrically throughout the part by having a global scale factor or deformed mesh to compensate to non-uniform shrinkage across the part.

The composition of laser scan vectors used in a slice or layer for L-PBF or Selective Laser Melting. The boundary is offset multiple times, with the interior or core filled with hatch vectors.
The typical composition of a layer used for scanning in exposure based processes. This consists of outer and inner contours, with the core interior filled with hatches.

This first initial offset is the outer-contour which would be visible on the exterior of the part. This contour will have a different set of laser parameters in order to optimise and improve the surface roughness of the part obtained. A further offset is applied to generate a set of inner-contours before hatching begins.

Depending on the orientation of the surface (e.g. up-skin or down-skin), the boundary and interior region may be intersected to fine-tune the laser parameters to provide better surface texture, or surface roughness – typically varying between Ra = 3-13 μm [12] primarily determined by the surface angle and a combination of the process variables including,

  • the powder feedstock (bulk material, powder size distribution)
  • laser parameters
  • layer thickness (pre-dominantly fixed or constant for most AM processes)

Overhang regions and surfaces with a low overhang angles tend to be susceptible to high surface-roughness. Roller re-coater L-PBF systems – available only on 3DSystems or AddUp system,, tend to offer far superior surface quality on low inclined or overhang regions. Additionally, progressive advancement and maturity of laser parameter optimisation, and those computationally driven using part geometry [13] are able to further enhance the quality and potentially eliminate the need for support structures. Depending on the machine platform, these regions are identified by sampling across two-three layers. Overhang regions obviously require support geometry, which is an entirely different topic discussed in this post.

Laser parameters in SLM (L-PBF) can be optimised based on the adjacent surface regions. Special regions, include the upskin, downskin and overhang regions
Laser parameters can be optimised based on the adjacent surface regions. Special regions, include the upskin, downskin and overhang regions needed to improve the surface roughness and reduce density in regions.

Following the generation of the contours, the inner core region requires filling with hatches. Hatches are a series of parallel scan vectors placed adjacent at a set hatch distance, h_d. This parameter is optimized according to the material processed, but is essentially related to the spot radius of the exposure point r_s in order to reduce inter-track and inter layer porosity. Across each layer these tend to be placed at a particular orientation \theta_h, which is is then incrementally rotated globally for subsequent layers, typically 66.6°. This rotation aims to smooth out the build process in order to minimise inter-track porosity, and generate homogeneous material, and in the case of SLM mitigate the effects of anisotropic residual stress generation.

The composition and terminology (hatch distance, hatch spacing, hatch angle) used in L-PBF. The Layer Geometry objects used to scan across a Layer in Selective Laser Melting (L-PBF). The various parameters such as the hatch distance and hatch angle are shown.
A general composition of the various LayerGeometry objects used to scan across a Layer. The various parameters such as the hatch distance, spacing and hatch angle are shown.

The distribution (position, length, rotation) of these hatch vectors are arranged using a laser scan strategy. The most common include a simple alternating hatch, stripe and island or checkerboard scan strategy.

Each set or group of scan vectors is stored together in a LayerGeometry, depending on the type (either a set of point exposures, contour or hatch vectors). These LayerGeometry groups usually share a set of exposure parameters – power, laser scan speed (point exposure time, point distance for a pulsed laser), focus position).

Some systems offer a greater degree of control and can control individual power across the scan vectors. Other can fine tune the acceleration and modulate the power along the scan vectors to support techniques known as ‘skywriting‘. For instance in SLM, it has been proposed that careful tuning of the laser parameters towards the end of the scan vector, i.e. turning can reduce porosity by preventing premature collapse of key holing phenomena [14]. In theory, PySLM could be extended to provide greater control of the electro-optic systems used in the process if so desired.

Hopefully, this provides enough background for those who are interested and engaged in working with developing scan strategies and material development using PySLM.


1 Plotkowski, A., Ferguson, J., Stump, B., Halsey, W., Paquit, V., Joslin, C., Babu, S. S., Marquez Rossy, A., Kirka, M. M., & Dehoff, R. R. (2021). A stochastic scan strategy for grain structure control in complex geometries using electron beam powder bed fusion. Additive Manufacturing46. https://doi.org/10.1016/j.addma.2021.102092
2 Ghouse, S., Babu, S., van Arkel, R. J., Nai, K., Hooper, P. A., & Jeffers, J. R. T. (2017). The influence of laser parameters and scanning strategies on the mechanical properties of a stochastic porous material. Materials and Design131, 498–508. https://doi.org/10.1016/j.matdes.2017.06.041
3 Open Beam Path – Freemelt, https://gitlab.com/freemelt/openmelt/obplib-python
4 Zavala Arredondo, Miguel Angel (2017) Diode Area Melting Use of High Power Diode Lasers in Additive Manufacturing of Metallic Components. PhD thesis, University of Sheffield.
5 Seurat AM. https://www.seuratech.com/
6 https://www.theengineer.co.uk/holographic-additive-manufacturing-lasers/
7 Kelly, B., Bhattacharya, I., Shusteff, M., Panas, R. M., Taylor, H. K., & Spadaccini, C. M. (2017). Computed Axial Lithography (CAL): Toward Single Step 3D Printing of Arbitrary Geometries. Retrieved from http://arxiv.org/abs/1705.05893
8 Yigit, I. E., & Lazoglu, I. (2020). Spherical slicing method and its application on robotic additive manufacturing. Progress in Additive Manufacturing, 5(4), 387–394. https://doi.org/10.1007/s40964-020-00135-5
9 Yadroitsev, I., & Smurov, I. (2010). Selective laser melting technology: From the single laser melted track stability to 3D parts of complex shape. Physics Procedia, 5(Part 2), 551–560. https://doi.org/10.1016/j.phpro.2010.08.083
10 Maamoun, A. H., Xue, Y. F., Elbestawi, M. A., & Veldhuis, S. C. (2018). Effect of selective laser melting process parameters on the quality of al alloy parts: Powder characterization, density, surface roughness, and dimensional accuracy. Materials, 11(12). https://doi.org/10.3390/ma11122343
11 Ferro, P., Meneghello, R., Savio, G., & Berto, F. (2020). A modified volumetric energy density–based approach for porosity assessment in additive manufacturing process design. International Journal of Advanced Manufacturing Technology, 110(7–8), 1911–1921. https://doi.org/10.1007/s00170-020-05949-9
12 Ni, C., Shi, Y., & Liu, J. (2019). Effects of inclination angle on surface roughness and corrosion properties of selective laser melted 316L stainless steel. Materials Research Express, 6(3). https://doi.org/10.1088/2053-1591/aaf2d3
13 Velo3D Sapphire Printer – SupportFree Technology. https://blog.velo3d.com/blog/supportfree-what-does-it-mean-why-is-it-important
14 Martin, A. A., Calta, N. P., Khairallah, S. A., Wang, J., Depond, P. J., Fong, A. Y., … Matthews, M. J. (2019). Dynamics of pore formation during laser powder bed fusion additive manufacturing. Nature Communications, 10(1), 1–10. https://doi.org/10.1038/s41467-019-10009-2