Statistical and Metrological Framework for Evaluation of micro-CT 3D Data with Application to Dental Implant Components

A. Components and μCT measurement

The assembled implant and its three separately reconstructed components are shown in Figs. 1–4. Panels (a)–(c) of each component figure show the nominal CAD model, the μCT-derived STL surface, and a color-coded μCT-to-CAD deviation rendering, respectively. These renderings provide spatial context; the color scales are component-specific and do not permit direct numerical comparison between maps. The 3 mm scale bars apply to the component renderings. Because each component was aligned separately to its CAD reference, these comparisons do not preserve the relative positions of the components in the assembly.

Micro-CT measurements of the titanium implant components were performed using a Nanotom 180 X-ray microtomography system (GE Phoenix|X-ray, Germany). The scans were acquired with a voxel size of 13.75 μm and 1800 projections over a full 360° rotation. The X-ray tube was operated at 160 kV and 90 μA in measurement mode 0, without detector binning. A 0.2 mm copper filter was used to attenuate low-energy X-rays and limit beam-hardening effects. The detector exposure time was 750 ms per projection, and three-frame averaging was applied to reduce image noise. These acquisition settings do not, by themselves, establish dimensional accuracy or uncertainty [1, 5].

The projection data were reconstructed into three-dimensional volumes using datos|x 2.1 software (GE Sensing & Inspection Technologies GmbH, Germany) and the Feldkamp cone-beam reconstruction algorithm. The volumes were processed in VGSTUDIO MAX 2.0 (Volume Graphics GmbH, Germany). The same software isosurface threshold setting of 1/2 was applied in VGSTUDIO MAX to the reconstructed volumes of D1, D2, and D3 before STL export. This is a segmentation setting, not a dimensional calibration or an estimate of surface-location uncertainty; the corresponding absolute grey-level cut-off and the surface-determination implementation are not documented in the available exports. No systematic perturbation of the threshold was performed, so the sensitivity of the measured deviations to its choice remains unknown. Applying the same setting ensures a consistent processing choice but does not yield equal segmentation bias across components with different geometries [1, 5].

Fig. 2.

Fig. 3.

Fig. 4.

B. Registration and definition of vertex deviations

Each muCT-measured STL surface was registered independently to its corresponding manufacturer-controlled, component-specific CAD geometry. Its revision was fixed for the comparison and remains identifiable in the project records; the proprietary revision identifiers are not disclosed. The initial alignment used selected nominal reference features for each component, followed by a surface-based best-fit registration. The identities of those features, their construction datum, and the numerical initial and final transformations have not been released with the vertex exports. Consequently, a third party cannot reconstruct the registration uniquely from these tabular data.

The rigid transformation maps an original measured vertex x_ij of component j to

An automated coarse-to-fine registration procedure was developed to minimize operator-dependent errors and to accommodate large initial translations, substantial relative rotations, and geometrically repetitive features. First, the centroids of the μCT and CAD surfaces were aligned. To avoid bias caused by non-uniform STL mesh density, the centroids were estimated from STL surface points sampled randomly but uniformly with respect to triangle area rather than from a simple arithmetic mean of the STL mesh vertices, so that they are not biased by non-uniform mesh density.

Principal component analysis (PCA) was subsequently used to determine the principal axes of both surfaces. Because PCA does not uniquely determine the signs and, in geometrically ambiguous cases, the correspondence of the principal axes, 24 proper three-dimensional orientation hypotheses satisfying det(R) = 1 were considered.

To resolve ambiguous orientations, feature matching using Fast Point Feature Histograms (FPFH) combined with random sample consensus (RANSAC) was used to achieve a robust global alignment. This step was particularly important for component D2, exhibiting relatively large μCT-to-CAD geometric differences, and for component D3, for which a substantial initial relative rotation and symmetry had to be resolved.

The resulting transformation was refined using robust multi-scale point-to-plane iterative closest point (ICP) registration with Tukey loss function to reduce the influence of local geometric discrepancies and outliers. A numerical target distance of 0.005 mm was used as a stringent ICP convergence criterion. It was selected to keep numerical contributions associated with the iterative registration, including finite convergence and rounding effects, small relative to the spatial sampling scale of the μCT data. The criterion itself is not an estimate of registration or measurement uncertainty.

For each measured vertex $x_{ij}^{*}$ of the registered μCT surface, the software determines the closest point c_ij on the corresponding CAD reference surface S_CAD,j,

The corresponding signed nearest-surface deviation (signed nearest-surface Euclidean distance) is defined as

The sign s_ij follows the inside–outside convention to interpret d_ij > 0 as the reconstructed micro-CT surface lying outside the corresponding oriented CAD surface, whereas d_ij < 0 denotes an inward deviation. This interpretation presumes consistent, correctly oriented surfaces. Near edges, openings, or geometrically ambiguous regions, the closest point may not define a unique functional feature deviation. Equation (3) is a signed nearest-surface distance, not necessarily the projection of a displacement onto a CAD normal, a radial deviation, or a fitted-dimension error.

The registration analysis presented in this report uses three component-specific exports. For each component, we retain the signed nearest-surface distance in the d_mm column and the measured-vertex coordinates without filtering, and convert distances from mm to μm. The resulting numbers of vertices are N₁ = 75 812, N₂ = 46 954, and N₃ = 16 416. Each vertex has equal weight; because triangulation density may vary across a surface, these distributions are not area-weighted estimates. The separately generated area-sampled normal deviations, ray-based radial differences, and repeated-registration distances represent different data products and are not pooled with the three vertex exports.

C. Empirical distributions, summaries, and KS comparison

For component j, the empirical cumulative distribution functions (ECDFs) of the signed and absolute deviations are

We report the mean d̅_j, median, standard deviation with divisor N_j − 1, root mean square (RMS), mean absolute deviation from CAD (MAE), and empirical quantiles of |d|. Specifically,

The percentile calculation uses linear interpolation between the sorted observations at position 1 + (N_j − 1)p. All reported quantities refer to the exported finite vertex sets; standard deviation across vertices is not μCT measurement uncertainty.

The signed histograms use common 20 μm bins from −500 to 900 μm, encompassing all exported deviations. For bin B_k of width h = 20 μm, its density is

Thus ∑_k hf̑_jk = 1 for each component. The diagnostic fraction of vertices inside a symmetric band is C_j(Δ) = 100G̑_j(Δ) in percent for half-widths Δ = 20, 50, 100, 200, 300 μm. These are shared comparison thresholds, not verified drawing tolerances or conformity probabilities.

For each component pair (j, k), the two-sided two-sample Kolmogorov–Smirnov (KS) distance is

D. Normality diagnostics and exploratory parametric fits

For each full signed vertex-deviation export, we calculate descriptive skewness and excess kurtosis and the Anderson–Darling (AD) statistic for a normal distribution fitted to that component’s deviations. The AD statistic is evaluated using logarithmic lower- and upper-tail probabilities to avoid numerical saturation. We also fit a normal distribution (two parameters) and a location-scale Student’s t distribution (three parameters) by vertexwise maximum likelihood and compare their Akaike information criteria, AIC = 2k − 2ℓ, where k is the number of fitted parameters and ℓ the maximized log-likelihood. These are exploratory fits to the marginal distribution of vertices on each particular mesh. Spatial dependence and nonuniform vertex density preclude the interpretation of conventional AD p-values or likelihood-based criteria as independent-part inference; the fits do not capture CT measurement uncertainty.

E. Detailed axial analysis of D1

The coordinates accompanying the D1 vertex deviations identify the y direction as its axial coordinate. To localize discrepancies without changing the observation set, we divide the observed y range into 12 adjacent equal-width intervals (the last includes its upper endpoint). Each bin contains the original vertices in that interval, with no resampling. We calculate its count, signed mean, RMS, 95th percentile of |d|, fraction with |d| ≤ 100 μm, and number with |d| > 200 μm. Count-weighted bin means and squared RMS values must reconstruct the global mean and squared RMS. D1 skewness and excess kurtosis are computed as descriptive moments of its full signed vertex distribution. Bins overlap in neither the vertices nor the coordinate range; their statistics remain sensitive to regional mesh density.

F. D1 CAD-normal projection sampling and radial diagnostics

A separate D1 export contains signed CAD-normal projection differences at approximately 80 000 area-uniformly sampled surface locations. We calculated its descriptive statistics directly from the point-wise data. Both the sampling locations and the deviation definition differ from the vertex-based nearest-point distances in (3). The signed CAD-normal projection deviation is defined as

where x_meas,i is the sampled point on the registered μCT surface, x_CAD,i is the corresponding point on the CAD reference surface, and n_CAD,i is the local unit normal of the CAD surface. We therefore use these data only for a D1 sensitivity comparison in the Discussion; they are excluded from the three-component ECDFs, histograms, KS comparisons, and axial analysis. The uniformity of the area sampling was not independently assessed in the present analysis; it can be checked against the available mesh and sampling procedure in a follow-up study.

The separate rotational-region export groups D1 along its recorded rotation axis into 497 adjoining intervals, each approximately 0.0244 mm wide. Its upper rendering assigns colors to 75 812 measured STL vertices according to interval; its lower profiles summarize radial differences at 79 841 valid locations from the separate area-sampled analysis. For each such location, the supplied radial quantity is Δr = r_CT − r_CAD relative to the recorded axis, and regional mean Δr, mean absolute |Δr|, and 95th percentile of |Δr| are shown. The 159 locations without a valid radial result are omitted. Where a ray intersects a threaded surface more than once, the correspondence and chosen radial intersection affect this diagnostic.

A further, separate ray-casting export samples 360 directions at 1° intervals in each of the 242 axial sections, spaced 0.05 mm apart. It retains 86 902 valid paired outer radial intersections of the CAD and aligned CT meshes. Its axis definition and ray correspondence differ from the preceding export; no additional registration was performed for this calculation. These radial outputs are exploratory measurements of radial correspondence, not estimates of the same signed distance or repeated-scan uncertainty.

A. Global vertex-based comparison

Table 1 summarizes the vertex-deviation exports. D2 has the largest RMS (207.72 μm), MAE (114.17 μm), and 95th percentile of |d| (526.78 μm). D1 has the smallest values of these three measures. The D2/D1 RMS ratio is 4.23, and D3/D1 is 1.57. These compare different component geometries with different mesh densities; they are not estimates of production capability. The means are positive across all three exports, but their magnitudes cannot be interpreted as calibrated material excess without verification of the CAD and sign conventions.

The ECDFs in Fig. 5 compare both the direction-sensitive distribution and the magnitudes |d|. The histogram densities in Fig. 6 use the same bin boundaries and vertical density unit for each component. D2 shows a broader distribution and a long positive tail. The two views agree on this descriptive contrast: histogram peaks depend on the bin width, whereas the empirical cumulative fractions do not.

Fig. 5.

Fig. 6.

Table 2 gives 100G̑_j(Δ) at selected diagnostic widths. At 20 μm, D3 covers 44.17% of its vertices, more than D1 (38.60%). D1 has higher coverage at 50 μm and every wider reported band. D2 has the lowest coverage at each threshold. This change in ordering between D1 and D3 illustrates why central concentration and upper tails should be shown together. These are fractions of mesh vertices; they are not area fractions or toleranced pass rates.

B. Normality diagnostics and parametric approximation

Table 3 reports additional diagnostics for the same complete signed vertex exports as Table 1. All three distributions are positively skewed and have positive excess kurtosis, most prominently D1 in the latter measure. The fitted-normal AD statistics are 795.21, 4414.05, and 833.92 for D1, D2, and D3, respectively. The conventional independent-observation AD calculation returns p ≤ 0.0005 (its lower reporting bound) for each component. This nominal result is consistent with the visible departures from a Gaussian shape, but its usual inferential calibration does not apply to spatially correlated vertices.

The location-scale Student’s t fit has a lower AIC than the normal fit for each component. The improvement is descriptive evidence that allowing heavier symmetric tails better represents these vertexwise marginals than a normal density. It does not establish that a Student’s t model adequately reproduces their positive skewness or spatially localized extremes, nor does it validate draws from that fit for uncertainty propagation.

C. Pairwise empirical-distribution comparisons

Table 4 presents the maximum vertical differences between the component ECDFs. For the signed deviations, the KS distance is 0.2474 for D1–D2, 0.1111 for D1–D3, and 0.2092 for D2–D3. Comparing absolute deviations gives 0.2908, 0.1175, and 0.2878, respectively. The signed and absolute comparisons answer different questions: a positive and negative departure of the same magnitude contribute differently to F̑ but identically to G̑.

Table 4 contains the observed distances only. Their magnitudes indicate how far the empirical vertex distributions separate on the displayed scale; they have no sampling-based significance interpretation for a population of parts. The three different components are each represented by one reconstructed mesh, so vertices cannot be treated as independent production replicates.

D. Spatial structure within D1

The 12 D1 axial intervals contain exactly its 75 812 exported vertices. Count-weighted pooling reconstructs the global mean and RMS in Table 1, confirming consistency between the global and regional computations. Their common width is 1.013 mm over the observed y range −0.036 to 12.119 mm. Fig. 7 shows that the largest bin RMS, 81.18 μm, occurs in interval 6, centered at y = 5.535 mm. Its Q_0.95(|d|) is 152.08 μm; the interval contains 6 307 vertices (8.32% of D1).

All 267 D1 vertices with |d| > 200 μm belong to interval 6, including the global extreme. Within a narrower ±100 μm diagnostic band, interval 4 has the lowest coverage (85.22%), whereas interval 6 has the highest (93.82%). The location that requires closer inspection, therefore, depends on the size of the discrepancy of interest. The signed D1 distribution has a descriptive skewness of 1.34 and an excess kurtosis of 12.05. These moments summarize its asymmetric tail and localized extremes; they are not estimates of an underlying independent measurement-error distribution.

Fig. 7.

Fig. 8 presents an axis-resolved view of D1, where the colored bands indicate uniform intervals along the rotational axis s rather than deviation magnitude. The lower panel shows the regional mean signed radial difference, the radial MAE, and the 95th percentile of the absolute radial difference. Across its 79 841 valid radial values, the supplied statistics give a mean Δr = 27.17 μm, an RMS of 80.49 μm, and Q_0.95(|Δr|) = 165.86 μm. Local 95th percentiles can exceed 0.6 mm at around s = −5.25 mm and −3.64 mm; these are sharply localized.

Fig. 8.