FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 16, N o 3, 2018, pp. 285 - 296 https://doi.org/10.22190/FUME170618018R © 2018 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper A CAD-BASED CONCEPTUAL METHOD FOR SKULL PROSTHESIS MODELING UDC 004.922.8:616-089.843 Marcelo Rudek 1 , Yohan B. Gumiel 1 , Osiris Canciglieri Jr 1 , Naomi Asofu 1 , Gerson L. Bichinho 2 1 Production and System Engineering Graduate Program – PPGEPS, Brazil 2 Graduate Program in Health Technology – PPGTS, Pontifical Catholic University of Parana – PUCPR, Brazil Abstract. The geometric modeling of a personalized part of the tissue built according to individual morphology is an essential requirement in anatomic prosthesis. A 3D model to fill the missing areas in the skull bone requires a set of information sometimes unavailable. The unknown information can be estimated through a set of rules referenced to a similar yet known set of parameters of the similar CT image. The proposed method is based on the Cubic Bezier Curves descriptors generated by the de Casteljou algorithm in order to generate a control polygon. This control polygon can be compared to a similar CT slice in an image database. The level of similarity is evaluated by a meta-heuristic fitness function. The research shows that it is possible to reduce the amount of points in the analysis from the original edge to an equivalent Bezier curve defined by a minimum set of descriptors. A study case shows the feasibility of method through the interoperability between the prosthesis descriptors and the CAD environment. Key Words: Prosthesis Modeling, CAD, Bezier Curves, 3D Image 1. INTRODUCTION The aging of the world population and thereby caused increasing demands for medical services set new challenges in the field of biomedical engineering. This is especially true in the field of image processing, design of personalized prosthesis and automated manufacturing. In the context of machining process, the additive technologies are capable of building complex structures in different materials geometrically compatible with the Received June 18, 2017 / Accepted November 21, 2017 Corresponding author: Marcelo Rudek Pontifical Catholic University of Parana, PUCPR/ Production and System Engineering Graduate Program – PPGEPS, Imaculada Conceicao, 1155, 80215-030, Curitiba, Brazil. E-mail: marcelo.rudek@pucpr.br 286 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR, N. ASOFU, G. L. BICHINHO tissues of human body. The reason why the additive technologies are so much applied in the personalized prosthesis production is the fact that it is possible to build a part of arbitrary complexity once you have its 3D geometric model. The congenital failure or trauma in the skull bone requires surgical procedures for prosthesis implant as functional or esthetical repairing. In this process, a personalized prosthesis built according to individual morphology is an essential requirement. Normally in bone repairing, the geometric structure is unrepeatable due to its “free form” [1]. Due to the complexity in geometry, free form objects do not have a mathematical expression in a close form to define their structure. However, numerical approximations are a feasible way to the geometric representation. The link between the medical problem and the respective manufactured product (i.e. prosthesis) is the geometric modeling; thus different approaches in the bone modeling have opened new research interests as in [2, 3, 4]. In the prosthesis modeling, we face different levels of information handling from a low level of the pixel analysis in image to the automated production procedure. In general, there are following levels: a “preparation level” (containing CT scanning, segmentation, feature extraction, i.e. entire image processing) and a “geometric modeling level” (containing the polygonal model, curve model, extraction of anatomic features, i.e. entire CAD based operations) [2]. The CAD systems are important tools in the design of these complex products because they can be used for three-dimensional (3D) modeling of the shapes of bones and respective scaffolds to machining [5]. In our strategy, we need to generate geometric representation of the bone without enough information (e.g. neither mirroring nor symmetry applicable). A common image segmentation procedure is executed as pre-processing at the preparation level. Moreover, from the segmented images we define a set of descriptors based on Bezier Curves [6] in order to describe the skull edge geometry on a CT image. This approach is applied in order to reduce the amount of points capable of representing the skull bone curvature. It was adapted from the method of [4] now using the de Casteljau algorithm to define the Bezier parameters. The paper explores the accuracy of the prosthesis modeled through the balanced relationship between curve fitness versus number of descriptors. It extends the work of [3] to demonstrate the relationship of descriptors within a CAD system in order to build a 3D prosthesis piece. 2. THE PROPOSED METHOD 2.1. The conceptual model In our study, the main question is related to the information recovering for automation of the prosthesis modeling process. Sometimes it is possible to reconstruct a fragmented image by using information of the same bone structure, e. g. by mirroring, that is, using body symmetry from the same individual. However, in many cases, there is not enough information to be mirrored. A handmade procedure can be performed by a specialized doctor by using a CAD system [7, 8, 9]. In order to circumvent mirroring limitations and the user‟s intervention, we are looking for an autonomous process of geometric modeling of skull prosthesis. Thus, the basis of our hypothesis is to find compatible information from different healthy individuals from image database. The problem addressed here is the method for finding a compatible intact CT slice to replace the respective defective CT slice. When working with medical images [10], a lot of A CAD-Based Conceptual Method for Skull Prosthesis Modelling 287 information is needed to be handled mainly after image segmentation and edge detection, where the total of pixels in edge are still too much information to be processed. Our approach is a content-based retrieval procedure contrary to the pixel-by-pixel comparison that is a hard processing task that we need to avoid. In order to optimize the search by similarity, we propose to define shape descriptors by Cubic Bezier Curves. In this way it is possible to reduce the amount of data-to-process to a few parameters. The next important issue is to find the descriptors capable of describing the edge shape as well as possible. Thus, we also look for a balance between accuracy and the minimum quantity of information. The next section will explain our approach in curve modeling. 2.2. The curvature representation The curve modelling adopted in this research is based on the de Casteljau algorithm applied in calculation of the points of a Bézier Curve [3, 6]. The de Casteljou method [11, 12] operates by the definition of a “control polygon” whose vertices are respective “control points” (or “anchor points”) used to define the shape of the Bezier Curve. A Bezier curve of degree n is built by n+1 control points. The Cubic Bezier Curve has two endpoints (fixed points) and two variable points. They define the shape (flatness) of curve. Figure 1 shows an example of a Cubic Bezier Curve where {P0, P1, P2, P3} are the vertices of the control polygon. Points {P0, P3} are fixed and they are the beginning and the ending of the curve, respectively; - these points belong to the curve. {P1, P2} are variable points occupying any random position in  2 . Fig. 1 Graphical Representation of de Casteljau method [11] According to Fig. 1, for all points i r i PtP )( , we have a )( 0 tP n as a point on the Bezier curve. Bezier curve B n (t) with degree 'n' is a set of points )(0 tP n , t[0, 1] , i.e. ]}1,0[);({)( 0  ttPtB nn . Then the polygon formed by n vertices {P0, P1,…, Pn} is so called “control polygon” (or Bezier polygon) [12]. Through the de Casteljau algorithm each line segment results in (n1) baselines as 10 PP , 21PP , 32 PP which are recursively divided to define a new set of control points. By changing 't' value as defined in Eq. (2) we obtain the position of the point in the curve:          rni nr tPttPttP r i r i r i ,...,2,1 ,...,2,1 )()()1()( 1 1 1 (1) 288 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR, N. ASOFU, G. L. BICHINHO 02 01 tt tt t    (2) The control points for P[t0 t1](t), are n PPPP 0 2 0 1 0 0 0 ,...,,, , and the control points for P[t1 t2](t) are 02 2 1 10 ,...,,, n nnn PPPP  . In order to avoid misunderstanding in representation, Fig. 2 shows the control points, and the recursive subdivision of the de Casteljau algorithm labelled as P, Q, R and S, where S is the final position of a point in the curve for different values of t. In Fig. 2a the value of t = 0.360 and in Fig. 2b the value of t = 0.770. a) b) Fig. 2 Position of the control points and its respective Bezier Curve adapted from [11]1 As presented in literature, for practical applications, the most common one is to apply the Cubic Bezier Curves (n=3) due to the large possibilities in adapting the shape (flatness) according to our necessities. Also in our proposal, the Bezier with n=3 is more suitable to fit the skull contour in tomographic cuts. An example of a segmented CT slice is shown in Fig. 3. a) b) c) d) Fig. 3 A CT slice sample and respective Bezier representation: a) A quadratic Bezier Curve (n=2); b) a quadratic Bezier Curve on small region; c) a Cubic Bezier Curve (n=3); d) The Cubic Bezier Curve on small region In Fig. 3a a Quadratic Bezier Curve (B n (t) with n = 2) adjusted on the skull edge is presented. In this case we have three control points and only two baselines. Note that the adjustment in the outer edge seems satisfactory but in the inner edge the result is poor. In 1 under an Attribution-NonCommercial-ShareAlike CC license (CC BY-NC-SA 3.0) A CAD-Based Conceptual Method for Skull Prosthesis Modelling 289 the same way, Fig. 3b shows the de Casteljau algorithm applied to the smallest segment of the inner edge; in this case, then, the curve representation is improved. In Fig. 3c a Cubic Bezier Curve is presented. Now more control points exist and the resulting adjustment looks very good for both the outer and the inner edge. Also, in Fig. 3d the method applied in a small section (the inner edge) is more accurate. The question that we intend to discuss in the next section concerns the similarity measurement. In other words, how good is the quality of a Bezier curve that represents a CT skull edge? This is essential for our approach because we need to define the best curve based descriptor. Good descriptors will permit us to retrieve compatible CT images to produce the skull prosthesis. 3. APPLICATION OF THE METHOD The aim of this research is to define a small set of descriptors to represent the bone curvature. The strategy is to use the Cubic Bezier Curve method calculated through the de Casteljou algorithm. In our previous section, we state that the accuracy of curve fitting in our approach by the Bézier depends on its degree n value and the length of the edge section. The edge sectioning is defined as follows: 3.1. The sectioning of the edge As presented above in Fig. 3, the curve generated on the edge seems to fit better to the smallest length region (i.e., the shape of the curve looks similar to the original edge shape). The first question is about the best number of sections to produce the best-fitted curve. As an example, the edges of a CT image can be sectioned as in Fig. 4. a) b) Fig. 4 A CT slice sample with: a) k=10 sections; b) k=20 sections Figure 4a shows the total of k=10 cuts (with P0 to P10 fixed points) whose section edges lengths are bigger than sections of Fig. 4b with k=20 cuts (with P0 to P20 fixed points). For each section, fitness value F is calculated using Eq. (3), defined in [3] as: 290 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR, N. ASOFU, G. L. BICHINHO    n i BB iyiyixixkF 1 2 0 2 0 ))()(())()(()( (3) where F(k) is the fitness value for each section k. Fitness F calculates the error between the Bezier coordinates (xB, yB) and the original edge coordinates (x0, y0) for each pixel „i‟ in the edge. The sectioning procedure and the control points calculation are fully covered in [3]. Table 1 shows the average of fitness (error) to respective 5, 10, 15 and 20 sections cuts. Table 1 Relationship between number of cuts and respective fitness (error) # of Section F(5) F(10) F(15) F(20) 1 61.81640 26.03820 8.96430 6.22810 2 87.16330 22.63830 17.94760 9.27160 3 99.14170 21.17040 20.07230 9.34800 4 78.18280 25.94340 16.91680 10.25800 5 68.10270 26.46990 18.00790 10.09840 6 - 30.86040 17.19840 9.31330 7 - 28.55470 19.37670 14.06260 8 - 24.36930 22.17160 9.27060 9 - 36.52540 17.10820 9.50420 10 - 48.82190 19.46810 11.26640 11 - - 19.99220 12.80600 12 - - 21.37690 7.48890 13 - - 18.15550 10.65490 14 - - 27.96400 12.36890 15 - - 40.85040 8.37890 16 - - 9.40440 17 - - - 9.74180 18 - - - 17.47880 19 - - - 15.90830 20 - - - 25.4706 Σ (error) 394.40690 291.39190 305.57090 228.32270 Fitness (avg.) 78.88138 29.13919 20.37139333 11.416135 Table 1 shows the cumulative error evaluated by Eq. (3) and the average of fitness for different values of sectioning. As expected, the error is minimized with larger values of k. The graph in Fig. 5 presents the relationship between the number of sections and the calculated error (difference between the original edge and the calculated Bézier). As presented in Fig. 5, the average of error calculated from the fitness equation goes down as the number of sections is increased. Then, in this condition, maybe we could define the k value as the maximum possible, i.e. the length of total of pixels of edge. However, the computational cost of the Cubic Bézier Curve calculation for hundreds of sections also increases. The same proportion of error occurs for all CT slices from different images. From the graph, selecting the value of k=20 is enough to match a relatively good fitness with a small error and give us an adequate balance between precision and computational cost. Thus, for k=20 we have in the de Casteljau algorithm, 20 “fixed points” and another 20 “variable points”, (i.e. 4 points per section) calculated as in [3]. Now, it is possible to represent the total length of each edge (inner and outer) in a A CAD-Based Conceptual Method for Skull Prosthesis Modelling 291 CT slice with 80 points descriptors each instead of ≈ 1250 in the original edge (around 15 times information reduced). 5 10 15 20 25 30 35 40 0 20 40 60 80 Sections X Error Number of Cut Sections F it n e s s V a lu e ( a v e ra g e ) Fig. 5 The relationship between the number of cuts in the edge and respective error from the fitted Bézier curve in each section 3.2. The curve fitting procedure The curve fitting procedure is applied to each CT slice of defective skull. The same procedure is also applied to each searching image on database. The compatible answer image is retrieval as in the example presented in Fig. 6. Fig. 6 The defective set of slices and respective compatible CT recovered from medical image database In Fig. 6 some samples are shown of the defective CT from the original dataset and respective retrieved CT with compatible descriptors (i. e., minimum error in descriptors). Fig. 6 also shows error value (E) for each images pairs. The error is the cumulative difference between the original bone and the calculated Bezier curve by applying Eq. (3). 4. EVALUATION OF THE METHOD A handmade testing failure is built-in in a skull through the FIJI software [13]. It is an open source Java suitable for a medical image analysis. A set of toolboxes permits us to handle CT slices from the DICOM file [10]. The edge from individual slices can be cut in 292 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR, N. ASOFU, G. L. BICHINHO sequence in order to build a failure in a region. Thus, after 3D reconstruction, we obtain a synthetically built failure in the skull as in the example in Fig. 7a. a) b) Fig. 7 Testing image: a) a handmade testing failure built on original image; b) region filled with prosthesis modeled Fig. 7b shows a failure region filled with compatible CT slices from the medical database. The piece is cut from different slices where Bezier descriptors are compatible, i.e. all those CTs retrieved with minimum error. The retrieved slices numbers and respective patient are shown in Table 2. Table 2 Retrieved set of CT slices Original Slice # 280 282 284 286 288 290 292 Compatible Individual # 7 6 6 6 5 7 6 Compatible CT image # 278 285 282 285 284 286 290 Fitness Value 130.2566 143.4175 128.9885 139.1415 189.7201 226.1416 140.9704 Note that the slices retrieved are never from the same patient. The set of retrieved slices (good slices) are recovered from healthy individuals (intact skull) whose descriptors match with the original image in each CT slice (defective slice). From Table 2 it is possible to see that many CTs are coming from the individual #6. In fact, the individual #6 has similar morphological characteristics with the testing patient of the same gender and of similar age. The filled region is evaluated through the Geomagic® software [14]. The differences between the original bone and the prosthesis piece are presented in Fig. 8. The software permits overlapping of 3D structures and provides a colored scale to show the spatial difference whose lower error values are represented in green while the higher error tends to red. As shown in Fig. 8, the region A008 has an error closest to zero because it is the original skull (skull of patient). The other colored identifications are in the prosthesis area like A001, A003, A005 and A006; they are below 1mm difference and the maximum error is in the region A004 with the value of about 1.7087 mm. A CAD-Based Conceptual Method for Skull Prosthesis Modelling 293 Fig. 8 3D evaluation of filled region 5. THE GEOMETRIC MODELING IN THE CAD SYSTEM A CAD system can be used for modeling of the bone shape as well as for generating prosthesis profiles used in machining preparation. The main objective of the presented approach is the verification of each CT slice by slice in order to create curvature descriptors; then the slices of this new 3D image are the profiles exported to a CAD for the creation of the virtual (missing) bone piece. Due to the geometrical complexity of the individual human bones, the model cannot be automatically solved in CAD systems, but a preliminary step is needed to bring guaranties that the geometry stays totally closed in order to apply CAD functions. Only after this preliminary step, the superposing of the created surfaces can be transformed into a solid geometry [15, 16, 17]. This model is what a medical team could examine in order to check whether it is suitable to be implanted in the patient and thus be translated into a CAM system for manufacturing. Assuming that the solid geometry can be generated by the CAD system, based on the tomography image, the dimensions of the raw material needed for the CAM process (machining) can be determined. A SolidWorks® [16] example is presented in order to demonstrate the method‟s feasibility. A cut region is handmade and built on a testing skull image. After the method of descriptors is applied, a similar set of CTs are retrieved from image database; the cloud of points is imported to SolidWorks and plotted in 3D as presented in Fig. 9a. Despite the fact that the cloud of points is enough to set machining procedures, a visualization resource by surfaces can improve a visual evaluation of the piece. Then, in the next step, we apply for each 2D plane its respective mesh and surface as presented respectively in the two samples in Figs. 9b and 9c. They are the top and bottom limits of the prosthesis surface. Also, another viewing possibility is the shape contour as a convex hull in each area generated by spline curves in Fig. 9d. 294 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR, N. ASOFU, G. L. BICHINHO a) b) c) d) Fig. 9 The imported calculated data to CAD system: a) cloud of points; b) mesh for each CT layer; c) surface based on imported data; d) splines applied in edge limit in each CT slice The spline curve around each 2D slice contour is the reference to the “loft” instruction in CAD. The Loft is a tool to create a 3D solid from cross sections, i.e., in our case it generates the 3D visualization based on the superimposed Splines from all CT surfaces. The resulting image of this process is presented in Figs. 10a and 10b. The resulting final image of prosthesis piece is presented in Figs. 10c and 10d. a) b) c) d) Fig. 10 The views of 3D modeled piece after loft in the CAD system A CAD-Based Conceptual Method for Skull Prosthesis Modelling 295 6. CONCLUSIONS The paper presents a method for generating skull shape descriptors based on the Bezier Curves whose parameters are generated by the de Casteljou algorithm. Edge sectioning in k=20 sections with the same length permits us to define two markers as respective “fixed points” in the Bezier curve generator. Two more “variable points” calculated by the de Casteljau define the total of 4 descriptors for each section. Thus, it is possible to reduce all edge size in CT to be represented by a set of 80 descriptors. The descriptors are used to look for compatible CT images whose bone edge shape versus Bezier curve calculated by their descriptors have a minimum error. The example shows the result with maximum error in image around 1.7mm. We show it is possible to represent a missing region of the patient‟s skull by a set of similar CTs from healthy individuals selected by a reduced descriptors group. All those descriptors can be exported to a CAD system to build a prosthesis piece. An example within the cloud of points, mesh and 3D surface is presented to illustrate the integration between data and its respective geometric model. In addition, the example shows that the retrieved slices are from individuals with similar characteristic as to age and gender. In a future work, the database searching engine can group individuals with these characteristics before proceeding to calculating descriptors. Acknowledgements: The authors would like to thank the Production and System Engineering Graduate Program – PPGEPS and the Graduate Program in Health Technology – PPGTS from Pontifical Catholic University of Parana – PUCPR by collaboration in providing all CT image data. Also, the authors would like to thank to the CNPQ Brazilian Grant Program for providing the financial support to young researchers. REFERENCES 1. Trifunovic, M., Stojkovic, M., Trajanovic, M., Manic, M., Misic, D., Vitkovic, N., 2015, Analysis of semantic features in free-form object reconstruction, Artificial Intelligence for Engineering Design, Analysis and manufacturing, 30(1), pp.1-20. 2. Majstorovic, V., Trajanovic, M., Vitkovic, N., Stojkovic, M., 2013, Reverse engineering of human bones by using method of anatomical features, CIRP Annals – Manufacturing Technology, 62, pp.167-170. 3. Rudek, M., Gumiel, Y. B., Canciglieri Jr.O., Bichinho, G. L., 2016, Optimized CT Skull Slices Retrieval based on Cubic Bezier Curves Descriptors”, Proceedings of ICIST 2016, Kopaonik, pp. 1-5. 4. Rudek, M., Gumiel, Y.B., Canciglieri Jr.O., 2015, Autonomous CT replacement method for the skull prosthesis modelling, Facta Universitatis-Series Mechanical Engineering, 13(3), pp. 283-294. 5. Liulan, L., Qingxi, H., Xianxu, H., Gaochun, X., 2007, Design and fabrication of bone tissue engineering scaffolds via rapid prototyping and CAD, Journal of Rare Earths, 25, pp. 379-383. 6. Shao, L.J., Zhow, H., 1996, Curve Fitting with Bezier Cubics, Graphical Models and Image Processing, 58(3), pp.223-232. 7. ***, 2016, Materialize Mimics, available at http://www.materialise.com/en/medical/software/mimics (Last access: 15.5.2016) 8. Osirix, Osirix Imaging Software – Advanced Open-Source Pacs Workstation DICOM Viewer, available at: http://www.osirix-viewer.com/ (Last access: 15.1.2016) 9. Nasr, E.A., Al-Ahmari, A., Moiduddim, K., Al Kindi, M., Kamrani, A., 2015, A Study on The Evaluation and Accuracy of Anatomic and Mirror Image Reconstruction Design Technique, Proceedings of 45 th Computer in Industry - CIE45, Metz, pp.1-8. 10. Dicom, Digital Imaging and Communications in Medicine Part 5: Data Structures and Encoding, available at: www.medical.nema.org/dicom. (Last access: 5.3.2015) 11. Christersson, M., 2014, De Casteljau's Algorithm and Bézier Curves, available at: http://www.malinc.se/ m/DeCasteljauAndBezier.php (Last access: 12.6.2015) 296 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR, N. ASOFU, G. L. BICHINHO 12. Simoni, R., 2005, Teoria Local das Curvas, undergraduation monograph, UFSC University, 2005, (in portuguese), 96p. 13. Schindelin, J., Carreras, I. A., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.-Y., White, D. J., Hartenstein, V., Eliceiri, K., Tomancak, P., Cardona, A., 2013, Fiji: an open-source platform for biological-image analysis, Nature Methods, 9(7), pp. 676-682. 14. Geomagic, 2015, Modeling Easter Island’s Moai with Geomagic 3D Scan Software, available in http://www.geomagic.com/en/ (Last access: 10.12.2016) 15. Greboge, T., Rudek, M., Canciglieri, Jr.O., 2011, Conceptual Geometric Model For Prosthesis Modeling In CAD System: A Case Study To Skull Repairing With Asymmetric Defect, 21st International Conference on Production Research - ICPR 21, Stuttgart, pp. 1-6. 16. SolidWorks, http://www.solidworks.com/ (Last access: 10.4.2015) 17. de Troyer, O., Billie, W., Kleinermann, F., 2009, Defining the Semantics of Conceptual Modeling Concepts for 3D Complex Objects in Virtual Reality, Journal on Data Semantics XIV, pp.1-33.