Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 13, N o 3, 2015, pp. 283 - 294 AUTONOMOUS CT REPLACEMENT METHOD FOR THE SKULL PROSTHESIS MODELLING  UDC 519.6:617.3 Marcelo Rudek, Yohan Boneski Gumiel, Osiris Canciglieri Jr Production and System Engineering Graduate Program – PPGEPS, Pontifical Catholic University of Parana, PUCPR, Brazil Abstract. The geometric modeling of prosthesis is a complex task from medical and engineering viewpoint. A method based on CT replacement is proposed in order to circumvent the related problems with the missing information to modeling. The method is based on digital image processing and swarm intelligence algorithm. In this approach, a missing region on the defective skull is represented by curvature descriptors. The main function of the descriptors is to simplify the skull’s contour geometry; and they are defined from the Cubic Bezier Curves using a meta-heuristic process for parameter’s estimation. The Artificial Bee Colony (ABC) optimization technique is applied in order to evaluate the best solution. The descriptors from a defective CT slice image are the searching parameters in medical image databases, and a similar image, i.e. with similar descriptors, can be retrieval and used to replace the defective slice. Thus, a prosthesis piece is automatically modeled with information extracted from distinct skulls with similar anatomical characteristics. Key Words: Prosthesis Modeling, Computed Tomography, Cubic Bezier Curves, Content Retrieval, ABC Optimization 1. INTRODUCTION Skull prosthesis is a mechanical piece artificially produced with similar geometric characteristics of the natural bone structure. Its production process has been modernized over time due to being supported by digital image processing and CAD tools. From the viewpoint of engineering, prosthesis can be considered as a complex product-feature because of its own geometry as well as because of all the processes required for its manufacture. Despite specific anatomical characteristics of individuals, it is possible that its production follows standard steps as long as it is characterized by clearly defined Received October 15, 2015 / Accepted November 20, 2015 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 Original scientific paper 284 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR methods. The basis of the production process is the geometric modeling in the early stage of design. The construction of a virtual 3D model allows us to evaluate the prosthesis from medical perspective in order to plan a surgical procedure as well as from product development perspective in order to plan its manufacturing process. This process involves different actions starting from the CT image acquisition and respective problem identification (failure in the bone) via the virtual model generation to, finally, the prosthesis production (real model). According to [1] the bones have the similar semantic representation of free-form objects, where the semantic features can contribute in the parameterization of bones‟ models. This approach solves the subjectivity to achieve custom pieces and overtakes the prosthesis modeling limitations as before addressed in [2, 3]. At the level of implementation, the key point is the parameterization of shapes. Recent researches have proposed the automation of bone‟s modeling process in some steps, as example, the sub pixel cubes applied in [4], the Ellipse Adjustment Algorithm (EAA) proposed by [5] and its modified approach by SuperEllipses in [6], as well as the anatomic features generation experienced by [7] through a reverse engineering approach. The well known methods like Splines [7, 8, 9] and Active Shape Models [10, 11] are frequently used to describe regions, and they can be, for instance, adapted to represent bones. Additional difficulties in the bone modeling are present over the open edges problems due to the missing information. In order to overcome limitations while dealing with open contours the new approaches have been proposed such as the ellipse adjustment by [5], which is focused on finding an ellipse with a similar curvature to the bone edge. Thus, a circumference fits the bone edge over a specific CT slice in order to represent an incomplete edge. The main limitation is that the bone contour is not always a regular circumference and, for some slices, the adjustment is weak. The recent research about prosthesis modeling is oriented in the same direction with the aim of finding new methods for the generation of the anatomically-adjusted prosthesis models. Among these computational processes, it is possible to classify the methods due to strategy in 3D or 2D approaches. From 3D approach, the virtual model is built straight on 3D image, through manual intervention where the prosthesis is drawn on its own 3D surface like in [12, 13, 14] or generated by specialized algorithms that mathematically create a model of implant as in [4]. From 2D image, the virtual model is created from the tomographic cuts. A failure in the skull surface generates an open contour in the slices of respective failure region. In this case, all the open sections on the slices are closed one- by-one through image processing methods in order to complete the open curvature in those bone images, as in [5, 6, 16]. Beyond these different approaches in modeling, we can also classify the reconstruction process according to the symmetry criteria. The failure position can be symmetric or non-symmetric. For the symmetric case, the most evident hypothesis is that the healthy side of skull can be mirrored [16] to the appositive side in order to cover the failure region by using the symmetry of body [17]. However, this is only possible when the failure region occurs in the right or left lateral position considering a sagittal cut. In the case of 3D surface reconstruction tools, as in [12], for instance, it is possible to add or remove image fragments from the mirroring region to build the prosthesis piece by a handmade management of the Boolean operations; surely the finest adjustment is also required in all the link positions between prosthesis and skull image. The mirroring can be automated as experienced by [16] through a slice by slice procedure in order to avoid manual intervention, but it also addresses problems in the Autonomous CT Replacement Method for the Skull Prosthesis Modelling 285 linkage frontiers between the model and the skull. In the non-symmetric case, the failure sometimes occurs in the frontal region as well while the symmetry-based process cannot be applied because the same proportion of the known image to mirroring does not exist. In this case, the closing of contour can be performed by computational methods as curve adjustment applied over the CT images. For this situation, some different techniques as Fourier Curves from [18, 19], and Active Shape Models [10, 11], can be used to graphically represent a closed contour region but they have restrictions when it comes to representing open sections because they cannot be adjusted on the missing area. Our proposal is to perform an autonomous content retrieval search in a lot of medical image database in order to find healthy bones with shape-similarity to replace missing regions. The method is based on a CT-by-CT approach. The problem concerned with content retrieval is the amount of data to be compared; then our approach is to define a small set of parameters, i.e. descriptors of the curvature shape. We are applying the Cubic Bezier Curves because its principle is the basis for other curve adjustment methods as B- Splines for instance [8, 20]. The adjustment is evaluated by ABC optimization algorithm [21-25] as a useful tool for evaluating the best of the descriptors parameters. 2. THE PROPOSED METHOD The method implies the determination of the control points of a Cubic Bézier Curve fitted on the bone edge curvature [25]. The control points are a set of shape descriptors because they permit us to build a synthetic edge similar to the original contour. The whole method‟s representation is presented in Fig. 1. modelling of missing region Failure Position Identification CT Slices Segmentation Parametric Modeling of curvature Curvature descriptors Contour Sectioning ABC optimization (Fitness Evaluation) Bezier curve adjustment Search on Image database create the 3D model create the STL file Find compative skulls with similar descriptors Fig. 1 Flow of operation of the proposed method 286 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR The operation flow in Fig. 1 shows the method as having three main layers: image processing, curve descriptors calculation and missing region modeling. This research addresses the curve descriptors calculation. 2.1. Curve descriptors calculation The process starts with the image processing. All the CT slices must be extracted from DICOM file [26]; this is followed by a threshold process [27]; it is an image processing method to perform the segmentation, i.e. transformation of the original grey level image into an equivalent binary (black & white) image. Thus, it enables us to separate the interesting objects of image from its background. Within this procedure, we are able to access the bone information only, meaning separated from the remaining tissues. Once the bone edge has an outer side and an inner side with their respective shapes we need to apply an edge detector to separate them. After this procedure all the selected CT slices contain binary data whose value „1‟ is a pixel of inner and outer edges and „0‟ is the background. The curve calculation is compounded by two processes, where the first is the Failure Position Identification whose aim is look for the slices with an open contour; and the second is the Parametric Modeling of Curvature whose objective is to determinate its descriptors. 2.1.1. Failure position identification Once these images are processed, we need to find only the CT slices whose bone edge comprehends the failure region. An automatic procedure is proposed by [5] in order to perform the selection of the set of slices with uncompleted edge. The method operates through the polar coordinates scanning the edge‟s pixels on each CT slice, and it produces a relationship between an angle theta versus the radius (distance from center to pixel) as presented in Fig. 2.a. Thus, when a pixel exists in that theta direction we get its polar coordinate, and otherwise when no pixel is present (gap position) an infinite radius value is found; then it is possible to mark the initial and final position of failure in each image. Following this, from the set of open slices, we proceed to the modeling for each one. 2.1.2. Sectioning of the edges The Bézier method requires a set of control points to produce an adjusted curve; in this case, some of them are known while the others must be calculated. The known control points P belong to the edge. Our method applies the edge sectioning procedure to sub-divide the edge in small sections Ck. Each section defines two initial parameters for curve generation. Fig. 2.b presents the sectioning on edge where each pair [Pk Pk+1], with k=1,..,n, defines the range of n divisions. The value of n is experimentally determined by Gumiel in [25] and in this example, the n value represents 10 cut sections where the sections are C0=[P0 P1], C1=[P1 P2], …, C9=[P9 P10]. The number of descriptors necessary to represent the curvature implies an equivalent number of sections. The difficulty is to determine the balance between n value and the precision required because a higher number of divisions (n) cause a big set of parameters and, consequently, a big amount data for processing. Otherwise, small n causes weak adjusted curves due to the lack of information. Thus, the main problem is to find equilibrium between the best accuracy of Autonomous CT Replacement Method for the Skull Prosthesis Modelling 287 the results with the minimum quantities of descriptors. This optimization is presented in the following section 2.2.4. Fig. 2 a) Polar coordinate mapping to identification of the limits in missing region b) Sectioning of contour for fixed control points (Pk) determination 2.1.3. Parametric modeling of curvature The main objective in this phase is to generate a synthetically built curve over original bone edge in order to describe the bone curvature in each CT slice. Thus, we get a small set of points to represent the same shape of bone curvature instead those all of edge pixels. The parameters to be calculated by the Bézier Method are presented in Fig. 3. Fig. 3 a) A sample of a Bézier Method showing the curve defined by its respective control points; b) The sample of a curve generate on a piece of bone edge and respective position of control points The aim is to define the “variable control points” that lie between two known fixed points P1 and P2. As presented in Fig. 3.a, the fixed control points P1 and P2 define the beginning and the end of the Bezier curve in a section, and its shape is given by p3 to p10. By general form we have: 1 2 3 4 5 6 7 8 9 {[ ];[ , , , , , , , ]} Ck k k k k k k k k k k B P P p p p p p p p p           (1) a) b) a) b) 288 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR Where BCk is the respective Bezier curve for section Ck. Ps are known and ps are unknown. The Bézier technique needs a “control polygon” in order to generate the curvature whose vertices are the set of control points. As presented in the example in Fig. 3.a, the Bezier curve is built from four control points so called P1, P2 (as the fixed points lying on edge) and p3, p4 (as variable points) whose coordinates define a “control polygon”. The polygon lines determine a base support line for the variable control points (p5, p6, p7). Changing those points, p8, p9 are moved along the support line and define p10. This point is the resulting point and it defines the pixel position to draw a curve. By changing all variable p points it is possible to obtain all curve points. Fig. 3.b presents an example of the Bezier curve adjusted (blue line) on the bone edge to the outer contour (green line). The quality of adjustment depends on the position of all control points, and the challenge is to calculate the possible coordinates of the “variable points” (red points) once the fixed points are known because they are on the bone edge. In Fig. 4.a there is a good Bezier adjustment and in Fig. 4.b a sample of a wrong adjustment. Fig. 4 Examples of Bezier Curves adjusted on edge Fig. 4 shows in (a) fixed points P1 and P2 and respective variable points p3 and p4. In these positions, the Bezier curve is on the original edge. In Fig. 4.b, P1 and P2 are in the same position as (a) but p3 and p4 are not correctly positioned, and then by consequence the Bezier curve does not represent the original edge curvature. The process of finding the best position for the variable points for each BCk curve is determined by application of the Artificial Bee Colony Optimization (ABC) as described further on. 2.1.4. The Artificial Bee Colony optimization (ABC) The ABC optimization algorithm proposed by Karaboga [21, 22] or its modified approach by Abro [23, 24] is based on concept of collective behavior of insects in nature (bees, for instance), modified to solve optimization problems. It is an algorithm of swarm intelligence (SI) where the individuals interact locally among them in order to solve a problem (global objective) [28]. In this case, we are applying the ABC algorithm in order to determinate the coordinate of the variable control points (p). Karaboga et al [22] have compared the ABC performance face the GA (Genetic Algorithm), DE (Differential Evolution), PSO (Particle Swarm Optimization) and demonstrated its applicability. Because of the experience acquired in this research, the ABC approach is selected as more suitable for our optimization problem. a) b) Autonomous CT Replacement Method for the Skull Prosthesis Modelling 289 The ABC algorithm is based on the bees searching for a source of food [21], and in our context, the “source of food” location is a possible solution of the variable control points. In addition, the amount of nectar corresponds to the solution‟s quality, i.e. the quality of the calculated Bezier curve. The quality of the calculated curve is evaluated by fitness function, f, as presented in Eq. 2. 2 2 1 ( ) ( ( ) ( )) ( ( ) ( )) n B O B O i f k X i X i Y i Y i      (2) Fitness function f(k) is calculated for each section Ck=1,2,...,L with L as the total number of sections of the edge, and n as the maximum number of control points for comparison. The value of f defines the difference between calculated values XB and YB (values of coordinates X and Y from the Bezier curve) and those of all XO and YO (values of original contour pixels). Following the ABC algorithm from [22], the initial population is a set of the possible solutions randomly generated and in our notation Pk and Pk+1 are respectively the starting point (first point) and ending point (second point) of the section on the bone edge. Then, they form a rectangular polygon that limits a searching area in order to look for the variable points (p3 and p4) as in Fig. 5.a. P1 and P2 are known because of lying over the edge at the cut position; the challenge is to find the best p3 and p4 candidate points whose coordinates are more suitable to produce the Bezier curve in each section. In Figs. 5.b and 5.c, they present a sample of the optimization result applied on contour through Ck=10 sections. In total, we have about 20 variable points from 11 fixed points. Fig. 5 Descriptors calculation a) Searching area for section Ck; b) Fixed and variable control points; c) Descriptors to the inner and outer contour The fitness of curvature is evaluated by a function f that measures the difference between the original edge and the generated curve. The best values after the convergence of method are selected as answer points. Then, the total of points {PBCk} is the curvature descriptors. a) b) c) 290 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR 3. RESULTS AND DISCUSSION For demonstrating our purpose, a synthetic irregular failure is hand-made built through ImageJ plugging from [16], [30]. Thus, the original cut piece is to be compared with modeled one. From image of Fig. 9, a total of 26 CT slices with their respective open contours are extracted as samples in Fig. 9.b. Now, n=10 sectioning of the edge gives us 20 set of descriptors. The process is repeated for images in database and the parameters are calculated again for each one. The compatible CT slices are extracted from the image database within CTs set about 25 patients. The selected are all whose descriptors which match with searching parameters as shown in Table I. Fig. 6 Defective skull and respective failure region As presented in Fig. 6, the skull‟s failure region can be decomposed by its respective CT slices. In Fig. 6.a an example of a frontal failure is presented. Note that in this case the modeling does not apply symmetry by mirroring, like show images I, II and III in Fig. 6.b. Otherwise, only the image IV could be solved by mirroring. Then the proposed method attends both symmetric and non-symmetric cases. Table I shows the recovered compatible slices for each testing CT image from patient#1. Note that the prosthesis range can be replaced by compatible information from 3 distinct patients (patients #05, #06, and #10). These three have similar morphological characteristics in the specific tomographic cut. The best CT is selected within minimum error for each tested slice. Due to the fact that the ROI position is known on the original image, the same part of CT can be extracted from compatible Cts. A sample of some compatible pieces is presented in Fig. 7. For example, the first slice number #249 from patient #06 can replace the slice number #249 in patient #01, the second slice number #250 in patient#01 can be replaced by the slice number#245 from patient#05, and so on. All of the recovered slices can be superimposed to form a 3D piece as presented in Fig. 8. The new 3D piece fills the failure region as presented in Fig. 9. a) b) Autonomous CT Replacement Method for the Skull Prosthesis Modelling 291 Table 1 Compatible slices whose descriptors match with testing images Testing Image (CT-patient #1) Compatible CT Image (CTs recovered) Fitness (similarity measure) Slice # 249 Patient n o 06, slice n o 249 31,0711 Slice # 250 Patient n o 05, slice n o 245 63,5262 Slice # 251 Patient n o 05, slice n o 247 73,5553 Slice # 256 Patient n o 06, slice n o 246 76,5050 Slice # 257 Patient n o 06, slice n o 245 81,3127 Slice # 258 Patient n o 05, slice n o 247 98,6153 Slice # 259 Patient n o 10, slice n o 249 92,8940 Slice # 267 Patient n o 10, slice n o 248 166,4766 Slice # 268 Patient n o 06, slice n o 244 174,0587 Slice # 269 Patient n o 10, slice n o 244 171,3661 Slice # 270 Patient n o 05, slice n o 274 81,0014 Slice # 271 Patient n o 06, slice n o 271 47,4940 Slice # 272 Patient n o 06, slice n o 273 44,2646 Slice # 273 Patient n o 05, slice n o 275 46,2553 (249, 6) (247, 5) (246, 6) (249, 10) (248, 10) (244,6) (274,5) (273, 6) Fig. 7 The samples of compatibles pieces that match the original missing region, labeled here as („slice#’, „patient#’) Fig. 8 Three views of modeled prosthesis 292 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR Fig. 9 The piece anatomically adjusted on the ROI Fig. 10 Measurement of error between original skull and the replaced region evaluated through Geomagic® system [29] In Fig. 9, there is a visual sense of the prosthesis piece covering the missing region on the testing skull image. The linkage between the prosthesis piece and the skull seems as discontinuity because the 3D view tool (from ImageJ/Fiji plugging [30]) caused shadows on the edges of these two objects. The better visualization is given in Fig. 10 by presentation of a colored scale for 3D error measured by Geomagic® system [29]. The values show the differences among ≈1mm until ≈2mm between the original skull and the replaced region inside the prosthesis range. 4. CONCLUSION AND FUTURE WORK The paper addresses the prosthesis modeling problem. It presents a method to fill a defective region on the skull by a CT replacement strategy. The proposed method is compounded by activities involving image processing, optimization and curves adjustment. The aim of proposal is the creation of the contours‟ descriptors based on the Cubic Bezier Autonomous CT Replacement Method for the Skull Prosthesis Modelling 293 Curves. The concept behind the method is to perform a content retrieval from a medical image database in order to look for compatible CTs slices from other patients with similar morphological characteristics. The good compatible slices with a closed contour are used to replace the defective slices with an open contour. Thus, a piece of image from a good bone edge is superimposed on the position of the missing area of testing image. All the processes for descriptor‟s generation and its respective geometric prosthesis modeling are performed automatically. A testing image is prepared for demonstration. An irregular hole is handmade built on the frontal position of a skull image. In this way, we have enough information to compare the modeled prosthesis with the real bone image. The Geomagic software is applied to evaluate the difference in measure. The presented result shows that it is possible to create a respective 3D prosthesis model with minimum error without the user‟s intervention. From the engineering point of view, in terms of geometry, we have demonstrated that it is possible to create anatomical bone models. However, for future work, it is necessary to complete the other phases of process such as printing as well as all medical evaluation. This proposal is an important step to prosthesis modeling in our ongoing study, but in order to comply with further practical application, it is also necessary to enlarge the testing image database to improve the selection of more compatible answers. Acknowledgements: The author would like to thanks the Graduate Program in Health Technology – PPGTS from Pontifical Catholic University of Parana – PUCPR by collaboration in providing all CT image data. 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, pp.1-20 (DOI:10.1017/S0890060415000153). 2. Huang, G. Y., Shan, L.J.; 2011, Research on the Digital Design and Manufacture of Titanium Alloy Skull Repair Prosthesis. Proceedings of the IEEE 5th Internacional Conference on Bioinformatics and Biomedical Engineering (ICBBE), Wuhan, China, pp.1-4. 3. Saldarriaga, J. F. I., Vélezes, S. C., Posada, A. C., Henao, B. B., Valencia, C. A. T., 2011, Design and Manufacturing of a Custom Skull Implant. American J. of Engineering and Applied Sciences, 4(1), pp. 169-174. 4. You, F.; Hu, Q.; Yao, Y.; Lu, Q., 2010, A new modeling method on skull defect repair. IEEE International Conference on Measuring Technology and Mechatronics Automation, 1, p. 568-572. 5. Greboge, T., Rudek, M., Jahnen, A., Canciglieri, O., 2013, Improved Engineering Design Strategy Applied to Prosthesis Modelling, Product Service Engineering in a Dynamic World. Led. Amsterdam: IOS Press, 1, pp.60-71. 6. Rudek, M.,Canciglieri, Jr. O., Jahnen A., Bichinho, G. L., 2013, CT slice Retrieval by Shape Ellipses Descriptors for Skull Repairing. The IEEE Internacional Conference on Image Processing (ICIP 2013), pp. 761-764. 7. 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. 8. Goldman, R., Lychie, T., 1993, Knot Insertion and Deletion Algorithms for B-splines Curves and Surfaces. Philadelphia: SIAM, 197 p. 9. Egerstedt, M., Martin, C., 2010, Control Theoretic Splines: Optimal Control, Statistics, and Path Planning, Princeton University Press, 217 p. 10. Cootes, T.F., Edwards, G.J., Taylor, C.J., 1998, Active Appearance Models. H. Burkhardt and B.Neumann, editors, 5th European Conference on Computer Vision, Springer, Berlin v.2, pp. 484-498. 294 M. RUDEK, Y. B. GUMIEL, O. CANCIGLIERI JR 11. Van Ginneken, B., Frangi, A.F., Stall, J.J., Romeny, B.M., Viergever, M.A, 2002, Active Shape Model Segmentation with Optimal Features. IEEE Trans. Med. Imag., 21(8), pp.924-933. 12. Blender, 2015, Creating a Craniofacial Prosthesis with InVesalius and Blender, available at: https://www.youtube.com/watch?v=_PUzeToMVoU (last access: Nov. 2015.) 13. Oliveira, M.A.F., Foggiato, J.A., 2012, Modelagem de Superfícies a Partir De Tomografias Computadorizadas de Partes Anatômicas Para Uso em Programas de CAD-3D. SICITE XVII Seminário de iniciação científica e Tecnológica da UTFPR, (in portuguese). 14. Kazhdan, M., Bolitho, M., Hoppe, H., 2006, Poisson Surface Reconstruction, Eurographics Symposium on Geometry Processing, pp.61-70. 15. Rouhani, M., 2015, Surface Reconstruction using Implicit B-Splines, available at: http://www.mathworks.com/ (last access: Nov. 2015.) 16. Rudek, M., Mendes, G. C., Jahnen, A, 2015, Failure-Correction Simulation Tool Applied to Skull Prosthesis Modelling, Proceedings of 5th International Conference on Information Society and Technology - ICIST 2015, Kopaonik – Serbia, pp. 1-6. 17. Li, H., Xie, Z., Ruan, S.,Wang, H., 2009, The Measurement and Analyses of Symmetry Characteristic of Human Skull Based on CT images. Biomechanics and Vehicle Safety Engineering Centre, Tianjin University of Science and Technology, China, 26(1), pp. 34-37. 18. Zahn, C. T.; Roskies, Ralph Z., 1972, Fourier Descriptors for Plane Closed Curves, IEEE Transactions on Computers, C-21(3), pp.269-281. 19. Yoshinori, U., 1984, A new Fourier descriptor applicable to open curves, Electronics and Communications in Japan (Part I: Communications), 67(8), pp. 1-10. 20. Shao, L.J. , Zhow, H., 1996, Curve Fitting with Bezier Cubics, Graphical Models and Image Processing, 58(3), pp.223-232. 21. Karaboga, D., Akay, B., 2009, Algorithms Simulating Bee Swarm Intelligence, Artificial Intelligence Review, 31(1), pp.68-85. 22. Karaboga, D., Gorkemli, B., Ozturk, C., Karaboga, N., 2012, A comprehensive survey: artificial bee colony (ABC) algorithm and applications, Artificial Intelligence Review, 42(1), pp. 21-57. 23. Abro, A.G., Mohamad-Saleh, J., 2012, An Enchanced Artificial Bee Colony Optimization Algorithm, Recent Advances in Systems Science and Mathematical Modelling, Ed. WSEAS Press, pp.222-227. 24. Sulaiman, N., Mohamed-Saleh, J., Abro, A.G., 2013, A Modified Artificial Bee Colony(JA-ABC) Optimization Algorithm, Proceedings of The International Conference on Applied Mathematics and Computational Methods in Engineering, Rhodes Island, Greece, pp.74-79. 25. Gumiel, Y. B., 2015, Modelagem de Próteses Anatomicas baseada em Descritores Definidos por Curvas de Bézier Cúbicas, Master Thesis, (in Portuguese) Pontifical Catholic University of Parana, Brazil, 108p. 26. Dicom, 2004, Digital Imaging and Communications in Medicine (DICOM) Part 5: Data Structures and Encoding, available at: http://dicom.nema.org/dicom/2004/04_05pu.pdf, (last access: Nov. 2015.) 27. Gonzalez, R. C., Woods, R. E., 2008, Digital Image Processing. 3ª ed. Prentice Hall, 954 p. 28. Jevtic, A., La Fuente, D., 2007, Swarm intelligence and its applications in swarm robotics, The 6th WSEAS Int. Conference on Computational Intelligence, Man-Machine Systems and Cybernetics, Tenerife, Spain, pp. 41-46. 29. Geomagic, 2015, Modeling Easter Island’s Moai with Geomagic 3D Scan Software, available at: http://www.geomagic.com/en/) (last access: Nov. 2015.) 30. Schindelin J., Arganda-Carreras I., 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.