Evaluation of geometric deviations in rapid prototyped three-dimensional models created from computed tomography data

Full Text

BACKGROUND

Medical prototyping technologies, which combine computer-aided design (CAD) and three-dimensional (3D) printing methods, have been increasingly used in various fields of clinical and basic medicine in recent decades [1–4]. Particularly, this was facilitated by the widespread use of additive manufacturing facilities since the 2000s, following the launch of the replicating rapid prototypes and the patent expiration for the fused deposition modeling technology.

The medical prototyping success in the clinical practice depends on the maximum compliance of the resulting products with the patient’s anatomy [5]. Thus, in some cases (printing personalized implants and creating guiding resection templates), the degree of congruence between the printed medical product and the bone surface may directly influence the surgical treatment outcome [6]. In other situations (preoperative assessment and creating training phantoms), accurate reproduction of geometric relations between the normal and abnormal tissues predicts possible complications during surgery [7].

Obtaining the reference values for the dimensions of the studied anatomical structures is difficult due to the barrier in the accuracy assessment of the resulting models in vivo. Direct measurement of the object of interest is not always possible even during the surgical intervention (due to organ deformations), and obtaining 3D models and assessing their inaccuracies are based on the medical imaging methods used for the noninvasive collection of morphological characteristics of organs. Alternatively, a comparison of models derived from ex vivo scans of bone structures may be used; however, this approach uses a relatively small number of checkpoints, and the complexity of their interpretation by different specialists may lead to errors [8]. Tests on 3D models of large mammals (e.g., pigs) give good results with their subsequent introduction into the clinical practice [9].

Physical models were selected as objects of study in earlier works devoted to the accuracy assessment of model reproduction by additive manufacturing tools. These models were measured by several checkpoints using direct means: ruler, caliper, and coordinate measuring machine [10–12]. The present study proposes to use a fundamentally different approach to assess geometric deviations, which consists of all polygonal mesh points comparison obtained using additive technologies with a reference representing the original 3D model. Therefore, a model of complex shape and predetermined dimensions designed using parametric modeling was used as a reference. Conversely, the compared object was a virtual model of a reference passed three main stages of medical prototyping (3D printing, scanning, and modeling). This approach simultaneously assessed the deviations over the entire surface of the product, without the need for numerous linear measurements, while reducing the human factor influence on the measurement process.

This study aimed to assess the geometric deviations that occurred when reproducing objects imitating bone tissues using 3D modeling (based on computed tomography [CT] data) and additive manufacturing.

MATERIALS AND METHODS

The study design is shown in Fig. 1.

 

Fig. 1. Study design.

 

In the FreeCAD¹ parametric modeling program (FreeCAD Community, Germany), a solid model based on a cuboid with dimensions of 20 × 20 × 12 mm (l × w × h) was designed, which contained five end-to-end parallel channels with diameters of 1, 2, 3, 4, and 10 mm. Two recesses and two eminences in the form of a hemisphere and a cone (imitating rounded and sharpened notches and protrusions on the bone surface) at 4 mm in height were created on the upper edge of the model (hereinafter—reference 1). The diameter of the channels was selected to simulate the different foramina of the human skull on one product, which are visualized by modern CT methods. A 1-mm step for smaller diameter channels and irregularities of a given shape on the product surface was applied to manually assess the quality during model printing. Using the Blender² software package (Blender Foundation, the Netherlands), volumes were calculated and a copy of this reference, doubled in length, width, and height (hereinafter—reference 2), was created to check the effect of size increase on linear deviations.

Parametric models saved in Standard Triangle Language (STL) format were uploaded to the Repetier Host³ program (Hot-World GmbH & Co. KG, Germany), where the G-code file for the 3D printer was generated using the CuraEngine⁴ slicer (Ultimaker, Netherlands) with the following printing parameters: 0.2 mm layer height, 0.8 mm wall thickness, 33% filling (selected empirically), 50 mm/s speed, 210°C nozzle temperature, 50°C platform temperature, forced model blowing, 5 mm retract, and 100% filament flow. The printing was performed on an Ender-3 3D printer (Creality3D, China) with 0.4 mm nozzle diameter using the BFCopper PLA plastic (Best Filament, Russia) containing copper particles to simulate the X-ray density of bone tissue (average X-ray density of plastic at 100% filling was +1762 HU, σ = 172 HU).

The resulting products were scanned using a 128-slice Ingenuity CT scanner (Philips, the Netherlands) in the air with channels oriented perpendicular to the gentry plane. The X-ray tube voltage, current, slice thickness, and pixel size of the reconstructed slices were 120 kV, 117 mA, 0.625 mm, and 0.43 × 0.43 mm, respectively.

Series of tomographic images in Digital Imaging and Communication in Medicine (DICOM) format were loaded into 3D Slicer⁵ software (3D Slicer Community, US). The images were used to create STL models of reference 1: four by automatic voxel selection with the threshold values of +500 HU, 0 HU, −500 HU, and −750 HU, respectively (Threshold Paint tool), and one by manual slice voxel tracing (Paint tool). In addition, one model of reference 2 was created by automatic selection of voxels with values > −500 HU. The threshold values for reference 1 were selected empirically based on the X-ray density values in the outer layer of the model (approximately 0.9 mm thick) as −1000…+500 HU (−1000 HU was replaced with −750 HU to exclude the ambient air from the model). Low parameters (−500 HU and −750 HU) were deliberately chosen due to the pronounced defects when using positive values of the segmentation threshold. The density threshold for reference 2 was randomly selected from the thresholds used for reference 1.

The resulting models were loaded into the Meshmixer⁶ program (Autodesk, USA), where polygonization errors were analyzed and corrected, structures not in contact with the outer shell of the models were removed, and the model mesh was rebuilt with a fixed polygon edge length (Remesh-Target Edge Length tool) equal to 0.25 mm. The appearance of the model at each of the listed stages is shown in Fig. 2.

 

Fig. 2. Appearance of the model of reference 1 after the stages: a - parametric modeling; b - 3D printing; c - CT scan (axial slice, window level +805 HU, window width 3718 HU, higher density of protrusions due to closer filament placement in the horizontal plane); d - Creates a polygon mesh based on CT data.

CT: computed tomography.

 

At the final stage, the models obtained from CT data, together with their references, were loaded in pairs into CloudCompare⁷ software (CloudCompare Project, France), where linear deviations of the final CAD model were calculated based on the iterative closest point (ICP) algorithm. The calculation was performed from each vertex of the obtained polygonal mesh along the normal to the nearest surface of the reference. Statistical analysis was performed using the GNU PSPP⁸ program (Free Software Foundation, US).

RESULTS

Model volumes were determined for a preliminary assessment of their shape distortion. Particularly, if the volume of the rebuilt polygonal mesh was less/more compared to that of the reference CAD model, the linear dimensions of the analyzed model were expected to decrease/increase. In the case of equal volumes, either a size match or a compensated distortion was expected. The volumes of references and models obtained by different methods of segmentation are shown in Table 1.

 

Table 1. Volumes of virtual models

Model name, cuboid dimensions, mm (segmentation type and threshold)	Volume, mm3	Differences with the reference, mm3
Reference 1	3576	-
20×20×12 (auto +500)	3607	31 (0,9%)
20×20×12 (auto 0)	3901	325 (9,1%)
20×20×12 (auto -500)	4255	679 (19%)
20×20×12 (auto -750)	4480	904 (25,3%)
20×20×12 (manual)	4538	962 (26,9%)
Reference 2	28 608	-
40×40×24 (auto -500)	31 140	2532 (8,9%)

 

Polygonal meshes processed with the Remesh tool have an approximately equal density of polygonal vertices distribution per surface unit. When comparing the models (combining cuboidal vertices) with the references, the resulting data set was the number of polygonal mesh nodes removed from the surface of the CAD model at a certain distance (mm) in the outer (positive values) or inner (negative values) directions. An example of aligning the model and the resulting histogram of surface point deviation is shown in Fig. 3.

 

Fig. 3. Alignment of polygonal meshes of models (a) and a histogram of the calculated deviations of linear dimensions (b).

 

The additional peak in the region of positive values is due to the anisotropy of the voxels along the Z-axis, as well as a slightly excessive material application on the side faces of the cuboid.

Obtained values of distorted geometric dimensions for each model are summarized in Table 2.

 

Table 2. Linear deviations of model dimensions

Model name	Minimum value, mm	Maximum value, mm	Average value, mm	Standard deviation, mm	Integral sum, mm
20×20×12 (auto +500)	-0,58	1,146	0,026	0,214	1904
20×20×12 (auto 0)	-0,533	1,019	0,141	0,223	9532
20×20×12 (auto -500)	-0,421	1,129	0,296	0,211	20 756
20×20×12 (auto -750)	-0,675	1,107	0,373	0,197	27 179
20×20×12 (manual)	-0,809	1,068	0,411	0,253	18 190
40×40×24 (auto -500)	-0,862	1,353	0,37	0,275	50 213

 

The number of measurements performed by the program corresponded to the number of vertices of the polygonal mesh (approximately 100 thousand), and the distribution of deviations of linear dimensions were close to normal. Based on obtained values, the program built a Gaussian (real Gaussian function), with the argument of maximization, which was used as the average deviation of linear dimensions of this model. All linear deviations calculated by the program were divided by their values into intervals of equal width (classes) to build a histogram. For each model, the value of the integral sum (hereinafter—Sum) of linear model deviations was calculated as follows:

$\mathrm{Sum} = \sum _{\mathrm{i}}^{\mathrm{n}} = ({\mathrm{d}}_{\mathrm{i}} \times {\mathrm{q}}_{\mathrm{i}})$         (1)

Where d_i is the minimum value of the linear deviation in the i-th class; q_i is the number of representatives in the i-th class; and n is the total number of classes.

Linear normalization by formula (2) was applied to each of the values (sum and volume change) for their bringing to a dimensionless form.

$\widetilde{{\mathrm{\chi }}_{\mathrm{i}}}=\frac{{\mathrm{\chi }}_{\mathrm{i}} ={\mathrm{\chi }}_{\mathrm{i}}– {\mathrm{\chi }}_{\mathrm{i} \min } }{{\mathrm{\chi }}_{\mathrm{i} \max }-{\mathrm{\chi }}_{\mathrm{i} \min }}$               (2)

The linearly normalized sum values were compared with the volume change indices of the corresponding models to check the shape distortion direction (Fig. 4).

 

Fig. 4. Linear normalized values. Red color: differences in the volume of models with the reference; blue color: differences in the integral sum of linear deviations.

 

A comparison of the sums of linear deviations and relative changes in the model volumes using Spearman’s rank correlation coefficient showed a high level of correlation (ρ = 0.83, t_emp = 5.27, and p = 0.05). Therefore, the change in the reproduced model volume corresponds to the linear deviations measured by the used software tools.

DISCUSSION

The process of creating a medical prototype includes three main stages: data aquisition (scanning), data processing (creating a virtual model), and 3D printing [13].

The first stage is a radiological examination, which obtains data on the 3D structure of the region of interest with a high spatial resolution (CT or magnetic resonance imaging or 3D ultrasound scanning), results of which are saved as discretized images in DICOM format. The change in geometry at this stage may be due to the specifics of obtaining and processing diagnostic information by the selected imaging method.

The second stage involves segmentation (selecting voxels of medical images related to the created model), voxel mesh into a polygonal mesh conversion, and resulting 3D model editing. Segmentation may be of three types: manual (fully performed by the operator), semiautomatic (performed by the computer and corrected by the operator), and automatic (fully performed by the computer) [14]. Automated methods are more attractive due to lower labor costs; however, because of characteristics of used computer algorithms, the geometry of the final product may be severely distorted, and their accuracy requires a separate study [15]. At this stage, a specialist with knowledge in anatomy (particularly, radiological anatomy) should be involved to correctly select and edit the objects of interest. Thus, geometry distortion at this stage mainly results from the human factor.

The third stage is the 3D printing of the resulting virtual model. For its successful implementation, the model is pre-processed (slicing) and post-processed (removal of supports and surface handling) after printing (if necessary). The contribution of this stage to the distortion of the final product geometry primarily depends on the type of used additive equipment.

Therefore, to assess geometric distortions in the prototyping process, all three stages must be performed on some physical phantom with known linear dimensions and defined design elements. Scanning conditions and printing settings may affect the accuracy of the resulting product [10, 16]. Thus, the relevant parameters were kept at the same level for all produced models.

Many different software packages can be used to segment medical images. The use of 3D Slicer at the stage of virtual model creation from the DICOM data was due to its accessibility (distributed as open-source software) and a large number of additional modules and extensions that make it an ideal tool for preoperative planning [17].

The results confirm the “dumbbell” effect described earlier, according to which a decreased automatic segmentation threshold leads to an expansion of resulting model outlines [8]. Thus, when the threshold was decreased from +500 to −750 HU, the average value of the surface point displacement from the reference consistently increased from 0.026 to 0.373 mm, respectively. These distortions depend more on the segmentation method compared to the size of the scanned object, since the average linear deviations were 24% higher in all dimensions when the model was proportionally doubled compared to the original size models with the same threshold, and 10% lower as opposed to the original size of manually segmented models.

The relatively high values of geometric deviations in the manually segmented model (by 0.41 mm for a 20 × 20 mm sample) may result from the operator’s aligning the side faces of the cuboid along the edges, which were somewhat displaced outward (during 3D printing). Moreover, similar deviations were observed during semiautomatic segmentation. Remarkably, these models reproduced the original shape better due to the preservation of all control elements and the absence of wall defects (Fig. 5).

 

Fig. 5. Appearance of models segmented by semi-automatic with a cutoff threshold of 0 HU (a - with a map of deviations from the standard, b - general view) and manually (c - with a map of deviations from the standard, d - general view) method.

 

The visual assessment of models that are automatically segmented with a threshold of 0 and +500 HU observed 1.2 mm thick wall defects, whereas models with −750, −500, and 0 HU thresholds lacked created end-to-end holes with a diameter of 1 mm. The defects with a width exceeding the module of minimum geometry deviations are explained by the thickness measurements performed from two surfaces, each with appropriate deviations (Fig. 6).

 

Fig. 6. Measurement of linear deviations from the standard (blue lines) in the defect area of the compared model (red lines).

 

Thus, the linear displacement summation of the vertex relative to the outer surface (d₁ segment) with the nearby vertex displacement relative to the inner surface (d₂ segment) may correspond to a wall defect of l thickness exceeding the maximum linear deviation value. Similarly, the positive deviation summation leads to the “obliteration” of holes that exceed the maximum size deviation calculated from the ICP algorithm results. Thus, the linear dimensions of the reproduced models include two values of linear displacement.

In addition, channels with a diameter of 1 mm, which were not detected on the mentioned copies during the visual assessment, were not evaluated by the algorithm since the polygonal vertices corresponding to these channels were absent on the models.

Study limitations

The limitations of this study could be the relatively small sample size and the use of materials that did not correspond to the true bone tissue composition.

CONCLUSIONS

Distortion of the original shape inevitably occurs in the creation of a product using medical prototyping (CT scanning, 3D modeling, and 3D printing). The change in the model volumes and the average linear deviations of their surface points from the reference, determined using the ICP algorithm, have positive values. The distortion depends more on the selected segmentation method compared to the overall scale of the model or its parts. The use of the manual segmentation compared to the semiautomatic method leads to a slightly higher distortion of linear dimensions resulting from aligning the model to the unreliable landmarks; however, it reproduces all elements of the original sample. Thus, for necessary visualized anatomical structure preservation, the resulting virtual models should be corrected by a specialist. In this case, the use of different thresholds in areas with dense and sparse structures may be sufficient for segmentation. Using the semiautomatic segmentation is more preferable if several small structures can be neglected and the main task is to match the overall product dimensions to the simulated organ; however, the thresholds of this method should be selected experimentally depending on the tissue composition and scanning parameters.

ADDITIONAL INFORMATION

Funding source. The authors declare that there is no external funding for the exploration and analysis work.

Competing interests. The authors declare no obvious and potential conflicts of interest related to the publication of this article.

Authors' contribution. Aleksandr V. Shirshin — search for relevant publications, literature analysis, research design development, data processing, writing; Igor S. Zheleznyak — determination of the main focus of the review, expert evaluation of literature review, research design development; Vladimir N. Malakhovsky — expert evaluation of literature review, systematization and final editing of the review; Sergey V. Kushnarev — search for relevant publications, determination of research materials; Nataliya S. Gorina — literature analysis, determination of research methods. All authors made a substantial contribution to the conception of the work, acquisition, analysis, interpretation of data for the work, drafting and revising the work, final approval of the version to be published and agree to be accountable for all aspects of the work.

¹ FreeCAD Your own 3D parametric modeler. Available at https://www.freecadweb.org/. Accessed on 05/15/2021.

² Blender 2.93.2 LTS. Available at https://www.blender.org/. Accessed on 05/15/2021.

³ Repetier. Available at https://www.repetier.com/. Accessed on 05/15/2021.

⁴ Ultimaker Cura. Available at https://ultimaker.com/software/ultimaker-cura/. Accessed on 05/15/2021.

⁵ 3D Slicer image computing platform. Available at https://www.slicer.org/. Accessed on 05/15/2021.

⁶ Autodesk Meshmixer free software for making awesome stuff. Available at https://www.meshmixer.com/. Accessed on 05/15/2021.

⁷ CloudCompare: 3D point cloud and mesh processing software. Open Source Project. Available at http://www.cloudcompare.org/. Accessed on 05/15/2021.

⁸ GNU PSPP. Available at https://www.gnu.org/software/pspp/. Accessed on 05/15/2021.

About the authors

Aleksandr V. Shirshin

Kirov Military Medical Academy; ITMO University

Author for correspondence.
Email: asmdot@gmail.com
ORCID iD: 0000-0002-1494-9626
SPIN-code: 4412-0498
Russian Federation, 6G, Akademika Lebedeva street, Saint-Petersburg, 194044; 49, Kronverksky pr., St. Petersburg, 197101

Igor S. Zheleznyak

Kirov Military Medical Academy

Email: igzh@bk.ru
ORCID iD: 0000-0001-7383-512X
SPIN-code: 1450-5053

MD, Dr. Sci. (Med.), Assistant Professor

Russian Federation, 6G, Akademika Lebedeva street, Saint-Petersburg, 194044

Vladimir N. Malakhovsky

Kirov Military Medical Academy

Email: malakhovskyvova@gmail.com
ORCID iD: 0000-0002-0663-9345
SPIN-code: 2014-6335

MD, Dr. Sci. (Med.), Professor

Russian Federation, 6G, Akademika Lebedeva street, Saint-Petersburg, 194044

Sergei V. Kushnarev

Kirov Military Medical Academy

Email: S.v.kushnarev@yandex.ru
ORCID iD: 0000-0003-2841-2990
SPIN-code: 5859-0480

MD, Cand. Sci. (Med.)

Russian Federation, 6G, Akademika Lebedeva street, Saint-Petersburg, 194044

Natalia S. Gorina

Kirov Military Medical Academy

Email: natali_bgmu@mail.ru
ORCID iD: 0000-0002-6220-8195
SPIN-code: 8175-6746
Russian Federation, 6G, Akademika Lebedeva street, Saint-Petersburg, 194044

References

Supplementary files