Wolfram Cloud Document

The Wolfram Emerging Leaders Program: High School (WELP:HS) is a remote, semester-long program designed for exceptional graduates from the Wolfram High School Summer Research Program to continue to develop their skills. Students work in a small group over the course of a semester to do either a research-based project or a product development–based project under the guidance of expert mentors. Projects are centered around a topic at the intersection of the group’s interests.

Authors:

◼

Austin Jiang

◼

Kevin Toren

◼

Ritvik Gupta

Abstract

Traditional infill structures commonly used in 3D printing applications are often based on human assumptions about the ideal shapes for strength, flexibility, and efficient use of material. They take the form of predefined geometric patterns (gyroid, cubic, etc), which are inherently limited in behavior by their complicated rules. We propose the use of Cellular Automaton based algorithms for the generation of optimized 3D printer infill. Unlike traditional algorithms, Cellular Automata (CA) displays emerging complexity from simple rules, allowing for natural patterns unattainable with human-induced constraints. On a 1-dimensional scale, it is known that just 255 elementary rules display a wide range of behaviors (including Turing completeness). In 3 dimensions, there are trillions of possible CA rules, signifying a high probability that there exist rules that provide an infill to conform to specific constraints such as lightness, flexibility, and strength. In the WELP: HS program, we created and tested our CA infill on 3D prints using a discretized voxel representation and used the Marching Cubes algorithm to export to an STL that can be printed with ubiquitous slicing programs.

Generating 3D CA

Just as a single evolution of a 1D cellular automata can be represented as a line, and a 2D cellular automata by a plane, we can use a volume to represent the current state of a 3D cellular automata. For each “voxel” (3D cell) in the state, we will consider neighbors in a 1-cell radius. Based on the provided rule, each voxel will update, producing another 3D state.

Representing a Rule

The type of CA we will consider in this project are called totalistic. In a totalistic CA, rules depend only on the count of active cells in a neighborhood, not their specific spatial arrangement. For example, consider binary neighborhood patterns represented as strings where ‘0’ indicates an inactive (white) cell and ‘1’ indicates an active (black) cell. The patterns “0100”, “1000”, “0010”, and “0001” would all be treated identically by a totalistic rule since they each contain exactly one active cell. By treating the majority of 3D states as symmetries of each other, we dramatically reduce the amount of possible rules we have to search through while still maintaining the complexity contained in the more repetitive set of non-totalistic CA. For our 3D implementation, we use two states (black and white for filled and empty space, respectively) and consider all cells within a radius of 1 as neighbors (27 cells total in the 3D neighborhood).

Evolving States

Calling the CellularAutomaton function with a rule and an initial state will apply the rule once. A common initial state is one black cell in the center of the matrix. One way to see how the CA will “evolve” is to recursively call CellularAutomaton and plot them in 3D space (after all, they will end up becoming real printed things).

In[]:=

rule={14,{2,1},{1,1,1}};initial=CenterArray[1,{5,5,5}];ArrayPlot3D/@NestList[CellularAutomaton[rule,#]&,initial,4]

Out[]=





Wolfram allows us to “look ahead” to whatever evolution we want. But we will continue with this step-by-step method because it allows us to iteratively test each evolution of a certain rule.

Creating Boundaries

Ultimately, we want our infill patterns to go inside of a pre-existing 3D print. So, we have to find a way to limit the “growth” of a CA within certain bounds. The CellularAutomaton function already limits CA into a cubic space using a built-in “boundary condition” which links opposite faces together (like folding a cube into a torus, or running off the right-side of the screen to reappear on the left-side in a videogame).

For any shape more complicated than a cube, a custom boundary condition test needs to be implemented. Our solution is to apply a “bitmask” operation that zeros out cells lying outside the target shape. While this masking operation may disrupt the perfect symmetry that CAs often display in a cubic domain, it preserves the structurally supportive patterns we’re interested in for infill applications. So, while the periodic boundary conditions handle containment of the CA computation, the bitmask operation prepares the pattern for insertion into an arbitrary 3D model. However, before implementing this masking, we need a way to represent arbitrary 3D shapes in the same discrete matrix format as our CA patterns.

Sphere

Creates a binary 3D matrix where cells within radius r of the origin are 1, effectively “chopping” a CA pattern into a spherical shape after applying this bitmask. The fundamental shape of the sphere always stays the same, but a larger radius makes it more detailed.

In[]:=

sphereBitmask[r_]:= Table If[Norm[{x, y, z}] <= r, 1, 0], {x, -r, r}, {y, -r, r}, {z, -r, r};

In[]:=

ArrayPlot3D[sphereBitmask[5]]

Out[]=

Cylinder

Creates a binary 3D matrix where cells within radius r of the z-axis and between 0 and h in height are 1, effectively “chopping” a CA pattern into a cylindrical shape after applying this bitmask. The shape of the cylinder can be changed by altering the ratio between height and radius. The detail can be increased by increasing the radius and height in proportion.

In[]:=

cylinderBitmask[r_, h_]:= Table If[Norm[{x, y}] <= r, 1, 0], {x, -r, r}, {y, -r, r}, {z, 0, h};

In[]:=

ArrayPlot3D[cylinderBitmask[5,30]]

Out[]=

Rectangular Prism

A rectangular prism does not “chop off” any voxels since it fully fills the available volume. This will be the default whenever a bitmask shape is required (equivalent to letting a CA evolve without specific boundary constraints).

In[]:=

prismBitmask[bounds_]:= ConstantArray[1, bounds];

In[]:=

ArrayPlot3D[prismBitmask[{10,10,10}]]

Out[]=

Voxelizer

We started with simple polyhedra, but for more complex 3D models, we need to work with meshes. To convert a mesh into a matrix, we use the voxelizer function, which discretizes a continuous shape into 3D voxels. Increasing the resolution makes the voxels smaller relative to the model, allowing for finer detail. Since we are operating on meshes which can represent any 3D shape within Wolfram, this function is compatible with any 3D model we decide to use.

The voxelizer function takes a 3D model and converts it into a grid of small cubes (voxels), representing the shape with a given resolution. Here’s how it works:

First, we determine the boundaries of the model in each dimension.

bounds=RegionBounds[model];

This gives the minimum and maximum x, y, and z coordinates that enclose the model.

We calculate the voxel size based on the dimension of the bounding box, dividing it by the amount of voxels as specified by resolution.

d=Min[(bounds[[All,2]]-bounds[[All,1]])]/resolution;

This means the model is divided into small cubes of side length d.

For each axis we create a list with all of the coordinates representing filled voxels in the infill pattern.

xRange=Range[bounds[[1,1]]+d/2,bounds[[1,2]],d];yRange=Range[bounds[[2,1]]+d/2,bounds[[2,2]],d];zRange=Range[bounds[[3,1]]+d/2,bounds[[3,2]],d];

We loop through all (x, y, z) combinations and determine if they are located inside the model.

Table[If[RegionMember[model,{x,y,z}],1,0],{x,xRange},{y,yRange},{z,zRange}]

For each position in the grid, RegionMember checks if the point belongs to the model. If it does, the voxel is marked as 1, otherwise, it’s 0. The end result is a discrete matrix of voxels which follow the general shape of the mesh.

Here is the full version of the voxelizer function that takes a 3D model and resolution as input and discretizes it into a matrix of particular size as determined by the resolution.

In[]:=

voxelizer[model_, resolution_]:= Module {bounds, d, xRange, yRange, zRange}, bounds = RegionBounds[model]; d = Min[(bounds[[All, 2]] - bounds[[All, 1]])]resolution; xRange = Range[bounds[[1, 1]] + d/2, bounds[[1, 2]], d]; yRange = Range[bounds[[2, 1]] + d/2, bounds[[2, 2]], d]; zRange = Range[bounds[[3, 1]] + d/2, bounds[[3, 2]], d]; Table If[RegionMember[model, {x, y, z}], 1, 0], {x, xRange}, {y, yRange}, {z, zRange} ;

Here is the visualization of the matrix created by voxelizer for a Parallelepiped.

Integrated Function

Now that we have declared all the methods for generating a bitmask, we need to integrate it with the built-in CellularAutomaton function. As stated previously, CellularAutomaton already acts as if it were using a rectangular prism bitmask, so this function is only relevant for more complicated boundaries (it will still be exclusively used when we need to generate a CA in this project). The “chop off” action is performed by performing bitwise AND between the bitmask matrix and the infill matrix, where only infill voxels that are also present in the bitmask matrix are kept (1 AND 1 --> 1). The rest of the voxels are “chopped off,” effectively limiting growth to within the bitmask shape and allowing for the infill pattern to be seamlessly combined with the actual 3D model in a later section. If we did not apply the bitmask, then it could be possible for an infill pattern to extend outside the bounds of the 3D model, making it unprintable.

This function serves as the bridge between our voxelized 3D models and the CA patterns we’ll use to fill them, ensuring our infill patterns remain perfectly contained within the target shape at each step of their evolution.

Creating Metrics

To evaluate the performance of 3D CA-based infills against real-world printing needs, we need to establish some metrics. Traditional 3D printing infills, such as gyroid and cubic patterns, follow predetermined shapes with known strength and flexibility. However, CA-based infills are highly customizable and can evolve some unique patterns that have never been seen. We decided to focus on infill percentage and connectivity to objectively compare CA-generated infills to each other to find the most suited structures to test physically, largely increasing the efficiency for testing.

Infill Percentage

Infill percentage is defined as the ratio of filled cells of the slice over the overall volume of the 3D model. This metric is crucial because it can determine the needed amount of material and the approximate efficiency of printing the model. Prioritizing the infill percentage directly correlates to the strength, time, and material cost for any 3D printed model.

Connectivity

Connectivity quantifies the degree to which cells in the structure are interconnected. A higher connectivity indicates a more cohesive and a potentially stronger structure. To calculate the connectivity of a model, we need to filter out all cells in the slice that are not connected to the bottom layer of the slice or suspended. We divide the total number of the non-filtered cells by the total volume of the model to get the connectivity index of a slice.

To remove cells in the slice that are either disconnected from the bottom layer or suspended, we’ll apply a Breadth-First Search (BFS) starting from the bottom as a filter. The bfsConnectedParts function returns a 3D array of 1s and 0s, marking only the disconnected or suspended cells.

The bfsConnectedParts function utilizes two key parameters.

visit is a 3D array matching the dimensions of the input slice. It serves as a flag system where 1 indicates a cell that is part of the object but hasn’t been checked for connectivity and 0 represents either empty space or a cell that has already been processed.

queue is an array storing cell positions as (x, y, z) triples. It contains cells confirmed to be connected to the current component and cells whose neighboring cells still need connectivity verification.

Initially, all cells at the bottom-most layer will be in the queue because these cells are the foundation of the whole model and are highly likely to be connected to the main body.

For each cell confirmed as connected and added to the queue, the algorithm checks the cells that are horizontally adjacent and directly above it. When an adjacent cell is part of the object and has not yet been visited, it is marked as connected and added to the end of the queue. This process recursively spreads through the structure, expanding outward from each connected cell. The search continues in a while loop until no new cells can be added to the queue, which indicates that the entire connected component has been fully explored and mapped.

In terms of adding more cells into queue, the function appends a number of valid cell positions that are adjacent to the current element cur to queue. These cells will be checked to be connected as the current cell is connected and they are adjacent to the current cell.

To have access to all of a cell’s adjacent cells, we’ll use the Map function to map each relative position of adjacent cells, dir, to its absolute position. The absolute position can be found by adding the relative position to the current position of a cell. dir is defined at the beginning of the module.

The reason we don’t check the cell located underneath {0,0,-1} as connected is because we don’t want suspension.

Not every adjacent cell can be added to the queue for checking. Two critical conditions must be satisfied: the cell must be within the slice’s range and must not have been previously checked. These criteria are implemented through a detailed If statement that carefully validates each potential adjacent cell before adding it to the queue. By enforcing these conditions, the function ensures an efficient and accurate exploration of the connected structure, preventing unnecessary computations and maintaining the search’s precision.

Once an adjacent cell satisfies our conditions, it will be checked to be connected in visit.

And its absolute position will be added to queue.

This is the full version of the bfsConnectedParts function:

Now, let’s look at the connectivity function. The only parameter taken by the function is slice, which is exactly the slice that we are trying to find the connectivity for.

Within the module, the function begins with an if statement that first checks if the slice is empty or not. If it’s empty, then its connectivity is defined as zero and the function returns zero.

Otherwise, bottomLayer will store the z-index (height) of the lowest occupied cell in the slice by iterating through all labeled positions, with z as the primary loop variable (same way it is applied to the isolatedParts function later in the essay).

This way we can ensure that the cell we are finding in the bottom layer is always the lowest.

Next, we need to initialize queue by putting the positions of all cells on the bottom layer into the array in triple while checking visit as it is marked to be connected:

Next, we initialize the queue by examining all cells on the bottom layer of the 3D structure. We iterate through each position (x, y) on this layer, checking if the cell is part of the object slice[[x, y, bottomLayer]] == 1. When a cell meets this criterion, it’s marked as visited visit[[x, y, bottomLayer]] = 0 and its coordinates {x, y, bottomLayer} are added to the queue. The Flatten function with level 1 ensures that the resulting queue is a list of triples, each representing a cell’s position. This process effectively seeds the queue with all connected cells from the bottom layer, providing a starting point for the breadth-first search to explore the entire connected structure.

If cells only exist on the bottom layer no cells are disconnected, and therefore connectivity will be 1.

Lastly, we utilize the bfsConnectedParts function to determine the connectivity of the structure. This function returns the number of disconnected cells. We then calculate the ratio of disconnected cells to the total number of cells in the slice. By subtracting this fraction from 1, we obtain the connectivity ratio.

This is the full version of connectivity function:

The connectivity function provides a crucial numerical measure of the structure’s connectivity, which is essential for evaluating the integrity and print-ability of the 3D object.

Examples

The first shape is hollow in the center (low infill), but has perfect connectivity. In the second shape, every other voxel is filled but only the bottom layer is directly connected to the ground (low connectivity). The third shape is both perfectly dense and connected. The fourth shape has very low infill, but since the ground is defined as the first layer with voxels, technically 100% of its voxels are touching the ground (high connectivity). The fifth shape has high infill, with half of its voxels on the ground and the other half floating.

Optimization

Now that we have a way to both evolve a CA rule and measure its strength/material usage at each step, we are ready to start searching for infill patterns on a large scale. Some CA rules die out immediately before doing anything useful, and promising rules are still only capable of producing a certain type of behavior. Leveraging the diversity and complexity of CA for optimizing infill means we have to search through different rules. The full brute force process involves:

Selecting a CA rule

Generating n evolutions and keeping ones which match infill/connectivity requirements

Incrementing the CA rule and repeating steps n more times

Options

Modern 3D printing software has a plethora of options for how infill is generated in a specific model. In order to make our methods work the best for the largest amount of use-cases as possible, we implemented a modular options system. Also, as a result of defining these fine-grained options, we are able to automatically filter out undesirable infill patterns later in our optimize function. Thus, these options serve to both satisfy user demands for the properties they want their infill to display and “tame” the vast complexity and diversity present in CA by filtering out the vast majority of possibilities we search.

Infill

This option is the equivalent for the “Infill Percentage” option in many slicers and is calculated using the infillPercentage function. This option is set to -1 by default in case the user has no preference for how dense their infill is, triggering a fall through in the optimize function. Since infill percentage directly correlates to how much a user will have to spend on filament material/how long the print will take, we added fine-tuning options. By passing a list in the form of {0.3, 0.1}, for example, a user is telling the function that they want patterns centered around 30% with a margin-of-error of 10%.

Connectivity

This option has no equivalent in a slicer, but directly correlates to the “viability” of a pattern. This option is set to -1 by default in case the user has no preference for how connected their infill is, triggering a fall through in the optimize function. Instead of providing an interval of desired values, users simply pass a proportion for the lower bound of connectivity. 0.95, for example, means that all patterns must have a connectivity percentage above 95%. Setting this value lower may result in more patterns to choose from, but could risk patterns which are weak or need too many “floating” voxels removed in order to print.

Shape

This option determines the bitmask that will be used when generating states within the main loop of the optimize function. Setting this option is equivalent to importing an STL into a slicer. For example, setting this option to “Sphere” will result in patterns that can fit inside 3D models of spheres, “Prism” for rectangular prisms, and “Cylinder” for cylinders. The default option is a prism. Setting this value to “Custom” will activate the “Mesh” option to allow for shapes besides the basic geometric options.

Bounds

This option is only relevant if the user selects “Prism” as their option for the Shape. This option is a positive integer triplet which defines the length, width, and height of the prism in terms of voxels. It is set to {0, 0, 0} by default.

Radius

This option is only relevant if the user selects “Sphere” or “Cylinder” as their option for the Shape. This option is a single positive integer that specifies the voxel radius for sphere and cylinder bitmasks.

Height

This option is only relevant if the user selects “Cylinder” as their option for the Shape. Along with a radius, this positive integer value describes the voxel dimension of the cylinder.

Rule

This option defines the starting point rule for the optimize search algorithm. For example, setting this option to {127, {2, 1}, {1, 1, 1}} would start searching at the 127th rule of this specification (3 dimensional totalistic). Advanced users can also change the rule specification to customize the way each infill pattern will evolve.

Resolution

This option is only relevant if the user selects “Custom” as their option for Shape. Given the custom mesh, this option defines how “finely” to split the mesh into discrete voxels. Higher numbers result in more detailed infill patterns which more closely follow the shape of the mesh. A value of zero would have the least detail possible, and is the default value.

Mesh

This option is only relevant if the user selects “Custom” as their option for Shape. Values for this option must have a Head of type MeshRegion. The overall shape of this mesh is fed into the voxelizer function in order to create a discretized bitmask boundary to contain the infill pattern.

Final Function

We have written our optimize function to find a set of Cellular Automata (CA) evolutions that align with user-specified parameters. This function balances the need for customization with efficiency in generating 3D printing infill patterns.

The optimize function takes three main arguments:

rulesMax: Specifies how many consecutive rules to search after the initial rule.

statesMax: Determines how many times to evolve each rule.

Options: A set of customizable parameters for fine-tuning the infill generation process.

The function works by first generating the correct bitmask based on the specified shape option. It then searches through possible CA evolutions, checking if each pattern meets the infill and connectivity specifications. The total number of evolutions checked is equal to rulesMax * statesMax, which directly impacts the function’s runtime.

To improve performance, the function utilizes parallelization, distributing the workload across multiple cores when available. This approach significantly reduces computation time, especially for complex or large-scale patterns.

Testing

Let’s see if we can find cubic patterns with very low infill percentage yet still maintain high connectivity.

Visual intuition shows us that the first option is the most desirable--other patterns have a low infill percentage because they are simply surrounded by empty space

Now, let’s make cylindrical patterns with higher infill. Let’s also decrease the amount of rules we search from 100 to 50.

Now, let’s use a custom shape.

These examples showcase the function’s ability to generate diverse infill patterns tailored to specific shape requirements and structural needs.

Visualization

In our system, we represent the slices of a 3D CA evolution as a 4D matrix, where the first three dimensions correspond to the 3D spatial coordinates (x, y, z) and the fourth dimension tracks the evolution of the CA across time steps. Each "slice" along the fourth dimension represents the state of the 3D CA at a specific point in its evolution.

The visualization function allows us to explore this 4D matrix interactively. By indexing into the fourth dimension, we can observe how the CA evolves over time, identifying whether a specific rule produces coherent and desirable infill patterns. With these functions, users can visually assess whether the generated pattern aligns with the desired characteristics, such as uniformity, density, and connectivity. The same visualization function integrates with metrics such as density and connectivity, giving immediate feedback on the quality of the generated pattern.

Detecting Unconnected Structures

isolatedParts[]

The first feature we implemented for our visualization code was the ability to single out unconnected structures. Later on, when we integrate this feature into the full visualization function, seeing the isolated parts is crucial for being able to intuitively filter out patterns that either have too many voxels floating (a low connectivity score) or where removing these voxels would cause the pattern to lose major structural integrity.

The isolatedParts function is used to identify and eliminate isolated voxels in the 3D matrix. Isolated voxels are those that are not connected to any other part of the infill structure. This function is done by the same source of code from the connectivity function but only returns the slice of those cells that are not being labeled or flagged. These disconnected elements do not contribute to the overall strength or flexibility of the model, potentially undermining its functional integrity, which is why we need to detect them before the actual printing.

The function is designed to identify and mark isolated components within a 3D structure represented by the input slice. It first checks if the slice is entirely empty (all zeros). If not, it initializes a visit array to track visited cells and a queue to manage the breadth-first search process. The function finds the lowest layer containing a filled cell as a starting point. It then uses the bfsConnectedParts function to perform a breadth-first search, marking all connected cells. The result is a 3D array where connected parts are marked differently from isolated parts, allowing for easy identification of disconnected regions in the structure. This function is crucial for analyzing the connectivity and integrity of 3D printed objects, ensuring that all parts are properly connected.

This function is especially useful for detecting large floating structures that are obscured by superficial detail. Consider this example:

Without identifying the smaller sphere enclosed within the larger one, this pattern would have been impossible to physically print.

connectedParts[]

To convert such patterns to printable versions, we can XOR the original slice and its isolated parts, so only the connected parts will be left.

Below, we can see that even though there are two masses that are connected, they are isolated from each other. We can use the connectedParts function to recognize this:

Only the voxels directly connected to the bottom layer remain.

With help from the isolatedParts and connectedParts functions, we can define our full visualization function, visualize3DMatrix, which will allow us to easily visualize matrix representations of our generated infill patterns, complete with metrics about connectivity.

We start by defining the options for the visualization function, which will be available in the final Manipulate that is outputted. These options are very similar to those that were used earlier for the optimization function:

The visualize3DMatrix function is designed for interactive 3D visualization of infill algorithms. By generating a bitmask based on specified shape options, the function creates multiple slices using cellular automaton rules, producing both regular and isolated part representations. It leverages Wolfram’s Manipulate function to create an interactive interface where users can adjust layer visibility, slice selection, and display modes.

The visualization renders selected voxels as blue cubes within a 3D graphics environment, providing immediate visual feedback on the generated pattern. Users can explore metrics like infill percentage and connectivity for each slice, making it an invaluable tool for researchers and engineers working on 3D printing infill algorithms. The function’s ability to switch between displaying all voxels and isolated parts enables detailed analysis and optimization of cellular automaton-generated structures.

Let's see what a sphere looks like using the default rule. The first argument means that we are only able to view the first 50 evolutions for this CA.

Note: the options we use in the above function make no difference in this cloud version of the notebook because the visualization function is precomputed. To run the full function with all available options, download the notebook.

Export to Physical Models

Let’s go through the entire workflow for CA infill patterns. First, we will specify the desired attributes for our pattern. In this case, we want to print a cube (rectangular prism with equal width, height, length), we want its infill to be in between 25% and 35%, and we want its connectivity to be above 95%

For memory purposes, we will only include a selected example produced by this function call.

Let’s check for any floating voxels before we print:

Nothing! But, remember, if there were unconnected parts we could have ran connectedParts[] to convert the pattern into something we can work with. Now, we are ready to turn this into a mesh and export. We’ll use the built-in Wolfram function ImageMesh, which uses an algorithm called Marching Cubes (the 3D version of Marching Squares). We must first convert our matrix into a graphics Image before applying this function.

Now, we simply use the Export function to turn this into an STL, import into our favorite slicer, drag it inside our cube 3D model, and we’re done!

Here is a screenshot using the BambuLab slicer. Looking inside of the print we can clearly see our pattern at work. Just to see how it stacks up against another common pattern, let’s find the infill percentage.

Next, we have a screenshot of an equivalently-sized and dense (33.7%) cube using a common infill pattern. Notice the regularity in the latter image and the complexity in the former.

For fun, here’s a couple of other meshes generated using our program:

Real World Testing

Aside from our metrics of infill percentage and connectivity (which were incredibly helpful in the brute force rule search stage), we looked to analyze the viability of our approach as a realistic option for cellular infill. In order to determine this, we wanted to make sure that it had similar print times to traditional infill patterns, and was able to print itself without being supported by any external supports on overhangs. To satisfy these conditions, we tested our cellular infill on a simple cube shape, the same one shown in the “Export to Physical Models” section of this notebook.

In terms of print time, our model was comparable to an identical shape using traditional infill at the same infill percentage (similar amount of material used), with both models showing an estimated print time of ~1 hour.

Due to the nature of our infill patterns, and the limitations of 3D printing, we were concerned that the infill patterns would have too much overhang to be self-supported during the printing process. However, this ended up not being the case. Even though our infill had overhanging regions, these overhangs were small enough that the printer could handle them without the need for any extra support structures. Through this finding, we show the viability of our infill pattern.

One result of our cellular infill approach is that the patterns look cool not only on the 3D model, but also on the final printed product. Here’s a picture of the model 30% of the way into the print, with the infill pattern on full display:

The resulting cube we were able to generate was incredibly solid, and yet noticeably light (due to the lower percentage of infill we were able to use without sacrificing on rigidity of the model). Here’s a picture of the final cube:

Impact

Proving the feasibility of Cellular Automata-based infill patterns has major implications for the 3D printing industry. Currently, an average slicer like PrusaSlicer only has about 20 infill patterns available. Only a few of these are responsible for a majority of the millions of prints yearly. Clearly, these patterns are doing a lot of heavy lifting trying to perform the best for a hugely diverse pool of prints. Certainly, there are some prints where these patterns simply underperform, but no better alternatives exist. Using CA immediately changes this dynamic. Now, we’re not trying to come up with one “generalist” pattern that magically works with millions of prints. Instead, we have a million options available, and we just have to filter the one that works best. If a popular pattern only provides desirable strength at >20% infill, a company could use our program to find patterns with 15% infill. By visually filtering candidate patterns and then running a few test prints, it is quite likely that the emergent complexity of CA will have “invented” an alternative pattern that works just as well at 15%, saving the company thousands in printing costs. In this case, using CA freed the company from only being able to use generalist patterns, instead allowing them to sample from beyond human ideas of what a “good” infill pattern looks like to find the truly best option. Creating a set of patterns that works optimally for millions of prints is an impossible task, but it seems like the 3D printing community has accepted this route. Our project shows that improvement in this domain is not likely to come from an expert human discovering another infill pattern with magical properties, but rather from an unthinking computer searching through deterministic processes that, by virtue of their sheer diversity and complexity, are statistically likely to beat generalist patterns in an extremely narrow problem (structurally supporting a specific 3D model).

Further Work

On a broad level, we achieved everything we set out to do. Using Cellular Automata, we were able to generate infill structures, filter possibilities with favorable characteristics, tailor the shape to a specific 3D print, export the infill structure as an STL, and create a physical print using the infill pattern. However, there are still many areas where our current methods struggle, improvements we haven’t implemented, and features that haven’t been fully utilized.

3D Boundary

So far, our only physical test prints have been using basic 3D solids as the boundary. Since the patterns are generated inside of a cube, where a bitmask cuts any “extras” off, shapes which have symmetry around a cube provide the most favorable prints. For example, if a pattern has a central column of support, a sphere would capture that, but a human figurine might not. Given the current resolution we are generating infill patterns, we are not confident that very detailed models would be successfully printed. The solution most likely lies in increasing the resolution of the prints–and thus creating more local complexity that is likely to be included in even asymmetric and detailed models–but we have run into some bottlenecks trying to visualize and export high-resolution patterns. More time would need to be dedicated to seeing if our current software is compatible with more complex prints (including a method for turning an imported STL into a bitmask). However, a modification to our current method might allow for more complex boundaries to be utilized.

Instead of applying a bitmask to the pattern after it has been generated, we should make the CA aware of the boundary during generation. Since CA is iteratively updated, and possesses emergent complexity, we believe that CA is able to respond to obstacles (like a boundary) in an intelligent manner. One possibility would be to map each vertex of a 3D boundary to another random vertex. If a CA cell tries to “escape” at any vertex, it will be regenerated at the corresponding vertex. This mimics the cubic boundary condition already present in CellularAutomaton[], but extends it to much more complex shapes. We hope that this method could allow patterns that are more tailored to their boundary, retain support even with complex boundaries, and require less computation/search overall.

More Metrics

Currently, our only metrics are infill percentage and connectivity. However, due to the extremely diverse nature of cellular automata, these offerings are under-utilizing the potential of the project. For example, while a human might not know how to make a part bend until exactly 20º, there is a good chance that a CA exists which could create a pattern that meets these requirements. The only thing we need to do is create a metric for measuring the properties of patterns in relation to such scenarios. Unfortunately, many metrics more complex than the current ones (e.g. flexibility and strength) would require solid-body simulations to accurately measure. This would greatly impact the performance of our program, since such simulations would have to be run at every iteration of every CA rule in order to check for matches. However, it is possible that we could find proxies for the true metrics, such as analyzing the connections between nodes in order to make assumptions about strength/flexibility.

Slicer Integration

Currently, our project requires users to import the infill STL and the model STL into a slicer and combine them. Also, even though we have metrics such as infill percentage, users still don’t know exactly how much infill their print will require until they load it into a slicer. Turning our project into a slicer plugin could allow us to solve both of these problems while allowing users to stay in one program for the entire experience. For example, the optimization algorithm could communicate with the slicer during the search for infill patterns, querying how much actual material each pattern would require and selecting accordingly. Also, with the ability to import and export STLs, the entire workflow would be integrated, from providing a model for the infill to printing the infill within that model.

Infill Percentage

Currently, the infill percentage is calculated as the ratio of occupied voxels in the pattern to total voxels in the bitmask (e.g. for a given shape, how much is actually filled vs. what could be filled). However, during the first few evolutions of the CA, the pattern might not “expand” to fill the space given by the bitmask. This results in underestimations for the total infill percentage, since the usable portion of the pattern only takes up a small portion of the calculated possible volume. One way to fix this could be to avoid performing metrics until we are certain that the CA has fully filled the bitmask. Another possibility would be to start with a different initial condition. Currently, the source of all the complexity in the patterns is a single filled voxel at the start of the evolutions. However, creating initial conditions that already fill a large volume of the bitmask could decrease the amount of evolutions needed for a large pattern to form. This might also have the effect of producing less symmetric patterns, since they are no longer evolving from a singular point of expansion.

Acknowledgements

This project was created under the Wolfram Emerging Leaders Program (WELP:HS), a 4-month opportunity for students to complete a detailed research program under the guidance of a mentor specialized in their field. Our mentor was Megan Davis, who has been with us ever since we were students at the Summer Research Program (a prerequisite for WELP:HS). They have kept us accountable through all stages of the development process, from creating a Sprint schedule to keeping a clean Git commit history. Also, through their experience working with large Wolfram codebases, Megan helped us all immensely with modularizing our code and working with Wolfram development tools. This transition from individual notebook development (from the Summer program) to package development on a shared repo was a large leap for us all, and we couldn’t have done it without Megan.

Rory Foulger, the Program Director, also provided important guidance to us as we participated in check-in meetings. His positive remarks about our project motivated us to keep a high standard of quality.

CITE THIS NOTEBOOK

Cellular automata generation of 3D printer infill
by Wolfram Education Programs
Wolfram Community, STAFF PICKS, February 4, 2025
https://community.wolfram.com/groups/-/m/t/3379925