Fractal algorithm for 3D‐printed continuous porous scaffold design

In practice, 3D modelling software, such as UG and Pro/E, are commonly adopted to design porous scaffolds within given contours. However, this manual method is quite time-consuming and complicated with poor adaptability. In this study, a novel facial algorithm is proposed for smart scaffold design. By forming a dendritic fractal network, the second-order uniformity and complete connectivity of pores are guaranteed while considering the axial symmetry of certain orthopaedic implants, and the algorithm is adapted for concave and hollow contour shapes. Through experiments, the stability of the algorithm and the characteristics of the two vital performance indicators, i.e. porosity P and surface area ratio R, are studied. The results show that P and R of the generated model are insensitive to the shape of the input model, which proves the stability of authors’ method, and the maximum available porosity reaches 78%, higher than the maximum effective porosity found in the literature.


Research background
During past decades, mechanical fixation methods had long been used to fix orthopaedic implants to bone tissue, which brought instability of the implant-bone tissue contact surface and detachment and peripheral fractures of orthopaedic implants [1]. In recent years, studies have found that implants made of certain porous materials, also known as tissue scaffolds, can allow remarkable stability due to the biological binding (i.e. synostosis) between the host bone tissue and the implant, by increasing the porosity, connectivity, and surface area ratio of the scaffold [2][3][4][5][6].
With the development of additive manufacturing (3D-printing) technology, its layered construction characteristics and highprecision features make it superior to conventional scaffold manufacturing methods (e.g. salt infiltration); thus, it has been gradually accepted for clinical usage, mainly with commercial 3D modelling software, such as UG and Pro/E. However, not only the software are difficult for doctors to pick up with, but also the algorithms commonly used by them are not flawless as well.
Methods of porous scaffold design should be roughly divided into two categories: (a) Accumulating certain structural units to form a meshed skeleton, such as various isotropic structures [7], adaptive optimisation stacking method [8], dual-cell stacking method [9], anisotropic stacking method [10], etc. (b) Constructing the continuous axis in a certain way and to generate pipelines on axes, such as the continuous deposition method proposed by Khoda et al. [11].
The problems are demonstrated as follows: as for the methods of simple isotropic stacking, the approximate axial symmetry, which commonly exists among bone implants, is not considered.
Regarding the methods of anisotropic stacking, the connectivity of the implant pores is gradually improved while it cannot guarantee that the implant has one and only one multi-connected area (i.e. there may be pores that are unconnected to the outside, which are impossible for tissue to grow into), due to the inherent contradiction between the (distant) invariance of structural elements and the variance of stacking pattern. In addition, Khoda's algorithm is more sophisticated, yet there are still some flaws with its adaptability to concave/hollow implants due to that it is based on an axis-aimed convergence method. Moreover, several aforementioned methods do not optimise surface area ratio. Therefore, it is worthwhile to design and implement a scaffold modelling algorithm that avoids these problems.

Problem statement
The algorithm input is the contour of the porous scaffold defined by an imperforated additive manufacturing model, shown in Fig. 1. In addition, the algorithm accepts several control arguments to construct the dendritic network and the scaffold, such as the maximum iteration order, the radius of pipelines, pipeline shape, etc. For the sake of clarity, the three most important arguments maximum iteration order, number of dividing points per cell, and number of segments per point will be referred to as I, J, K . Available output includes the porous scaffold (stored in a stereolithography file), slices in three directions, log recording performance indicators, and other variables. Fig. 2a shows the destination scaffold constructed from the contour, shown in Fig. 1, while Fig. 2b shows slices obtained from it.
The indicators that measure the effectiveness of this algorithm include porosity, surface area ratio, uniformity, adaptability, stability, and so on.

Process synthesis
The algorithm mainly operates as follows: with the contour given, the algorithm obtains its central axis and generates a dendritic fractal network within its bounding box by iteration and linear mapping method. On the dendritic network, a connected cylinderpolyhedron skeleton is formed; thereby it is converted into channels within the given contour by a Boolean method to generate the destination scaffold. A more specific flowchart is shown in Fig. 3. In the remainder of this section, the main parts of the algorithm are illustrated and discussed in terms of logical sequence. In addition, unless otherwise specified, all example pictures are generated from the same set of contour and arguments.

Pretreatment:
The algorithm's core is building the fractal network inside the input model. However, if the triangle mesh surface of the model itself is directly used as a boundary condition in the fractal network generating process, two problems will arise: (a) Such implant model is often concave or even hollow, i.e. the intersection(s) of any straight line intersecting the model surface might not be the only one or two, which leads to ambiguity in boundary conditions. (b) Even if the boundary ambiguity is minimised through technical means, when an irregular triangular mesh is used as the boundary to constrain the fractal network formation, it is necessary to calculate and use normal vectors for each triangular mesh to perform operations, which is not economical on neither time nor space complexity. In addition, this approach still requires Boolean operations to ensure that internal space of the model has a sufficient number of pipelines that connect to outside, so the former mentioned computational demands cannot be eliminated elsewhere.
Therefore, it is necessary to adopt an appropriate method to obtain a convex bounding surface that is easy to apply intersection operation on and close to the model surface, and to use it as an initial constraint for generating the fractal network. Such method is often referred to as finding the bounding box, and commonly used bounding boxes include the axis-aligned bounding boxes (AABB) (minimum hexahedrons surrounding the model with sides parallel to the coordinate axis), the oriented bounding boxes (AABB without the parallel constraint), and the fixed directions hulls. Although the latter two are closer to the model surface, both theory and experiments have shown that their complexity makes them unsuitable for our modelling algorithm [12], whose boundary conditions required a high demand to be simple and robust.
Before constructing the fractal network, it is also necessary to measure several geometric parameters of the input model, such as size, pointing, centre position, surface area, volume, and so on. Among them, the surface area and the volume data are obtained by an algorithm based on the divergence theorem proposed by Alyassin et al. [13].

Obtaining medial axis:
In order to iterate and generate the dendritic network, a 0th order segment set must be picked as the entrance and then stored into two datasets, CS and NS.
Naturally, due to its topological equivalence, the curve made-up of connected mid-axis segments is commonly picked (often with a layer-by-layer contraction method) for similar purposes [11]. However, such an approach is not suitable for the proposed algorithm for two main reasons: (a) The layer-by-layer contraction method meets difficulties while handling concave or hollow implants. (b) When fractal polylines are radially expanded, the curvature of the 0th order curve will increase the extent of the overlap of higher-ordered pipelines and reduce the destination porosity and symmetry.
Therefore, a more feasible way is to linearly fit the medial axis equation in the pointing direction of the model. Moreover, consider that the difference in medial axes caused by the variance of the number of fitting points is only significant for weak axisymmetric implants that are insensitive to the position variance of the 0th order axis, the fitting process is better to be simplified by fetching centre points of the two opposite faces in the model pointing direction with the bounding box, and such a hypothesis is confirmed by experiments.

Building fractal network:
After the medial axis segment is picked, the iterator traverses cells stored in the CS dataset for I + 1 times. In the ith step, on each ith-ordered segment, l i , a set of dividing points P i j is determined, which contains a total of J points. Moving forward, on each dividing point P i j there are K segments uniformly omitted on the plane perpendicular to l i , all marked as i + 1 th-ordered. All the former mentioned elements are appended to the NS dataset and along with each traversal, a NS − CS copying procedure is operated to enable the next round. Such a process has considered the approximate axial symmetry of some orthopaedic implants, yet to meet other needs J ⋅ K segments l 1 , …, l JK [ i + 1 th order] must be generated on l while guaranteeing three vital points: (a) l 1 , …, l JK are perpendicular to l. (b) l 1 , …, l JK are generated only inside the bounding box and evenly distributed around l (the first-order iteration constraints). (c) The distribution of l 1 , …, l JK ensures that all the i + 2 th segments generated by the next iteration step are still evenly distributed around l (the second-order iteration constraint), to reduce the overlap of subsequent scaffolds and make rational use of the space.
In order to satisfy these conditions, a linear mapping is used to convert the emitting operation to a two-dimensional problem in a relatively computational-friendly Euclid space. Let AB denotes the vector from the starting point to the ending point of l and e →S denotes the vector base used to define the world coordinate system. Under e →S , the coordinate matrix of AB is denoted by Define the vector base e →S′ that satisfies two properties: (i) The origin of e →S′ which is point A.
Under these definitions, the linear mapping ℳ: S′ → S that transforms the world coordinate system S to the new Euclid space S′ is determined, and ℳ is equivalent to a combination of the rotational ℛ and translational T transformations. Let θ and c denote the rotation angle and axis pointing matrix of ℛ, and t denotes the translation coordinate matrix of T, then Hereby, AB is equally divided in the S′ space to obtain J + 1 points (including the two endings). At this time, to ensure the establishment of condition (a), the problem only needs to be considered on the plane family N = N j : x = x j while emitting rays vertical to l at dividing points P 0 x 1 , 0, 0 T , …, P j x j , 0, 0 T , …, P J 1, 0, 0 T , i.e. the simplification of dimension reduction is achieved. For the convenience of narrative, the polar coordinate system ρ, φ is established on each plane, all their polar axes are along the direction of the vector denoted by 0, 1, 0 T in S′ space.
To ensure condition (b), ray family r j = r jk is emitted on P j to divide plane N j to K equal parts. In the polar coordinate system origins at P j , r j should be denoted by To ensure condition (c), the value of the above-mentioned initial angle φ j should be appropriately selected as Through this process, the ray network represented in the S′ space is obtained and mapped into the world coordinate system with ℳ.
After which the intersection of each ray and the bounding box is obtained, and segment cells, each one is constructed with one intersection point and the corresponding dividing point in P 0 x 1 , 0, 0 T , …, P j x j , 0, 0 T , …, P J 1, 0, 0 T , are marked as i + 1 th order and added to the NS dataset. At this point, a fractal network is generated inside the bounding box. Fig. 4 shows two structures with different arguments.

Building scaffold:
To build the pore area inside the scaffold, cylinders and polyhedrons are, respectively, built upon segments and nodes of the fractal network and then stored as the trianglepatched closed body. Taking into account the previous studies on the topic of suitable circumstance for bone ingrowth, a hexagonal prism with a 300-µm radius is picked for default usage [14,15]. In addition, the polyhedron radius is slightly larger than the cylinder radius by default to ensure connectivity at channel nodes. Fig. 5 shows what the cylinder-polyhedron skeleton should look like. As there still are boundary edges in the skeleton model, a mending filter provided by the Visualization Toolkit is applied to fill the possible loopholes. Eventually, the normal vector of every mesh on the skeleton surface is calculated and a Boolean operation [16] is then used to convert the skeleton to pipelines inside the input model. At this point, the destination scaffold is generated.

Experiments
Consider that the uniformity of pores and adaptability to concave models are guaranteed theoretically, this section aims to clarify and study the stability and other indicators of the algorithm.

On stability
In this section, the word stability is defined as the correlation between control arguments and outputs, e.g. if we abstract our algorithm and quantify it with the indicator of porosity, it should be like P = P A P , I M = P I, J, K, …, I M (4) where A P stands for input arguments directly entered by users and I M stands for the information gathered from the contour model. In order to improve the universality of the algorithm, the correlation between P and I M should be minimised and the nature of the algorithm would then mostly characterise the properties of the scaffolds -porous material that made-up of the implant -itself, rather than the shape or size of the implant. Additionally, the input arguments can be mostly referenced when generating scaffolds with the desired performance indicator. The experiment is designed to study the stability on R, P while changing the configuration of the models. Consider that several control arguments are extensive amounts to input model size, for each set of input parameters the model size is fixed.
The experiment operates as follows: Two sets of the fractal network with control arguments A 3, 4, 4 , B 3, 3, 6 are picked, and pipeline radius and shape are set to default. The input contour models are polyhedrons with resolutions of 4-9 for latitude and 2-9 for longitude, each one with a 5 mm radius. Obtain two sets of samples of the surface area ratio R and the porosity P.

On performance indicators
The fractal network constructed in this study is an exponential growth structure. For a network with arguments I, J, K , the number of segment cells that represent the complexity is Despite that the number of cells in this less-than-two-dimensional network can increase indefinitely, the space in the given contour is indeed limited; hence, blindly increasing N will cause overlap and variance of actual pipeline radius. It can be expected that there is a critical point, depending on the input model and arguments, which indicates the appearance of overlap. Consider that the porosity should not be, as tested in Section 3.1, a variable to geometric size while N is one, the critical point should be represented by the value of porosity in order to study properties of the scaffolds itself. Hereby, the critical porosity P cr equals the maximum porosity that scaffolds can have without affecting the pipeline radius. The experiment is designed for three purposes: (a) To find out whether such critical point exists and, if it does, the value of the maximum available porosity P cr . (b) To study the relationship between R, P and the complexity of the fractal network in a given contour.
(c) To give guidance to the appliance of the algorithm.
Therefore, the experiment operates as follows: A regular icosahedron with a radius of 5 mm is selected as the input model, and input arguments are set to I, J, K I = 2, 3, and 4 ≤ JK ≤ 18 with default pipeline shape and radius. Calculate N with 4 and obtain P, R from the program.

On stability
Assume that P and R are normally distributed, using the Student's distribution to estimate overall distribution parameters of the two sets of samples with a 95% confidence interval, as presented in Table 1.
The result indicates that the stability of R and P is very strong while changing the model shape under given network arguments and model size. Therefore, when generating scaffolds with the desired porosity, the network arguments only need to be determined with reference to the model size without paying much attention to model details.

On performance indicators
The obtained R − N and P − N scatters (destinations) and the fitted lines are shown in Fig. 6. (a) When the porosity is lower than 0.788, there is a good linear relationship between porosity and network complexity, and the partial derivative of porosity on the number of segments is 6.2 ±6 × 10 −4 . The situation corresponds to a large channel spacing and little overlap. (b) When the porosity is higher than 0.788, the partial derivative of the porosity to the number of segments mutates to 3.4 ±9 × 10 −5 .
The reduction in the magnitude of the partial derivative strongly indicates the overlap of channels, changes of radii, and appearance of isolated parts; thus, P cr is inferred to be 0.788. Considering that the maximum effective porosity 70% is lower than P cr , this maximum value meets practical needs.

Prospect
One feasible prospect can be the development of a novel and efficient Boolean operator to deal with triangle patches because the operator used in Section 2.2.4 leads to the consumption of plenty of computation time and memory space, sometimes even crashes. In addition, it does not work well when dealing with sharp cuts. This may cause a problem regarding the robustness of the algorithm.
Another prospect is to conduct an experiment on real additive manufactured scaffolds and actual tissue. Before this type of experiment is run through, a mechanical strength test should be conducted in advance.