Heterogeneous porous scaffold generation in trivariate B-spline solid with triply periodic minimal surface in the parametric domain
HHeterogeneous porous sca ff old generation in trivariate B-spline solid with triply periodicminimal surface in the parametric domain Chuanfeng Hu a , Hongwei Lin a,b, ∗ a School of Mathematics, Zhejiang University, Hangzhou, 310027, China b State Key Lab. of CAD & CG, Zhejiang University, Hangzhou, 310027, China
Abstract
A porous sca ff old is a three-dimensional network structure composed of a large number of pores, and triply periodic minimalsurfaces (TPMSs) are one of conventional tools for designing porous sca ff olds. However, discontinuity, incompleteness, and highstorage space requirements are the three main shortcomings of TPMSs for porous sca ff old design. In this study, we developed ane ff ective method for heterogeneous porous sca ff old generation to overcome the abovementioned shortcomings of TPMSs. The inputof the proposed method is a trivariate B-spline solid (TBSS) with a cubic parameter domain. The proposed method first constructsa threshold distribution field (TDF) in the cubic parameter domain, and then produces a continuous and complete TPMS within it.Moreover, by mapping the TPMS in the parametric domain to the TBSS, a continuous and complete porous sca ff old is generatedin the TBSS. In addition, if the TBSS does not satisfy engineering requirements, the TDF can be locally modified in the parameterdomain, and the porous sca ff old in the TBSS can be rebuilt. We also defined a new storage space-saving file format based on theTDF to store porous sca ff olds. The experimental results presented in this paper demonstrate the e ff ectiveness and e ffi ciency of themethod using a TBSS as well as the superior space-saving of the proposed storage format. Keywords:
Porous sca ff old, Trivariate B-spline solids, Triply periodic minimal surfaces, Parametric domain
1. Introduction
Porous structures are widely found in natural objects, such astrabecular bones, wood, and cork, which have many appealingproperties, such as low weight and large internal surface area.In the field of tissue engineering, which aims to repair dam-aged tissues and organs, porous sca ff olds play a critical role inthe formation of new functional tissues for medical purposes,i.e., they provide an optimum biological microenvironment forcell attachment, migration, nutrient delivery and product ex-pression [1]. To facilitate cell growth and di ff usion of both cellsand nutrients throughout the whole structure, high porosity, ad-equate pore size and connectivity are key requirements in thedesign of porous sca ff olds [2]. Therefore, it is important to beable to control the pore size and porosity when designing het-erogenous tissue sca ff olds.Recently, triply periodic minimal surfaces (TPMSs) havebeen widely employed in the design of porous sca ff olds [3]. ATPMS is a type of minimal surface with periodicity in three in-dependent directions in three-dimensional Euclidean space andis represented by an implicit equation [4]. Generally speaking,porous sca ff old design methods based on TPMS can be classi-fied into two categories. In the first class of methods, a volumemesh model is embedded in an ambient TPMS, and the inter-section of them is taken as the porous sca ff old [5–8]. In thesecond class of methods, a regular TPMS unit is transformed ∗ Corresponding author: phone number: 86-571-87951860-8304, fax num-ber: 86-571-88206681, email: [email protected] into each hexahedron of a hexahedron mesh model, thus gen-erating a porous sca ff old [9–11]. However, the first class ofmethods can generate incomplete TPMS units near the bound-ary of a volume mesh model, and the second class of methodsmay cause discontinuities between two adjacent TPMS units.Moreover, in both method classes, porosity and pore size aredi ffi cult to control.In this study, we developed a method for generating heteroge-nous porous sca ff olds in a trivariate B-spline solid (TBSS) withTPMS designed in the parametric domain of the TBSS. We alsodeveloped a porous sca ff old storage format that saves signifi-cant storage space. Specifically, given a TBSS, a threshold dis-tribution field (TDF) is first constructed in the cubic parameterdomain of the TBSS. Based on the TDF, a TPMS is generatedin the parameter domain. Finally, by mapping the TPMS in theparameter domain to the TBSS, a porous sca ff old is producedin the TBSS. In addition, the TDF can be modified locally inthe parameter domain to improve the engineering performanceof the porous sca ff old. All of the TPMS units generated in theporous sca ff old is complete, and adjacent TPMS units are con-tinuously stitched. To summarize, the main contributions of thisstudy are as follows: • A TBSS is employed to generate porous sca ff olds, whichensures completeness of TPMS units, and continuity be-tween adjacent TPMS units. • Porosity is easy to control using the TDF defined in theparametric domain of the TBSS.
Preprint submitted to Elsevier August 7, 2019 a r X i v : . [ c s . G R ] A ug A storage format for porous sca ff olds is designed basedon the TDF in the parametric domain, which saves signif-icantly storage space.The remainder of this paper is organized as follows. In Sec-tion 1.1, we review related work on porous sca ff old design andTBSS generation methods. In Section 2, preliminaries on TBSSand TPMS are introduced. Moreover, the heterogeneous poroussca ff old generation method using a TBSS and TDF in its pa-rameter domain is presented in detail in Section 3. In Section 4,some experimental examples are presented to demonstrate thee ff ectiveness and e ffi ciency of the developed method. Finally,Section 5 concludes the paper. In this section, we review some related work on porous scaf-fold design and TBSS generation methods.
Porous sca ff old design: In recent years, TPMS has been ofspecial interest to the porous sca ff old design community ow-ing to its excellent properties, and many sca ff old design meth-ods have been developed based on TPMS. Rajagopalan andRobb [12] made the first attempt to design tissue sca ff olds basedon Schwarz’s primitive minimal surface, which is a type ofTPMS. Moreover, the other two typical TPMSs (Schwarz’s di-amond surface and Schoen’s gyroid surface) are constructed byemploying K3DSurf software to design tissue sca ff olds [13],which achieve a gradient change of gyroid structure in terms ofpore size by adding a linear equation in z into the TPMS func-tion.Yoo [9] developed a method for generating porous sca ff oldsin a hexahedral mesh model, where the coordinate interpola-tion algorithm and shape function method are employed to mapTPMS units to hexahedron elements to generate tissue scaf-folds. To reduce the time consume in trimming and re-meshingprocess of Boolean operations, a tissue sca ff old design methodbased on a hybrid method of distance field and TPMS was fur-ther proposed in [5]. Moreover, to make the porosity easier tocontrol in the design a heterogeneous porous sca ff old, Yoo [6]introduced a method based on an implicit interpolation algo-rithm that uses the thin-plate radial basis function. Similar tothe method of Yoo [6], Yang et al. [7] introduced the sigmoidfunction and Gaussian radial basis function to design tissuesca ff olds. However, the hexahedral mesh based porous sca ff oldgeneration methods cannot ensure continuity between adjacentTPMS units.Recently, in consideration of the increasing attention towardsgradient porous sca ff olds, Shi et al. [11] utilized the TPMS andsigmoid function to generate functional gradient bionic poroussca ff olds from Micro-CT data reconstruction. Feng et al. [8]proposed a method to design porous sca ff old based on solid T-splines and TPMS, and analyzed the parameter influences onthe volume specific surface area and porosity. In addition, aheterogenous methodology for modeling porous sca ff olds usinga parameterized hexahedral mesh and TPMS was developed byChen et al. [10]. TBSS generation:
TBSS modeling methods are developedmainly for producing three dimensional physical domain in iso- geometric analysis [14]. Specifically, to analyze arterial bloodflow through isogeometric analysis, Zhang et al. [15] intro-duced a skeleton-based method of generating trivariate non-uniform rational basis spline (NURBS) solids. In [16], a tetra-hedral mesh model is parameterized through discrete volu-metric harmonic functions and a cylinder-like TBSS is gen-erated. Aigner et al. [17] proposed a variational frameworkfor generating NURBS parameterizations of swept volumes us-ing the given boundary conditions and guiding curves. Opti-mization approaches have been developed for filling boundary-represented models to produce TBSSs with positive Jacobianvalues [18]. Moreover, a discrete volume parameterizationmethod for tetrahedral mesh models and an iterative fitting al-gorithm have been presented for TBSS generation [19].
2. Preliminaries
Preliminaries on TBSS and TPMS are introduced in this sec-tion.
A B-spline curve of order p + p , and a B-spline curve has C p continuity at its breakpoints [20]. A knot vector U = { u , u , . . . , u m + p + } is defined by a set of breakpoints u ≤ u ≤· · · ≤ u m + p + . The associated B-spline basis functions N i , p ( u ) ofdegree p are defined as follows: N i , = , f or u i ≤ u < u i + , , otherwise , N i , p ( u ) = u − u i u i + p − u i N i , p − ( u ) + u i + p + − uu i + p + − u i + N i + , p − ( u ) (1)A TBSS of degree ( p , q , r ) is a tensor product volume definedas P ( u , v , w ) = m (cid:88) i = n (cid:88) j = l (cid:88) k = N i , p ( u ) N j , q ( v ) N k , r ( w ) P i jk (2)where P i jk , i = , , · · · , m , j = , , · · · , n , k = , , · · · , l arecontrol points in the u , v and w directions and N i , p ( u ) , N j , q ( v ) , N k , r ( w )are the B-spline basis functions of degree p in the u direction,degree q in the v direction, and degree r in the w direction.In this study, the input to our porous sca ff old generation algo-rithm is a TBSS that represents geometry at a macro-structuralscale. The TBSS can be generated either by fitting the meshvertices of a tetrahedral mesh model [19], or filling a closedtriangular mesh model [18]. A TPMS is an implicit surface that is infinite and periodic inthree independent directions. There are several ways to evalu-ate a TPMS, and the most frequently employed approach is to2pproximate the TPMS using a periodic nodal surface definedby a Fourier series [21], ψ ( r ) = (cid:88) k A k cos [2 π ( h k · r ) /λ k − P k ] = C , (3)where r is the location vector in Euclidean space, A k is ampli-tude, h k is the k th lattice vector in reciprocal space, λ k is thewavelength of the period, P k is phase shift, and C is a thresholdconstant. Please refer to [22] for more details on the abovemen-tioned parameters. The nodal approximations of P, D, G, andI-WP types of TPMSs, which were presented in [4], are listedin Table 1. The valid range of C guarantees that the implicitsurface is complete.In TPMS-based porous sca ff old design methods, the thresh-old value C (3) controls the porosity, and the coe ffi cients ω x , ω y ,and ω z (refer to Table 1) which a ff ect the period of the TPMS,are called period coe ffi cients . The e ff ects of the two types ofparameters in porous sca ff old design have been discussed in de-tail in the literature [8]. Additionally, in this study, the march-ing tetrahedra (MT) algorithm [23] is employed to extract theTPMS (shown in Fig. 1).
3. Methodology of porous sca ff old design The whole procedure of the heterogenous porous sca ff oldgeneration method based on TBSS and TPMS in its paramet-ric domain is illustrated in Fig.2. Specifically, given a TBSS asinput, we design a method for constructing the TDF in its cubicparameter domain. Based on the TDF, a TPMS is generated inthe parameter domain using MT algorithm [23]. Moreover, bymapping the TPMS in the parameter domain to the TBSS bythe TBSS function (2), a porous sca ff old with compeletenessand continuity is produced in the TBSS. If the porous sca ff oldin the TBSS does not meet the engineering requirements, theTDF can be locally modified in the parameter domain, and theporous sca ff old in the TBSS can be rebuilt. Finally, we developa storage format for the porous sca ff old based on the TDF in theparametric domain, which saves significant storage space. Thedetails of the porous sca ff old design method are elucidated inthe following sections. After the TDF is constructed in the parameter domain, andthe period coe ffi cients are assigned, the TPMS ψ = C in the pa-rameter domain can be defined such that it is polygonized into atriangular mesh. As a type of iso-surface, a TPMS can be poly-gonized by many algorithms such as iterative refinement [12],Delaunay triangulation [24], marching cubes [25] and MT [23].To avoid ambiguous connection problems [26] and simplify theintersection types, the MT algorithm is adopted to polygonizethe TPMS in the parametric domain of a TBSS. For this pur-pose, the cubic parameter domain is uniformly divided into agrid. In our implementation, to balance the accuracy and stor-age cost of the porous sca ff old, the parametric domain is di-vided into a 100 × ×
100 grid.Moreover, we define three types of volume TPMS structures (refer to Fig. 3): • pore structure represented by ψ ≥ C , • rod structure represented by ψ ≤ C , • sheet structure represented by C − (cid:15) ≤ ψ ≤ C .However, the triangular meshes of the three types of volumeTPMS structures, generated by the polygonization, are open onthe six boundary faces of the parameter domain, but they shouldbe closed to form a solid. Take the pore structure ( ψ ≥ C )as an example. In the polygonization procedure of the TPMS ψ = C by the MT algorithm, the triangles on the boundaryfaces of the parameter domain are categorized into two classesby the iso-value curve ψ = C on the boundary faces: outsidetriangles, where the values of ψ at the vertices of these trianglesare larger than or equal to C , and inside triangles, where thevalues of ψ at the vertices of these triangles are smaller thanor equal to C , Therefore, the pore structure can be closed byadding the outside triangles into the triangular mesh generatedby polygonizing ψ = C . The other two types of volume TPMSstructures can be closed in similar ways. In Fig. 3, the threetypes of volume TPMS structures are illustrated. The porosity is an important parameter in porous sca ff old de-sign because it has direct influences on the transport of nutritionand waste. The porosity and pore size of porous sca ff olds de-signed by TPMS units can be controlled by adjusting the thresh-old C (see Table 1). Moreover, the relationship between theporosity and the threshold C is illustrated in Figs.4-6. We cansee that, for the three types of volume TPMS structures, i.e.,pore, rod, and sheet, each has valid threshold range and theyare listed in Table 1.To design heterogeneous porous sca ff olds, we change thethreshold C to a trivariate function C ( u , v , w ) defined on theparametric domain of a TBSS, which is called threshold dis-tribution field (TDF). Then, the TPMS in the parameter domain is represented by the zero-level surface of, f ( u , v , w ) = ψ ( u , v , w ) − C ( u , v , w ) = . (4)Therefore, the TDF plays a critical role in the heterogeneousporous sca ff old generation, and designing the TDF becomes akey problem in porous sca ff old design.In this study, we developed some convenient techniques forgenerating a TDF in the cubic parameter domain of a TBSS. Inbrief, the cubic parameter domain of a TBSS is first discretizedinto a dense grid (in our implementation, it is discretized into agrid with a resolution of 50 × × parametric grid .Then, the threshold values at the grid vertices are assigned usingthe techniques presented later in this paper, which constitute a discrete TDF . Finally, to save storage space, the discrete TDFis fitted by a trivariate B-spline function, i.e., C ( u , v , w ) = n u (cid:88) i = n v (cid:88) j = n w (cid:88) k = N i , p ( u ) N j , q ( v ) N k , r ( w ) C i jk , (5)where the scales C i jk are the control points of the trivariate B-spline function.3 able 1: Nodal approximations of typical TPMS units. TPMS Nodal approximations Valid range of C Schwarz’s P Surface ψ P ( x , y , z ) = cos ( ω x x ) + cos ( ω y y ) + cos ( ω z z ) = C [ − . , . ψ D ( x , y , z ) = cos ( ω x x ) cos ( ω y y ) cos ( ω z z ) − sin ( ω x x ) sin ( ω y y ) sin ( ω z z ) = C [ − . , . ψ G ( x , y , z ) = sin ( ω x x ) cos ( ω y y ) + sin ( ω y y ) cos ( ω z z ) + sin ( ω z z ) cos ( ω x x ) = C [ − . , . ψ I − WP ( x , y , z ) = cos ( ω x x ) cos ( ω y y ) + cos ( ω y y ) cos ( ω z z ) + cos ( ω z z ) cos ( ω x x )] − [ cos (2 ω x x ) + cos (2 ω y y ) + cos (2 ω z z )] = C [ − . , . (a) (b) (c) (d)Figure 1: Four types of TPMS units. (a) P-type. (b) D-type. (c) G-type. (d) I-WP-type.
Figure 2:
The procedure of the heterogenous porous sca ff old generationmethod. Now, we elucidate the techniques for generating the discreteTDF.
Filling method.
Initially, all of the values at the parametergrid vertices are set to 0. Then, the geometric quantities atpoints of the boundary surface of TBSS, such as mean curva-ture, Gauss curvature, are calculated and mapped to the bound-ary vertices of the parametric grid. Furthermore, the quanti-ties at the boundary vertices are di ff used into the inner para- metric grid vertices by the Laplace smoothing operation [27].Thus, the entire parametric grid is filled, and a discrete TDFis constructed. In Fig. 7, the mean curvature distribution ofthe TBSS boundary surface is first calculated (Fig. 7(a)) andthen is mapped to the boundary vertices of the parametric grid(Fig. 7(b)). Layer method.
The parametric grid vertices are classified intolayers according to their coordinates ( u i , v j , w k ). For example,the vertices with the same w k coordinates can be classified intothe same layer. Vertices of the same layer are assigned the samethreshold values. The TDFs in Figs. 8(a) and 9(e) were gener-ated using the layer method. In Fig. 8(a), the vertices of the fourside surfaces are taken as the first layer, the vertices adjacent tothe first layer are taken being in the second layer, · · · , and soon. In Fig. 9(e), the vertices with the same w k coordinates aretaken as the same layer. Prescribed function method.
The threshold values at theparametric grid vertices can be assigned by a function pre-scribed by users. For example, in Fig. 17(b), the thresholdvalue at the vertex with coordinates ( u i , v j , w k ) is assigned bythe function, f ( u i , v j , w k ) = (cid:12)(cid:12)(cid:12) u i − v j (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) v j − w k (cid:12)(cid:12)(cid:12) + | u i − w k | . After the discrete TDF is generated, the values at the grid ver-tices are linearly transformed into the valid threshold range ac-cording to the type of volume TPMS structure being produced.Then, we fit the discrete TDF with a trivariate B-spline func-tion (5), using the least squares progressive-iteration approx-imation (LSPIA) method [28]. The subscripts of the grid ver-tices, i.e., ( i , j , k ), are the natural parametrization of the vertices.For the purpose of B-spline fitting, they are normalized into theinterval [0 , × [0 , × [0 , a) (b) (c) (d)(e) (f) (g) (h)(i) (j) (k) (l)Figure 3: Three types of volume TPMS structures for the four types of TPMS units. (a)(b)(c)(d) Pore structures for the P-type, D-type, G-type andI-WP-type TPMSs. (e)(f)(g)(h) Rod structures for the P-type, D-type, G-type and I-WP-type TPMSs. (i)(j)(k)(l) Sheet structures for the P-type,D-type, G-type and I-WP-type TPMSs.
Figure 4:
Relationship between the threshold C and the porosity of thefour types of TPMSs based on pore structures. of the trivariate B-spline function is taken as 20 × ×
20. Addi-tionally, the initial values for LSPIA iteration at the control gridof the B-spline function (5) are produced by linear interpolationof the discrete TDF.
Figure 5:
Relationship between the threshold C and the porosity of thefour types of TPMSs based on rod structures. Suppose the LSPIA iteration has been performed for l steps,5 igure 6: Relationship between the threshold C and the porosity of thefour types of TPMSs based on sheet structures. (a) (b)Figure 7: Filling method. (a) Mean curvature distribution on the TBSSboundary surface. (b) Mean curvature distribution on the boundaryvertices of the parametric grid. and the l th B-spline function C ( l ) ( u , v , w ) is constructed: C ( l ) ( u , v , w ) = n u (cid:88) i = n v (cid:88) j = n w (cid:88) k = N i , p ( u ) N j , q ( v ) N k , r ( w ) C ( l ) i jk . (6)to generate the ( l + th B-spline function C ( l + ( u , v , w ), the dif-ference vector for each parametric grid vertex is calculated, δ ( l ) α,β,γ = T α,β,γ − C ( l ) ( u α , v β , w γ ) , (7)where T α,β,γ is the threshold value at the vertex ( α, β, γ ), and( u α , v β , w γ ) are its parameters. Each di ff erence vector δ ( l ) α,β,γ isdistributed to the control points C ( k ) i , j , k if the corresponding ba-sis functions N i , p ( u α ) N j , q ( v β ) N k , r ( w γ ) are non-zero. Moreover,a weighted average of all the di ff erence vectors distributed toa control point is taken, leading to the di ff erence vector for thecontrol point , ∆ ( l ) i jk = (cid:80) I ∈ I αβγ N i , p ( u I ) N j , q ( v I ) N k , r ( w I ) δ ( l ) l (cid:80) I ∈ I αβγ N i , p ( u I ) N j , q ( v I ) N k , r ( w I ) , (8)where I αβγ is the set of indices such that N i , p ( u α ) N j , q ( v β ) N k , r ( w γ ) (cid:44) . Next, the ( l + th control points C ( l + i jk are formed by addingthe di ff erence vectors ∆ ( l ) i jk to the l th control points, C ( l + i jk = C ( l ) i jk + ∆ ( l ) i jk . (9)Thus, the ( l + th B-spline function C ( l + ( u , v , w ) is produced: C ( l + ( u , v , w ) = n u (cid:88) i = n v (cid:88) j = n w (cid:88) k = N i , p ( u ) N j , q ( v ) N k , r ( w ) C ( l + i jk . (10)The convergence of LSPIA iteration has been proved in [28].After the LSPIA iterations stop, the iteration result C ( u , v , w ) istaken as the TDF in the parametric domain. Local modification.
With the TDF C ( u , v , w ) in the paramet-ric domain, the TPMS in the parametric domain (4) can be gen-erated. By mapping the TPMS into the TBSS, a porous sca ff oldis produced in the TBSS. However, if the generated porous scaf-fold does not satisfy the practical engineering requirements, theTDF C ( u , v , w ) can be locally modified in the parametric do-main, and then the porous sca ff old can be rebuilt to meet thepractical requirements.To locally modify the TDF, users first choose some verticesof the parameter grid, and change their threshold values to theirdesirable values. Then, a local LSPIA iteration is invoked to fitthe changed values at the chosen vertices. In the local LSPIA it-eration, the di ff erence vector δ (7) is calculated only at the cho-sen vertices, and just the control points to which the di ff erencevectors δ are distributed are adjusted. The other control pointswithout distributed di ff erence vectors remain unchanged. Bylocally modifying the TDF in Fig. 8(a) using the method pre-sented above, the TDF is changed as illustrated in Fig. 8(b). (a) (b)Figure 8: Local modification of TDF. (a) TDF generated by the layermethod. (b) TDF after local modification. ff old in TBSS Until now, there was only one unknown in the heterogeneousporous sca ff old generation: the period coe ffi cients ω x , ω y , ω z (Table 1). It is worth noting that the internal connectivity of thesca ff old is crucial to the transferring performance of the scaf-fold. For a TPMS unit, the smaller the unit volume, the worsethe internal connectivity of the micro-holes, and the larger the6nit volume, the better the internal connectivity of the micro-holes. The period coe ffi cients ω x , ω y , ω z (Table 1) can be em-ployed to adjust the number of TPMS units, as well as the sizeof TPMS units in the three parametric directions. The values ofthe period coe ffi cients ω x , ω y , ω z used in the examples in thispaper are listed in Table 2.After the TDF in the parametric domain and the period coef-ficients are both determined (refer to Table 1), the TPMS in theparametric domain (Eq. (4)), i.e., f ( u , v , w ) = ψ ( u , v , w ) − C ( u , v , w ) = , can be calculated. For examples, the TPMSs (sheet structure)calculated based on the TDFs in Figs. 9(e)- 9(h) are illustratedin Figs. 9(a)- 9(d). It can be seen clearly from Figs. 9(a)- 9(d)that, the porosity is controlled by the TDF in Figs. 9(e)- 9(h).Finally, the heterogeneous porous sca ff old (Fig. 10(b)) in theTBSS can be generated by mapping the TPMS in the parametricdomain (volume TPMS structures) (Fig. 10(a)) into the TBSS,using the TBSS function. It should be noted that, to avoid fold-up, the Jacobian value of the TBSS should be positive. Becausethe TPMS in the parametric domain is unitary (Fig. 10(a)), ithas completeness and continuity between adjacent TPMS units.Therefore, the heterogeneous porous sca ff old in the TBSS isensured to be complete and continuous (Fig.10(b)). Owing to their complicated geometric and topological struc-ture, the storage costs for porous sca ff olds are very large, usu-ally requiring hundreds of megabytes (MB) (refer to Table 2).Therefore, the large storage cost becomes the bottleneck inporous sca ff old generation and processing. In this study, we de-veloped a new porous sca ff old storage format that reduces thestorage cost of porous sca ff olds significantly. Using the newstorage format, the space required to store the porous sca ff oldmodels presented in this paper ranges from 0 .
567 MB to 1 . .
402 MB to 1449 .
71 MB. Thus, the newstorage format reduces the storage requirement by at least 98%compared with the traditional STL file format. Moreover, thegeneration of heterogenous porous sca ff olds from the new fileformat costs a few seconds to a dozen seconds (Table 2).Specifically, because the TDF in the parametric domain andthe period coe ffi cients (Table 1) entirely determine the het-erogenous porous sca ff old in a TBSS, the new storage formatmust only store the following information: • period coe ffi cients ω x , ω y , ω z , • control points of the TDF C ( u , v , w ), • knot vectors of the TDF C ( u , v , w ), • control points of the TBSS P ( u , v , w ), • knot vectors of the TBSS P ( u , v , w ).Therefore, the new storage format is called the TDF format , andis summarized in Appendix for clarity.
4. Implementation, results and discussions
The developed heterogeneous porous sca ff old generationmethod is implemented in the C ++ programming language andtested on a PC with a 3.60 GHz i7-4790 CPU and 16 GB RAM.In this section, some examples are presented, and some im-plementation details are discussed. Moreover, the developedmethod is compared with classical sca ff old generation methods. ff old The heterogeneity of a porous sca ff old, i.e., its pore size andshape, is controlled by the TDF C ( u , v , w ) (4). The larger thevalue of C ( u , v , w ), the larger the pore size. In Fig. 11(a), theTDF is generated by the layer method, which takes the para-metric grid vertices with the same w coordinates as the samelayer. The small to large TDF values are visualized by blueto red colors. Fig. 11(b) illustrates the heterogeneous poroussca ff old (P-type, pore structure) generated based on the TDF inFig. 11(a). We can see that, with the TDF values varying fromlarge to small (Fig. 11(a)), the pore size of the porous sca ff oldalso changes from large to small (Fig. 11(b)). ff old generationmethods The heterogenous porous sca ff old generation method devel-oped in this study is compared here with two classical methodspresented in Refs. [8, 9]. First, because the method proposedin [9] generates a porous sca ff old by mapping a regular TPMSunit to each hexahedron of a hexahedral mesh model, it has thefollowing shortcomings: (1) Continuity cannot be ensured be-tween two adjacent TPMS units in the porous sca ff old. (2) Thethreshold values C of all TPMS units are the same. (3) The ge-ometric quality of the porous sca ff old is greatly influenced bythe mesh quality of the hexahedral mesh model. As illustratedin Fig. 12(a), the boundary mesh quality is poor, with manyslender triangles.Secondly, the method presented in [8] produces a poroussca ff old by first immersing a trivariate T-spline model in anambient TPMS and then taking the intersection of them as theporous sca ff old. Therefore, this method cannot guarantee com-pleteness of the TPMS units. As demonstrated in Fig. 12(b),many boundary TPMS units are broken.However, because our method generates a heterogeneousporous sca ff old by mapping a unitary TPMS in the paramet-ric domain to a TBSS, it avoids the shortcomings of the othertwo methods. The heterogeneous porous sca ff old generated byour method has the following properties (Fig. 12(c)): • Completeness of TPMS units and continuity between ad-jacent TPMS units are guaranteed. • The TDF can be designed by users to easily control theporosity. • Because the degree of a TBSS is relatively high, and aTBSS has high smoothness, the heterogenous porous scaf-fold generated by TBSS mapping is highly smooth. • The TDF file format saves significant storage space.7 a) (b) (c) (d)(e) (f) (g) (h)Figure 9:
Generation of the TPMSs (with their three-view drawing) (a,b,c,d) based on the corresponding TDF in the parametric domain (e,f,g,h).Note that the porosity of the TPMS is controlled by the TDF.
Table 2: Statistical data of the heterogenous porous sca ff old generation method developed in this paper.Model Type Structure Period coe ffi cients Run time(s) Storage space(MB) TDF TPMS Porous sca ff old STL format TDF formatBall joint P pore (16 , ,
18) 2.745 0.326 3.194 741.367 0.810rod 0.319 3.109 721.599sheet 0.658 6.651 1449.71Venus D pore (10 , ,
10) 2.728 0.293 3.024 701.682 0.947rod 0.286 3.007 695.564sheet 0.528 6.581 1399.55Moai I-WP pore (6 , ,
16) 2.736 0.290 2.554 557.324 0.824rod 0.284 2.555 554.906sheet 0.583 5.530 1105.41Tooth G pore (8 , ,
8) 2.732 0.238 1.974 394.402 0.567rod 0.237 1.956 394.422sheet 0.499 4.156 773.741Isis P pore (6 , ,
16) 2.754 0.238 2.089 473.401 1.355rod 0.234 1.994 454.359sheet 0.521 4.236 908.205 Run time (in second) for TDF construction, generation of volume TPMS in parametric domain and generation of heteroge-nous porous sca ff old. Storage space (in megabyte) of heterogenous porous sca ff olds using the traditional STL file format and the TDF file formatdeveloped in this paper. In this section, some heterogenous porous sca ff old resultsgenerated by our method are presented (Figs. 13-17) to demon-strate the e ff ectiveness and e ffi ciency of the method. InFigs. 13-17, (a) is the input TBSS, (b) is the TDF in the para-metric domain, and (c,d, and e) are the heterogeneous poroussca ff olds of pore structure, rod structure, and sheet structurewith di ff erent TPMS types. Moreover, statistical data are listedin Table 2, including, period coe ffi cients ( ω x , ω y , ω z ), run timesfor generating the TDF, TPMS in parametric domain, and het-erogenous porous sca ff old in the TBSS, and storage costs ofporous sca ff olds with the traditional STL file format and newTDF file format. In Table 2, the storage spaces required to store the poroussca ff old using the traditional STL file format and the new TDFfile format are listed. Using the TDF file format, storing theporous sca ff olds costs 0 .
567 to 1 .
355 MB, while using the STLfile format, it costs 394 .
402 to 1449 .
71 MB. Therefore, at least98% of storage space is saved by using the new TDF file format.Moreover, in Table 2, the time cost for generating the heteroge-neous porous sca ff old from the TDF file format is also listed,including the run time for generating volume TPMS structuresand porous sca ff olds. We can see that the time costs range from2 to 7 seconds, which is acceptable for user interaction.Finally, the TDF file format brings some extra benefits. Tra-ditionally, heterogeneous porous sca ff olds have been stored as8 a) (b)Figure 10: Generation of heterogeneous porous sca ff old. (a) TPMSin the parametric domain. (b) Heterogeneous porous sca ff old in theTBSS. (a) (b)Figure 11: Influence of threshold on the heterogeneity of sca ff old. (a)TDF generated by the layer method. (b) Heterogenous porous sca ff old(P-type, pore structure), produced based on the TDF in (a). linear mesh models. However, the TDF file format stores atrivariate B-spline function. Therefore, in theory, a porous scaf-fold can be generated to any prescribed precision using theTDF file format. In addition, the period coe ffi cients and controlpoints of the trivariate B-spline function, stored in the TDF fileformat, can be taken as some types of parameters . Therefore, aheterogeneous porous sca ff old can be changed by altering theseparameters, just like in the parametric modeling technique.
5. Conclusion
In this study, we developed a method for generating a het-erogeneous porous sca ff old in a TBSS by the TDF designed inthe parametric domain of the TBSS. First, the TDF is easy tobe designed in the cubic parameter domain, and is representedas a trivariate B-spline function. The TDF can be employed tocontrol the porosity of the porous sca ff old. Then, a TPMS canbe generated in the parameter domain based on the TDF andthe period coe ffi cients. Finally, by mapping the TPMS into theTBSS, a heterogeneous porous sca ff old is produced. Moreover,we presented a new file format (TDF) for storing the poroussca ff old that saves significant storage space. By the methoddeveloped in this study, both completeness of the TPMS unitsand continuity between adjacent TPMS units can be guaran- (a)(b)(c)Figure 12: Comparison with classical porous sca ff old generating meth-ods. The left-most column is the models input to the correspondingmethods. (a)Method in [9] with hexahedral mesh as input. (b)Methodin [8] with T-spline solid as input. (c)Our method with TBSS as input. teed. Moreover, the porosity of the porous sca ff old can be con-trolled easily by designing a suitable TDF. More importantly,the TDF file format not only saves significant storage space,but it can also be used to generate a porous sca ff old to any pre-scribed precision. In terms of future work, determining howto change the shape of a porous sca ff old using the parametersstored in the TDF file format is a promising research direction. Acknowledgements
This work is supported by the National Natural ScienceFoundation of China under Grant No.61872316.
References [1] D. W. Hutmacher, Sca ff olds in tissue engineering bone and cartilage.,Biomaterials 21 (24) (2000) 2529–2543.[2] B. Starly, C. Gomez, A. Darling, Z. Fang, A. Lau, W. Sun, W. Lau,T. Bradbury, A. Youssef, C. Gaylo, Computer-aided bone sca ff old de-sign: a biomimetic approach, in: 2003 IEEE 29th Annual Proceedings ofBioengineering Conference, 2003, pp. 172–173. a) (b) (c) (d) (e)Figure 13: Heterogeneous porous sca ff old of Ball-joint . (a) TBSS. (b) TDF in the parametric domain. (c) P-type pore structure. (d) P-type rodstructure. (e) P-type sheet structure. (a) (b) (c) (d) (e)Figure 14:
Heterogeneous porous sca ff old of Venus . (a) TBSS. (b) TDF in the parametric domain. (c) D-type pore structure. (d) D-type rodstructure. (e) D-type sheet structure. (a) (b) (c) (d) (e)Figure 15:
Heterogeneous porous sca ff old of Moai . (a) TBSS. (b) TDF in the parametric domain. (c) I-WP-type pore structure. (d) I-WP-type rodstructure. (e) I-WP-type sheet structure. [3] A. Y¢nez, A. Cuadrado, O. Martel, H. Afonso, D. Monopoli, Gyroidporous titanium structures: A versatile solution to be used as sca ff oldsin bone defect reconstruction, Materials & Design 140 (2018) 21–29.[4] H. G. V. Schnering, R. Nesper, Nodal surfaces of fourier series: Funda-mental invariants of structured matter, Zeitschrift Fr Physik B Condensed Matter 83 (3) (1991) 407–412.[5] D. Yoo, Porous sca ff old design using the distance field and triply periodicminimal surface models, Biomaterials 32 (31) (2011) 7741–7754.[6] D. Yoo, Heterogeneous minimal surface porous sca ff old design using thedistance field and radial basis functions, Medical Engineering & Physics a) (b) (c) (d) (e)Figure 16: Heterogeneous porous sca ff old of Tooth . (a) TBSS. (b) TDF in the parametric domain. (c) G-type pore structure. (d) G-type rodstructure. (e) G-type sheet structure. (a) (b) (c) (d) (e)Figure 17:
Heterogeneous porous sca ff old of Isis . (a) TBSS. (b) TDF in the parametric domain. (c) P-type pore structure. (d) P-type rod structure.(e) P-type sheet structure.
34 (5) (2012) 625–639.[7] N. Yang, K. Zhou, E ff ective method for multi-scale gradient poroussca ff old design and fabrication., Materials Science & Engineering C 43(2014) 502–505.[8] J. Feng, J. Fu, C. Shang, Z. Lin, B. Li, Porous sca ff old design by solidt-splines and triply periodic minimal surfaces, Computer Methods in Ap-plied Mechanics & Engineering 336 (2018) 333–352.[9] D. Yoo, Computer-aided porous sca ff old design for tissue engineering us-ing triply periodic minimal surfaces, International Journal of PrecisionEngineering & Manufacturing 12 (1) (2011) 61–71.[10] H. Chen, Y. Guo, R. Rostami, S. Fan, K. Tang, Z. Yu, Porous structure de-sign using parameterized hexahedral meshes and triply periodic minimalsurfaces, in: Proceedings of Computer Graphics International 2018, CGI2018, Bintan Island, Indonesia, June 11-14, 2018, ACM Press, 2018, pp.117–128.[11] J. Shi, L. Zhu, L. Li, Z. Li, J. Yang, X. Wang, A tpms-based method formodeling porous sca ff olds for bionic bone tissue engineering, ScientificReports 8 (1) (2018) 7395.[12] S. Rajagopalan, R. A. Robb, Schwarz meets schwann: design and fabrica-tion of biomorphic tissue engineering sca ff olds., Medical Image Analysis10 (5) (2006) 693–712.[13] F. P. Melchels, K. Bertoldi, R. Gabbrielli, A. H. Velders, J. Feijen, D. W.Grijpma, Mathematically defined tissue engineering sca ff old architec-tures prepared by stereolithography., Biomaterials 31 (27) (2010) 6909– 6916.[14] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad,finite elements, nurbs, exact geometry and mesh refinement, ComputerMethods in Applied Mechanics & Engineering 194 (39) (2005) 4135–4195.[15] Y. Zhang, Y. Bazilevs, S. Goswami, C. L. Bajaj, T. J. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of bloodflow, Computer methods in applied mechanics and engineering 196 (29)(2007) 2943–2959.[16] T. Martin, E. Cohen, R. M. Kirby, Volumetric parameterization andtrivariate b-spline fitting using harmonic functions, Computer Aided Ge-ometric Design 26 (6) (2009) 648–664.[17] M. Aigner, C. Heinrich, B. Jttler, E. Pilgerstorfer, B. Simeon, A. V.Vuong, Swept volume parameterization for isogeometric analysis, in:Mathematics of Surfaces XIII, Springer Berlin Heidelberg, Berlin, Hei-delberg, 2009, pp. 19–44.[18] X. Wang, X. Qian, An optimization approach for constructing trivariateb-spline solids, Computer-Aided Design 46 (1) (2014) 179–191.[19] H. Lin, S. Jin, Q. Hu, Z. Liu, Constructing b-spline solids from tetrahedralmeshes for isogeometric analysis, Computer Aided Geometric Design 35-36 (2015) 109–120.[20] A. Piegl, Les, W. Tiller, The NURBS Book, Springer Berlin Heidelberg,Berlin, Heidelberg, 1997.[21] P. J. F. Gandy, S. Bardhan, A. L. Mackay, J. Klinowski, Nodal surface pproximations to the ja:math and i-wp triply periodic minimal surfaces,Chemical Physics Letters 336 (3) (2001) 187–195.[22] W. Yan, Periodic surface modeling for computer aided nano design,Computer-Aided Design 39 (3) (2007) 179–189.[23] A. Doi, A. Koide, An e ffi cient method of triangulating equivalued sur-faces by using tetrahedral cells, Ieice Trans 74 (1) (1991) 214–224.[24] P. L. George, H. Borouchaki, Delaunay Triangulation and Meshing, Paris:´Edition Herm`es, 1998.[25] W. E. Lorensen, H. E. Cline, Marching cubes: A high resolution 3d sur-face construction algorithm, Acm Siggraph Computer Graphics 21 (4)(1987) 163–169.[26] T. S. Newman, H. Yi, A survey of the marching cubes algorithm, Com-puters & Graphics 30 (5) (2006) 854–879.[27] D. A. Field, Laplacian smoothing and delaunay triangulations, Commu-nications in Applied Numerical Methods 4 (6) (1988) 709–712.[28] C. Deng, H. Lin, Progressive and iterative approximation for least squaresb-spline curve and surface fitting, Computer-Aided Design 47 (1) (2014)32–44. Appendix: TDF file format ffi cients( ω x , ω y , ω z ) ω x ω y ω z n u + n v + n w + C , , C , , ... C n u , n v , n w u -direction of TDF u u · · · u n u + v -direction of TDF v v · · · v n v + w -direction of TDF w w · · · w n w + m + n + l + Px , , Py , , Pz , , Px , , Py , , Pz , , ... ... ... Px m , n , l Py m , n , l Pz m , n , l u -direction of TBSS U U · · · U m + v -direction of TBSS V V · · · V n + w -direction of TBSS W W · · · W l +4