Tetrahedron-Based Porous Scaffold Design for 3D Printing

Abstract

Tissue repairing has been the ultimate goal of surgery, especially with the emergence of reconstructive medicine. A large amount of research devoted to exploring innovative porous scaffold designs, including homogeneous and inhomogeneous ones, have been presented in the literature. The triply periodic minimal surface has been a versatile source of biomorphic structure design due to its smooth surface and high interconnectivity. Nonetheless, many 3D models are often rendered in the form of triangular meshes for its efficiency and convenience. The requirement of regular hexahedral meshes then becomes one of limitations of the triply periodic minimal surface method. In this paper, we make a successful attempt to generate microscopic pore structures using tetrahedral implicit surfaces. To replace the conventional Cartesian coordinates, a new coordinates system is built based on the perpendicular distances between a point and the tetrahedral faces to capture the periodicity of a tetrahedral implicit surface. Similarly to the triply periodic minimal surface, a variety of tetrahedral implicit surfaces, including P-, D-, and G-surfaces are defined by combinations of trigonometric functions. We further compare triply periodic minimal surfaces with tetrahedral implicit surfaces in terms of shape, porosity, and mean curvature to discuss the similarities and differences of the two surfaces. An example of femur scaffold construction is provided to demonstrate the detailed process of modeling porous architectures using the tetrahedral implicit surface.

1. Introduction

Tissue Engineering (TE), as an interdisciplinary discipline, is addressed to create artificial, functional three-dimensional (3D) tissues and organs. The superb porous scaffold plays a critical role in cell proliferation and differentiation phases. There are a few basic, widely accepted criteria for building high-quality scaffolds. Firstly, the internal porous structure is required for interconnection and high porosity, and also has proper pore size. Secondly, the materials of scaffolds have to be biocompatible and bioresorbable with a controllable degradation rate. Thirdly, a certain level of mechanical strength is necessary to maintain the structure and support the weight from other parts of the body. Fourthly, a large surface area is of great importance to facilitate the delivery of nutrients and metabolic waste. While replication of the tissue’s geometry is feasible when using advanced computer-aided design (CAD) and medical imaging techniques, the design of internal microstructures with the same bio-mechanical functionality as real tissue is more challenging.

A variety of fabrication methods have been proposed for scaffold manufacturing. Earlier on, Mikos et al. demonstrated a native scaffold design technique using salt leaching [1], while Mooney et al. proposed fabricating porous sponges of poly using gas-foaming [2]. Meanwhile, a few other conventional manufacturing methods, such as thermal-induced phase separation [3,4], porous ceramics [5], electrospinning [6,7], biomineralization [8,9], and phase-separation followed by freeze-drying [10] were also suggested for porous scaffold fabrication. Nevertheless, these methods have limitations with regard to accurately controlling pore morphologies and sizes. With the rapid development of additive manufacturing (AM) techniques and their applications in tissue engineering, a paradigm shift toward computer-aided design has occurred in tissue engineering. Common additive manufacturing (also known as 3D printing) technologies include fused deposition modeling (FDM), selective laser sintering (SLS), precision extrusion deposition (PED), and stereolithography (SLA) [11,12,13,14].

Various commercial CAD softwares, either based on constructive solid geometry (CSG) or boundary representation (B-Rep), have been commonly employed to assist in engineering design and modeling. Some widely used tools are CATIA (Dassault Systèmes, Vélizy-Villacoublay Cedex, France), NX (Siemens PLM Software, Plano, TX, USA), SolidWorks (Dassault Systèmes, Vélizy-Villacoublay, France), and Pro/Engineer (PTC, Needham, MA, USA) [15,16]. These CAD tools combine diverse standard solid primitives (cylinders, spheres, cubes, etc.) through Boolean operations to produce complicated 3D models. Since B-Rep methods describe the solid based on a combination of vertices, edges, and loops, a preliminary check is usually placed to ensure no gaps or overlaps among the boundaries of different cells [17]. In order to enrich the parametric library of scaffolds structures, a set of unit cells with bio-inspired features have also been revealed in some recent research [18,19]. Although the utilization of present CAD tools has largely reduced the time taken for scaffold design and automated the manufacturing process, the biomechanical properties of the resulting scaffolds are still relatively poor. As an alternative solution, a few image-based approaches that combine image processing and free-form fabrication techniques are also used for modeling 3D scaffold geometries. Due to their compatibility with real patient data, image-based methods can quickly create porous architectures by intersecting two 3D binary images, among which one depicts the outline of the defect and the other one consists of a stack of binary unit cells [14]. As an example, Hollister et al. defined the craniofacial scaffold’s topology by specifically setting the voxels’ densities within image design cubes [20]. Later, Smith et al. also used image-based design and CAD tools to generate an accurately sized and shaped scaffold for osseous tissue via selective laser sintering (SLS) [21]. Compared to other CAD scaffold design patterns, implicit surface modeling (ISM) methods allow for complex structures to be easily fabricated using a single mathematical equation. It has high flexibility in describing cellular architectures with the freedom to integrate diverse porous shapes with controlled porosity. One of the most representative ISM methods is the triply periodic minimal surface (TPMS) which has been a topic of great research interest over the past few years. TPMS has been widely used in the replication of natural structures, especially in biological tissues. More details about the TPMS are given in Section 2.

This paper will introduce a new, advanced ISM method using a tetrahedral implicit surface (TIS). A variety of tetrahedral implicit surfaces, including P-, D-, and G-surfaces, are constructed and demonstrated. The major contribution of our work is to address the problem of porous scaffold design from a different perspective based on tetrahedral basic elements, rather than the conventional hexahedrons. The remainder of the paper is organized as follows. In Section 2, we give a brief introduction to TPMSs. Section 3 focuses on the basic mathematical equations of the proposed TIS method and provides a proof to the continuity of TIS surfaces between neighboring tetrahedrons. A flowchart of generating porous scaffolds from a CAD surface model is elaborated on by using the proposed TIS method. In Section 4, various results are compared between different types of TISs and TPMSs. Additionally, some essential properties, such as porosity and mean curvatures, are presented and discussed. Conclusively, a summary of the advantages and deficiencies of TIS is given in Section 5.

2. Triply Periodic Minimal Surfaces

Among various state-of-the-art implicit surface modeling methods, the triply periodic minimal surface has attracted the most attention over the past few decades due to its inherent advantages of geometrical and biological arrangement of pore structure. A minimal surface is locally area-minimizing, which has the smallest area for a surface spanning the boundary of it. The properties of minimal surfaces are very similar to those of a natural biological structure. Many years ago, Galusha proved the existence of minimal surface geometries in intravital biological tissues, including beetle shells and weevils [22,23]. Kapfer also proposed that the TPMS-based sheet solids could provide relatively high stiffness with a high porosity, which creates the possibility of developing structural scaffolds using TPMSs [24]. Later, Rajagopalan first designed coterminous seeding-feeding networks based on Schawartz’s Primitive (P) surface to supply nutrients to the proliferating cells [25]. Melchels et al. used CAD softwares to generate TPMS diamond (D) and gyroid (G) architectures, and introduced the gradient in porosity and pore size by adding a linear term for z-values of the Cartesian coordinate system [26]. Based on these sample designs using a specific TPMS in simple cubes, Yoo summarized and presented a comprehensive CAD porous scaffold design methodology using TPMS libraries composed of diverse primitive units [27].

In mathematics, The Enneper-Weierstrass’s representation describes a global parametrization of TPMS as in [25]:

$$\begin{array}{cc}x& =Re(e^{i\theta }{\int }_{{\omega }_0}^{\omega }(1-{\tau }^2)R(\tau )d\tau )\hfill \\ y& =Re(e^{i\theta }{\int }_{{\omega }_0}^{\omega }i(1+{\tau }^2)R(\tau )d\tau )\hfill \\ z& =Re(e^{i\theta }{\int }_{{\omega }_0}^{\omega }2\tau R(\tau )d\tau )\hfill \end{array}$$

(1)

where $i^2=-1$, $\tau$ is a complex variable and $Re$ means the real part, and $\theta$ is the Bonnet angle. More generally, the approximation of TPMS can be expressed as a single-valued function of three variables from the x-, y-, and z-axis. The surface is the locus of points for which the function has a constant value C, which is defined as the threshold of TPMS. The most common TPMS units are listed in

Table 1. Mathematical definitions of some common TPMS surfaces.

TPMS Surface Type	Mathematical Definitions
P Surface	${\varphi }_P(r)=cos(x)+cos(y)+cos(z)=C$
D Surface	${\varphi }_D(r)=sin(x)sin(y)sin(z)+sin(x)cos(y)cos(z)+$
	$cos(x)sin(y)cos(z)+cos(x)cos(y)sin(z)=C$
G Surface	${\varphi }_G(r)=sin(x)cos(y)+sin(z)cos(x)+sin(y)cos(z)=C$
I-WP Surface	${\varphi }_{I-WP}(r)=2[cos(x)cos(y)+cos(y)cos(z)+cos(z)cos(x)]$
	$-[cos(2x)+cos(2y)+cos(2z)]=C$

.

Although TPMS methods can generate a continuous smooth porous network, the continuity and integrality problems may exist when irregular models directly intersect with the iso-surface. As shown in Figure 1, some isolated malformed curves appear on the boundary of the porous femur. Moreover, it is difficult and time-consuming to partition some complex domains into a set of regular or nearly regularly-shaped hexahedrons. To address these problems, it is necessary to introduce a few advanced mapping techniques to interpolate a regular TPMS matrix into a parameterized model mesh. Yoo proposed a heterogeneous porous scaffold design method by combining the distance field with TPMS [28]. Feng et al. presented a different approach to generating TPMS porous scaffolds based on T-spline [29]. More recently, Chen et al. attempted to map regular TPMS into parameterized hexahedral mesh using the Laplace’s equation [30].

Even though the benefits of TPMS and its applications have been widely discussed in the literature, it might not always be the best tool when manufacturing porous structures. The major limitations of hexahedron-based modeling can be summarized from two perspectives:

(1)

Quadrilateral and hexahedral meshes are usually less frequently and efficiently used than triangular and tedrahedral meshes in various image-processing areas.

(2)

More importantly, a hexahedron will introduce a distortion issue when deformed in 3D which is not seen in a tetrahedron.

3. Tetrahedral Implicit Surface Modeling

3.1. Mathematical Description of Tetrahedral Implicit Surface

TPMS is characterized as a minimal surface which is periodic in three perpendicular directions. Obviously, the same manner does not apply to the generation of periodic surfaces in a tetrahedron. Our goal is to make the surface periodic in certain particular directions, such as the vertical direction of each of the four faces of a tetrahedron. Here, we introduce a new coordinates system using point-to-plane distances. In this system, there are four axes, corresponding to the vertical directions of four faces of the tetrahedron. The coordinate of an arbitrary point inside the tetrahedron on a certain axis is determined by the ratio of the distance from the point to the corresponding face divided by the corresponding height of the tetrahedron. As shown in Figure 2, the coordinates of point P, ($\alpha$, $\beta$, $\gamma$$\sigma$), can be calculated by:

$$\begin{array}{cc}{P}_{\alpha }& ={\displaystyle \frac{d_{\alpha }}{h_{\alpha }}}\times 2\pi \times {\lambda }_{\alpha }\hfill \\ P_{\beta }& ={\displaystyle \frac{d_{\beta }}{h_{\beta }}}\times 2\pi \times {\lambda }_{\beta }\hfill \\ P_{\gamma }& ={\displaystyle \frac{d_{\gamma }}{h_{\gamma }}}\times 2\pi \times {\lambda }_{\gamma }\hfill \\ P_{\delta }& ={\displaystyle \frac{d_{\delta }}{h_{\delta }}}\times 2\pi \times {\lambda }_{\delta }\hfill \end{array}$$

(2)

in which $d_{\alpha }$, $d_{\beta }$, $d_{\gamma }$, and $d_{\delta }$ denote the projection distances from point P to plane $BCD$, $ACD$, $ABD$, and $ABC$, respectively, and $h_{\alpha }$, $h_{\beta }$, $h_{\gamma }$, $h_{\delta }$ are the height of the tetrahedron corresponding to each plane. The parameter ${\lambda }_i(i=\alpha ,\beta ,\gamma ,\delta )$ determines the periodicity of its corresponding direction.

To mimic similar properties of TPMS, we adopt the mathematical expressions of the TPMSs on each periodic component and derive the following TIS equations:

$$\begin{array}{cc}{\varphi }_P(r)=& cos(\alpha )+cos(\beta )+cos(\gamma )+cos(\delta )\hfill \\ & +cos(\alpha )cos(\beta )cos(\gamma )cos(\delta )=C\hfill \\ {\varphi }_D(r)=& sin(\alpha )cos(\beta )cos(\gamma )cos(\delta )+cos(\alpha )sin(\beta )cos(\gamma )cos(\delta )\hfill \\ & +cos(\alpha )cos(\beta )sin(\gamma )cos(\delta )+cos(\alpha )cos(\beta )cos(\gamma )sin(\delta )\hfill \\ & +cos(\alpha )cos(\beta )cos(\gamma )cos(\delta )=C\hfill \\ {\varphi }_G(r)=& cos(\alpha )sin(\beta )+cos(\gamma )sin(\delta )+cos(\alpha )sin(\beta )\hfill \\ & +cos(\gamma )sin(\delta )+cos(\alpha )cos(\beta )cos(\gamma )cos(\delta )=C\hfill \end{array}$$

(3)

It is worth noting that these TIS surfaces are expressed by implicit functions with a constant threshold C. With a proper threshold given, the TIS surface mesh can be generated as an iso-surface by using the Marching Cubes (MC) algorithm, which is one of the most famous approaches to extracting iso-surfaces from a three-dimension data field.

3.2. Continuity Analysis of Tetrahedral Implicit Surface

As mentioned in Section 2, TPMS couldn’t produce a porous structure fitting any irregular complex shape unless some mapping techniques, such as the distance field [31], T-spline [29], or hexahedral parameterization [30], is integrated. Unlike TPMS, a significant advantage of the proposed tetrahedral implicit surface method is that it can generate surfaces with same morphological features inside the tetrahedron of any shape. When multiple tetrahedra are considered, which is almost always true in modeling real-world shapes, the surfaces between neighboring tetrahedrons must be seamlessly connected at the joint. Otherwise, it may cause holes or gaps in the resulting surfaces. In this section, a mathematical geometry proof is provided to validate the $C^0$ continuity of TIS between two adjacent tetrahedrons.

As shown in Figure 3, $ABCE$ and $ADCE$ are two adjacent tetrahedrons of arbitrary shapes that share the triangular face $ACE$. F is an arbitrary point on the face $ACE$. I and J are the perpendicular projections of the point F onto the two faces $BCE$ and $DCE$, respectively. Similarly, H and K are the projections of the point A onto the two faces $BCE$ and $DCE$, respectively.

$$\begin{array}{c}becauseFI\perp BCEandAH\perp BCE\\ therefore,FI\parallel AH\\ and\Delta GFI\sim \Delta GAH\\ soweget:r_1=\frac{GF}{GA}=\frac{GI}{GH}=\frac{FI}{AH}\end{array}$$

(4)

The same holds for tetrahedron $ACDE$ and we have:

$$r_2=\frac{GF}{GA}=\frac{GJ}{GK}=\frac{FJ}{AK}$$

(5)

Combine Equations (4) and (5), then:

$$r=r_1=r_2=\frac{GF}{GA}=\frac{FI}{AH}=\frac{FJ}{AK}$$

(6)

In tetrahedron $ABCE$, $\frac{FI}{AH}\ast 2\pi$ and $\frac{FJ}{AK}\ast 2\pi$ is the coordinate of point F with respect to the faces $BCE$ and $CDE$, respectively. Similarly, we can conclude that the other two coordinates of point F in tetrahedron $ABCE$ with respect to faces $ABC$ and $ABE$ are the same as those of point F to faces $ADC$ and $ADE$ in tetrahedron $ADCE$. Meanwhile, the coordinate of F to face $ACE$ equals to 0 in both tetrahedrons $ABCE$ and $ADCE$, which means the four coordinates of point F are exactly the same in both of the two tetrahedrons. As a result, we can conclude that the TIS function values of all points lying on the shared face of two neighboring tetrahedrons will be the same, regardless of their shape.

3.3. Porous Structure Generation from Surface Moldes

Based on the proposed TIS method, we can achieve porous scaffolds by tetrahedralizing the tissue with a customized configuration. Biological tissues are inherently heterogeneous architectures. At the macrostructure level, the tissue exhibits great differences in both morphology and biofunctionability. Pore-intensive parts have higher structural strength and can support most of the weight. Conversely, the sparsely populated parts have better permeability and transport of nutrients and metabolic waste. Figure 4 shows a flowchart of the details for designing a 3D heterogeneous porous scaffold from the 3D anatomical shape.

Since the quality of the tetrahedral mesh, including the face angle and dihedral angle, has a decisive influence on the quality of the generated scaffold, our approach uses a tetrahedral mesh generation algorithm based on the body-centered cubic (BCC) tetrahedral lattice [32], followed by quality tetrahedral mesh smoothing via boundary-optimized Delaunay triangulation (B-ODT). Compared with the two current popular methods, Tetgen and Netgen, the BCC lattice tetrahedral meshing algorithm is more robust and general and has no quality requirement for the input surface mesh, and the sharp features are also reconstructed better. To tetrahedralize the input model, the octree of the input was firstly subdivided into 3D mesh based on the Euclidean distance transformation. Then, we computed the sign of each node in the BCC lattice and calculated the intersecting points where the edge crossed the input surface mesh. Next, we detected the cutting points that were “too close” to the original BCC nodes and snapped them to the corresponding nodes. Finally, the boundary polyhedra was decomposed into tetrahedrons and the entire generated tetrahedral mesh was post-processed by the B-ODT smoothing algorithm. For more details on the tetrahedral meshing algorithm, please refer to [33,34].

4. Results and Discussion

In this section, various types of TIS surfaces are compared with the TPMS surfaces from different perspectives. The TIS algorithm is implemented in C programming language running on an 8 × 2.20-GHz Intel(R) Core(TM) i7 computer with 16 GB of memory. A volume of $200\times 200\times 200$ voxels was used as the resolution for both TPMS and TIS to generate various types of surfaces. While different selections of the threshold C can significantly change the porosity and pore size, we chose the constant 0 as the threshold of the generation of P-, D-, and G-surfaces for all TPMS models in this paper. For TIS surface generation, however, we used a concept of relative threshold, which is related to the actual threshold as follows:

$$actual\_threshold=(max-min)\ast relative\_threshold+min$$

(7)

in which $max$ and $min$ are the minimum and maximum of the calculated TIS function values from Equation (3). The range of the relative threshold was rescaled into 0 to 1 according to the actual threshold, which is the real threshold used in the experiments. The periodicity parameter $\lambda$, as used in Equation (2), was chosen as follows. To stereoscopically illustrate the fabricated TIS surfaces with a moderate pore size, $0.16$ and $0.44$ were used as the thresholds for P-surface ($C=0.16$ when $\lambda =1$, $C=0.44$ when $\lambda =2$ for the generation of TIS surfaces). Also, we used the threshold of $0.67$ for the TIS D-surface and a threshold of $0.5$ for the TIS G-surface both for $\lambda =1$ and $\lambda =2$.

In Figure 5, Figure 6 and Figure 7, the surfaces in the two left columns are TPMS Schwarz P-, D-, and G- surfaces, while the right two columns contain the TIS surfaces of the same type, respectively. The meshes in the first row were constructed using a single period ($\lambda =1$), while the meshes in the second row were constructed using a double period ($\lambda =2$) on each axis instead. From these figures, it is not hard to figure out that the shapes of the surfaces and the locations of the holes are very similar between the same types of surfaces for different groups. For example, the P-surfaces in both TPMS and TIS have an opening on each face of the unit domain (either tetrahedron or hexahedron). The opening on the tetrahedral face is an ellipse of the triangle outline, while the opening on the hexahedral face is a regular circle. The interiors are also connected, which allows for cell growth through the holes on the surface. In addition, Figure 8, Figure 9 and Figure 10 also give closer views of the surfaces between two neighboring tetrahedrons.

To further analyze the characteristics of these tetrahedral implicit surfaces, we first used the well-known open source CAD software, Meshlab, to compute the average curvatures over the whole surface. Figure 11 demonstrates the distributions of curvatures of different surfaces, as well as their corresponding detailed histograms. The color map ranges from red to blue, indicating the curvature values from a minimum of $-30$ to a maximum of 30 (the values on the X-axis are curvatures in ascending order, while the values on the Y-axis are the number of points whose curvatures share the same value). Additionally, more information including the minimum/maximum curvature and the mean curvature, as well as the standard deviation of the curvatures are also provided for each of the TPMS and TIS surfaces in

Table 2. Statistical analysis to the curvatures of TPMSs and TISs.

Curvatures	TPMS Surfaces			TIS Surfaces
	P-Surface	D-Surface	G-Surface	P-Surface	D-Surface	G-Surface
Minimum	−12.66	−14.99	−14.34	−33.19	−14.88	−26.48
Maximum	12.59	13.83	18.32	21.09	12.16	15.17
Medium	−0.57	−1.18	−0.66	−1.68	0.45	−1.68
Average	−0.49	−1.29	−0.79	−2.88	−0.94	−2.79
Standard Deviation	3.78	3.32	2.85	4.11	4.32	4.51
Variance	14.28	11.04	8.11	16.88	18.63	20.36

. In the table, the minimum curvatures of the TIS P-surface and G-surface are almost double the values of the TPMS P-surface and G-surface. Only the maximum curvature of the TIS P-surface is significantly larger than the TPMS surface of the same type. Moreover, the standard deviation of the curvatures of TIS surfaces is relatively larger than that of the TPMS surface of the same type. From Figure 11 and Table 2, it is not hard to conclude that the distributions of curvatures of TPMSs are all symmetrical, while this symmetry does not apply to TISs. Nonetheless, one of the most important common attributes for TPMS and TIS is that the overall distribution of the curvatures tends to converge to the mean curvature, which approximates to 0.

Besides, we also realized that the porosity and pore size of the surfaces were quite sensitive to the selection of the threshold. Through the experiments, we found that not all thresholds in the range of 0 to 1 could generate a completely connected surface. To create the TIS surfaces of a single period from a regular tetrahedron, the ideal threshold range for a P-surface is 0.118–0.275, for the D-surface is 0.569–0.765, and for the G-surface is 0.431–0.627. A value out of these ranges will cause the generated corresponding type of the deformed surface, which means the hole structure will disappear. On the other hand, TPMS surfaces also have the same issues of surface discontinuity and deformation as TIS surfaces. Nonetheless, the difference is that the selective areas of the reasonable threshold for TPMSs all uniformly concentrate on the range from $-1$ to 1.

Table 3. Effect of diverse thresholds C on different TPMS and TIS surfaces.

TPMS P-Surface		TPMS D-Surface		TPMS G-Surface
Threshold C	Porosity (%)	Threshold C	Porosity (%)	Threshold C	Porosity (%)
−0.95	22.79	−0.95	10.28	−0.95	18.77
−0.5	35.71	−0.5	29.42	−0.5	33.81
0	49.98	0	49.94	0	50.00
0.5	64.28	0.5	70.57	0.5	66.19
0.95	77.21	0.95	89.71	0.95	81.22
TIS P-Surface		TIS D-Surface		TIS G-Surface
Threshold$\mathit{C}$	Porosity (%)	Threshold$\mathit{C}$	Porosity (%)	Threshold$\mathit{C}$	Porosity (%)
0.118	18.36	0.569	57.35	0.431	44.03
0.157	35.80	0.618	66.88	0.480	57.29
0.196	49.36	0.667	75.73	0.529	67.76
0.235	61.17	0.716	83.04	0.578	76.09
0.275	71.86	0.765	89.90	0.627	84.04

gives more information of the relationship between measured porosities and different thresholds, C. Additionally, some examples of the deformed TPMSs and TISs with a threshold that exceeds the proper range are demonstrated in Figure 12. It is noted that the hole size gets bigger and the porosity behaves more saturated when C gets increased in the given range for both TPMS and TIS. Meanwhile, the surface areas of different TPMS and TIS surfaces for a unit volume are measured by the sum area of the generated mesh triangles, and the results are presented in

Table 4. Surface areas of TPMSs and TISs for a unit volume.

Surface Type		Surface Area for Unit Volume
TPMS Surfaces	P-Surface	2.353180
	D-Surface	3.828864
	G-Surface	3.083422
TIS Surfaces	P-Surface	3.798055
	D-Surface	2.367160
	G-Surface	2.654556

as a reference. All related parameters, including the resolution and threshold, are the same as mentioned at the beginning of Section 4. In the end,

Table 5. Running time comparison between TPMS and TIS.

Resolution	TPMS P-Surface		TIS P-Surface
	Computational Time (s)	Number of Triangles	Computational Time (s)	Number of Triangles
100	0.108	74,216	0.083	27,379
200	1.212	299,672	0.674	112,381
300	3.895	675,848	2.930	254,074
400	9.749	1,203,704	7.361	454,630
500	18.524	1,882,376	14.708	708,986

lists the running time of the TPMS and TIS algorithms, as well as the number of generated triangles for the P-surfaces with different resolutions. The results indicate that the computation time for the TIS algorithm is relatively shorter and mesh size is much smaller than the TPMS when the same resolution is used.

Numerous experimental results implied that the TIS P-surface performs better in the construction of large models than the D- and G-surfaces. Below, we use the P-surface as an example to elaborate the processes of the generation of a porous femur scaffold. The input model is a 3D femur mesh composed of 7587 vertices and 15,174 faces as shown in Figure 13. Through the tetrahedral meshing operation, the femur model is firstly tiled into 989 tetrahedrons, in which over 95 percents of triangular angles lies in the range between 50 degrees and 100 degrees. Figure 14a shows the surface mesh of the tetrahedralized femur model and Figure 14b presents the cross-sectional view of the generated tetrahedral mesh. Secondly, by calculating the TIS function values for all grid points inside the tetrahedral network, we can easily extract the iso-surface through the Marching Cubes algorithm using the threshold $C=0.176$. Moreover, we applied Laplacian smoothing and curvature-based smoothing on the entire porous surface mesh to improve the $C^1$ continuity at the joints. Figure 15 depicts the smoothed porous TIS surface constrained by the input femur’s surface mesh. Lastly, we intersected the original input femur model with the constructed porous TIS volume to produce a porous femur scaffold, as shown in Figure 16.

5. Conclusions and Future Work

In this paper, an implicit 3D surface method for designing bio-scaffolds based on the tetrahedral implicit surface was proposed. To the best of our knowledge, this is the first attempt to construct porous structures based on tetrahedral elements. The TIS method demonstrated its great potential in the field of designing optimized porous structures. Compared to the well-established TPMS approach, this concise method brings several advantages to the tissue engineering scaffold design:

(1)

Generation of the TIS surface could avoid the restriction of modeling based on a regular unit domain, such as the regular hexahedron in TPMS. Even in the deformed tetrahedron, there can still be the creation of a characteristic tetrahedral surface, which can guarantee $C^0$ continuity between adjacent TIS surfaces.

(2)

Different from the TPMS method that needs to parametrize tissue surfaces into hexahedral meshes, it is more convenient to tetrahedralize the triangular surface mesh as in the proposed TIS method. Without the special mapping procedure, the entire process of generating pores is also simplified.

(3)

The strong interconnectivity and tetrahedron-based modeling grants the TIS more flexibility and creativity for complicated shape modeling.

Notwithstanding this, there are also several limitations of this method which needs further research. Though the TIS surface of multiple periods obtained from a single unit has shown the characteristic of approximately zero mean curvature, and though we also showed the $C^0$ continuity of the generated surfaces, the $C^1$ continuity of the surface at the junction of two neighboring tetrahedrons is not guaranteed and needs further investigation. Currently, the mesh smoothing algorithm is used to greatly improve the smoothness of the resulting surface. Moreover, some isolated closed surfaces may be produced in the TIS D-surface from a complex input model.