Application of an edge detection algorithm for surface determination in industrial X-ray computed tomography

Surface determination is an essential step of the measurement process in industrial X-ray computed tomography (XCT). The starting point of the surface determination process step is a single grey value threshold within a voxel volume in conventional surface determination methods. However, this value is not always found in the reconstructed volume in the local environment of the surface of the measurement object due to various artefacts, so that none or incorrect surfaces are determined. In order to find surfaces independently of a single grey value, a three-dimensional approach of the initial contour determination based on a Prewitt edge detection algorithm is presented in this work. This method is applied to different test specimens and specimen compositions which, due to their material or material constellation, their geometric properties with regard to surfaces and interfaces as well as their calibrated size and length dimensions, embody relevant properties in the examination of joining connections. It is shown that by using the surface determination method in the measurement process, both a higher metrological structure resolution and interface structure resolution can be achieved. Surface artefacts can be reduced by the application and it is also an approach to improved surface finding for the multi-material components that are challenging for XCT.


Introduction
Industrial X-ray computed tomography (XCT) is a method for non-destructive testing. With this procedure, components can be qualitatively inspected for internal defects such as cracks, cavities or pores. Industrial metrology also aims to gain quantitative results using such measurement systems. In sub-project C05 of the Collaborative Research Centre 285, the metrological applicability of XCT for use in joint inspection is being investigated [1]. Joining connections are a challenge for the XCT process due to several of their properties (Fig. 1). That is on one hand due to their material or material constellation, on the other hand due to their geometry and the associated location of the measurement points relevant for quality evaluation with interface and surface points. Strongly X-ray absorbing metal leads to metal artefacts and different material combinations lead to multi-material artefacts due to an indifferent grey value assignability. Both can have a negative effect on the surface identification. In the same way, absorbent materials require higher X-ray power, which is associated with a reduction in the measurability of smalldimensional features [2]. Irregularly shaped geometries as well as the component orientation in the measuring system can lead to insufficient penetration of X-rays and to a diversity of grey values and thus to the different assessment of volume components despite the same material properties and can also have a negative effect on surface identification [2]. The grey values result from the fusion of several projections of an energy integrating X-ray detector, based on the absorption of X-rays, to a volume model. However, the surface determination of a component is an essential step in XCT, because the surface represents the geometrical properties of the measurement object. Standard surface determination methods use single grey values or grey value intervals for surface determination or they use at least single grey values as starting value for surface determination based on a local grey value evaluation. If the used grey value is not in the immediate surroundings of a region of interest, the image or volume information contained in it is not taken into account from the start, which is why this starting contour in particular is important for surface determination. Due to the different interpretability of grey values and the proven improvement of the surface determination with local grey value consideration, it is obvious to carry out the entire surface determination procedure dynamically and thus also independently of static grey values for the starting contour. In this work, therefore, the starting contour is determined by using a known edge detection filter applied to the sectional images of the reconstructed volume dataset in three orthogonal directions beforehand. Because the edge detection filter is included as additional information on a voxel basis in the surface determination on a sub-voxel basis, the risk of information loss due to the filter application is much lower in contrast to the filtering which is performed directly and repeated on the reconstructed volume. Using different test specimens with their individual properties, which are also relevant in the long term for the assessment of joint connections, the method of surface determination is examined for its performance both on the basis of real part as well as simulated measurements. The results of the different test specimen measurements in combination with the modified surface determination show that the methodology used is not limited to joint connections, but can also be used for surface determination in general.

State of the art
The state of the art covers two topics. First, the concept of XCT for dimensional measurement with different kinds of surface determination and often occurring artefacts. The second focus is on the filter algorithms which are used and known in conventional image processing.

X-ray computed tomography
Industrial X-ray computed tomography provides volumetric imaging for non-destructive testing. In XCT a component is positioned in the measuring device between an X-ray tube and a detector (Fig. 2). After scanning the part with a high number of projections regularily distributed angular positions with respect to the rotation axis of the manipulator, a volume model can be reconstructed, consisting of volumes of pixels (voxels). To each voxel a grey value is assigned to. Then the models of the components can be examined for internal gaps, cracks and pores [3]. For dimensional metrology it is important to extract a surface from the volume on the basis of which measurements can be made with regard to the component characteristics. Conventional evaluation software provides surface detection methods that define surfaces based on grey values.

Methods of surface determination
A standard method of surface determination is the ISO50method [4]. It is the use of a single threshold grey value applied globally to the entire volume with possible grey values between 0 and 65.535. The boundary surfaces which mean the border between the material and the  background lies in the middle of the two grey value peaks in the histogram (Fig. 3). A more advanced method is the use of this grey value as a starting value for local adaptive thresholding (LA) [4]. This means the different grey values are used as threshold values depending on the surrounding voxel grey values. The number of voxels to be included can be varied. If the component consists of several materials that also differ in terms of their absorption, a further material peak appears in the grey value histogram (Fig. 4). To find the surface of each material, a further threshold is then defined between the grey value peaks. These can also be automatically set as local grey value peaks within an interval. Often, the histogram distributions overlap or a material peak is so small that it disappears within the grey-scale range, so that their separation is not trivial. If a component is standardised or known, a virtual surface model like STL can be used instead of a grey value as the starting contour for LA surface determination. The procedure is particularly suitable for comparing components that remain the same.
Another approach to determine surfaces is the search for contours. Well-known edge detection algorithms in image processing are for example Canny, Sobel, and Prewitt [5]. Due to its background, this is a two-dimensional procedure. Yagüe-Fabra et al. describe in [6] a 3D edge detection technique for surface extraction in computed tomography for dimensional metrology applications based on the Canny algorithm with sub-voxel resolution. Ontiveros et al. compare the 3D Canny and the local adaptive method with the result of a higher accuracy with regard to the local adaptive method while the Canny algorithm is attributed a high repeatability [7]. When evaluating interface structures of a component made of two different polymers, however, the edge detection methods are credited with the more accurate results by Jiménez-Pacheco et al. [8].
In [9] a 3D adapted edge detection algorithm based on Deriche [10] is used. After the voxel-based edge detection, the sub-voxel refinement is performed using the gradient information of the neighborhood voxels. Due to a significant reduction of the deviations obtained with the Derichebased algorithm compared to the Canny-based algorithm, the Deriche-based algorithm is attributed the potential prove pivotal in reducing systematic errors and uncertainty in CT measurements. The Deriche-based algorithm also showed improved-though not satisfactory-results for the challenging multi-material components [9].

Artefacts
XCT is as shown in numerous standards on computed tomography an artefact afflicted procedure [11]. Many of those artefacts are unavoidable in multi-material setups due to their different absorption characteristics [12]. After the surface of a component has been defined, artefacts are often part of this surface. A reason for this is the energydependent non-linear interaction. Lower energy photons are attenuated more than compared to higher energy photons, causing the average energy of the beam to increase while passing through the component. This beam hardening leads to different grey values in the reconstruction even if the same material properties such as electron and mass density are present in the component. If the most commonly applied Feldkamp algorithm is used for volume reconstruction and there is a parallel surface in the cone beam CT that is horizontal from the central beam path, its formation is influenced by the cone beam or feldkamp artefacts. Feldkamp or beam hardening artefacts are also mentioned in [2] as physical causes of structural resolution loss.

Filter application in image processing
An approach to determine quality parameters specifically for joints using micrographs based on a filter was recently Illustration of a grey value histograms with a second threshold value due to a further material influence for defining a surface or a start value for surface determination presented by Zirngibl et al. [13]. Applying mathematical filters is also a commonly used approach to reduce artefacts in CT. For example, an improved Sobel filter is used in [14] for better detection of internal crack defects in slice images. In [15] an improved segmentation method by an enhanced gradient norm of the contact zone of two plastic parts is presented. Araújo et al. used a combination of the smoothing anisotropic diffusion and edge enhancement unsharp mask filter for enhancement of micro CT images of steel cracks [16]. The use of filters carries the risk of destruction of information [17]. For example, the consequence of applying a Gaussian or Median filter to a dataset can generate a part which appearance is smoothed. This may be desirable in the case of reverse engineering. But in case of a component inspection the real surface and the examination of internal component changes are of interest. Nevertheless image smoothing is also common practice at Deriche [10] and in general a step before edge detection. Since this paper is mainly concerned with the Prewitt filter as a edge detection algorithm, it is presented in the following. The edge detection algorithm named after Judith M. S. Prewitt calculates two gradient images G x (Eq. 1) and G y (Eq. 2) in horizontal and vertical direction by convolution of a greyscale image A with a small, divisible and integer filter [18]: Compared to the Canny or Sobel algorithm, this structure leads in a high degree of sensitivity to horizontal and vertical edges [19]. The results are used to calculate the edge magnitude G (Eq. 3) and the edge direction (Eq. 4).

Methods
For the generation of the volume models on which the surface determination was applied, the CT measuring system Zeiss Metrotom 1500 and the simulation tool aRTist of the German Federal Institute for Materials Research and Testing (BAM) were used. Five test specimens or combinations of test specimens were used (Table 1). Each test specimen has properties that are also relevant when examining a joint. With the help of the hourglass specimen (Table 1a), it is possible to examine which dimensional measurement value the XCT assigns to the smallest geometrical properties of the test specimen. Another test specimen is a gauge block specimen (Table 1b) consisting of gauge blocks [20] made by steel which build a defined gap, like it is used in [21]. This can be used to practically investigate which distances between two components are recognised and which measured values are assigned to these distances. The next test specimen consists of two stacked universal test specimens (Table 1c) with steel at the top and aluminium at the bottom. It was used in this work to analyse the segmentability of components consisting of several materials. A STL file, a cylinder with bore holes (Table 1d) with diameters between 4 μ m and 40 μ m, was used for simulation purpose in aRTist. Similar to the gauge block specimen, this test specimen enables the investigation which distances or bore diameters can still be recognised. Last but not least a calotte cube developed by Physikalisch-Technische Bundesanstalt Braunschweig (PTB), [22] (Table 1e), a cube with spherical calottes, which sizes and distances are calibrated, is used to make a general statement about the trustworthiness of length and size measurements with regard to the compared surface determination methods. Further details on the test specimens as well as on the measuring system used and its process settings are given in Table 1. The experimental procedure for filter-based surface determination of the generated measurement volumes is shown in Fig. 5.
In a first step, the reconstructed volume is loaded into the measuring software VGSTUDIO MAX (VGS) (Version 3.5 from Volume Graphics). The second step is the export in the form of a slice image stack. In this case, the file format DICOM (Digital Imaging and Communications in Medicine) was used because it contains the parameters for the reconstruction. The image stack is imported into the numeric computing environment MATLAB (The MathWorks, Inc.) where it forms a 3D cubic voxel volume. Within the third step, the image stacks were processed with the filter functions (Sect. 4.1). Canny, Prewitt and Median were tested in different sequences on the universal test specimens stack coboid. In the case of using the Median filter, a size of 2 × 2 voxels was used.
Matlab provides the standard algorithms with the functions: PrewittImage = edge(Image,'Prewitt'), CannyImage = edge(Image,'Canny') and MedianImage = medfilt2(Image, [2 2]). Starting from the standard 2D image processing with a single stack of images, the filter operators are also applied to the other two image stacks resulting from the orthogonal Table 1 Test specimens used, measuring device, specifications and settings as well as steps and settings for surface determination in the software directions, which is understood in the following as a kind of 3D application. Each information about edges >0 (0 = black/no information) obtained from the filter applications was assigned a uniform grey value which are collected in a separate 3D volume according to their position-related information origin. This results were visually evaluated by a mean layered image with regard to image artefacts and edge detectability. Considering the two characteristics, only the sequential application of Prewitt and Median filtering was chosen, which was also used for the surface finding of all other test specimens (Sect. 4.2 et seqq.). All determined information was saved in a new DICOM-file with the already existing reconstruction data. In the fourth step, the new image stack was imported into VGS. The ISO50 method was used to form a starting contour from the edge information.
Afterwards in a fifth step the prior image stack was imported too, which ensures an exact local registration which means an exact overlay of the two data sets. The volume that is determined from the voxel edge dataset after surface definition using the ISO50 method was used as the starting contour for LA surface determination based on the prior image stack. Starting from this contour search beam lengths were varied according to Table 1. To reduce or completely eliminate the many small edge detection filter-based artefacts that arise after the surface was determined, which are illustrated exemplary in Fig. 5, in an optional final step a region of interest (ROI) was generated from the resulting volume. This ROI was minimally reduced at its border area by the VGS function erosion. Search beam extension when using the LA method is necessary by using higher erosion values. By splitting the ROI, it is possible to keep only the largest components in the measurement volume. The reduced ROI was then used as the starting contour for another LA surface determination.
For the comparison of the edge-based surface determination method with other conventional methods, the prior image stack was again imported into VGS. For the monomaterial test specimens, the LA surface determination method was used. For the universal test specimen stack coboid consisting of two materials, the multi-material (MM) method and the LA based on a STL start contour were used to determine the surface. For comparability of the edge-based surface determination with other methods, the final used search beam length of the edge-based surface determination always corresponds to the length of the comparison method ( Table 1). The local surface search by using defined search beam length has a more important role in the comparison of these surface finding methods, because we have different start values and surface is always searched and defined in a defined environment.

Results of the application of different filter operators
The results of the filter applications for a mean layered image of the universal test specimens stack coboid are shown in Fig. 6. The Canny algorithm (Fig. 6a) shows strong image artefacts compared to the Prewitt filter (Fig. 6b), but the edges in the area between the two test body elements are better represented. With filtering applied in all three orthogonal directions the filter artefacts increase (Compare Fig. 6c, d). When using the Prewitt filter, there is a slight improvement in edge detection between the two test piece components (Fig. 6d). By applying a Median filter before edge detection (Fig. 6e, f), an artefact reduction is shown, especially with the Prewitt filter (Compare Fig. 6d, f), but the edges between the components are harder to detect. By using the Median filter after the edge detection (Fig. 6g, h), a slight artefact reduction can be seen, especially with the Prewitt filter (Compare Fig. 6d, h), compared to the 3D application without Median filtering (Fig. 6d).
In the following only the 3D-Prewitt filter with subsequent median filtering is used as a start contour and will be referred to as 3DPM in short. A reason for this is the artefacts in the Canny algorithm increase the number of iterations required for edge detection. Furthermore, edges are recognised poorly overall in a non 3D filter application, which shows itself in a holey surface after the surface determination. Median filtering does not seem appropriate before edge detection due to a possible loss of edge information. However, it is used after edge detection to close the edge contour between the voxels and for a slight reduction of artefacts as well as a resulting reduction of computing time.

Structural resolution according to the hourglass specimen
The Houglass specimen is a test specimen which, according to Zanini et al., can be used to determine the metrological structural resolution (MSR) of a CT [23]. Two spheres stacked in a carbon tube form a theoretically pointshaped transition at their contact surface. In the reconstructed volume model, a specific neck thickness is formed depending on the selected device parameters. Depending on the ball diameter D used, the characteristic parameter h, which provides information about the still detected interface length, can be concluded via the diameter of the neck thickness d (Fig. 7) according to Eq. 5. Table 2 shows the diameters d of the extended contact point of the theoretically point-shaped transition determined via an fitted Gaussian circle as a dependency of the  search beam length and the surface determination method LA and 3DPM used. Overall, smaller diameters and thus smaller values of h can always be determined with the 3DPM method. Both methods show deviations of h of about 10 μ m with variation of the search beam length. With the LA method, only slight improvements in the value of h can be achieved with a search beam length of more than 10 voxel . Shorter search beam lengths lead to a lower metrological resolution. The best resolution is achieved with a search beam length of 4 voxel when the 3DPM method is applied, which also indicates that the starting contour is closer to the true value. Figure 8 shows that even smaller diameters d can be determined above and below the contact point of the balls using the 3DPM method.

Surface artefacts and structural resolution on the gauge block specimen
The gauge block specimen consisting of several gauge blocks as used in [21] is another approach to the study of structural resolution. The focus is on the smallest structure that can still be determined between two surfaces and is called interface structural resolution (ISR). But as described in [24] and as shown in [21], steel gauge blocks are difficult to measure because of the artefacts caused by the effect of scattered radiation and beam hardening in the volume data [25]. The result of a standard surface determination with the LA method with a search beam of 4 voxel at the 50 μ m gap gauge block specimen is a data model which shows strong artefacts, both at the outer regions and the transitions of the individual specimen components (Fig. 9a) as well as at the two gap-defining surfaces (Fig. 10a). The reason for the external artefacts is the distance of the starting contour of the LA method, based on the ISO50 value, from the real surface. The artefacts between the gap are mainly due to the fact that the ISO50 value is only partially represented in the gap, which leads to a tapering of the measurement surface and a lower informative value with regard to the ISR [21]. Extending the search beam to 15 voxels results in a minimal expansion of the gapdefining surfaces (Fig. 10b), but most significantly, outer surface artefacts are reduced (Fig. 9b). Using the 3DPM surface determination method in combination with erosion and also with search beam of 15 voxels, it shows a significant reduction of the described artefacts (Fig. 9c) as well as a higher resolution with regard to the ISR (Fig. 10c). When examining the gauge block specimen with different gap sizes, a gap of 40 μ m no longer be detected using the LA method. The 3DPM method still detects a gap-size of 30 μm.

Interface and surface detection on the multi-material specimen
In the case of the multi-material universal test specimens stack coboid, the aluminium body, which absorbs less radiation than the steel, is a challenge in regard to surface determination in the region of transition to the steel body. The use of the LA method with a single ISO-threshold as the starting grey value makes no sense with the present constellation of two materials that differ strongly in their absorption properties. For a qualitative comparison of nominal and real values with regard to surface recognitions on the aluminium component, the LA method using an STL start contour, the MM surface recognition (VGS) as well as the 3DPM method were used (Fig. 11).
For the one-time registration of the nominal STLobject, which was also used as the starting contour for the method LA, it was fitted to the surface of the volume data defined with the method MM on the basis of two lateral and the lower surface which are far away from the steel and therefore little influenced by it. The surface of the nominal object was used as a reference for the nominalactual comparison and the qualitative values in Fig. 11 are displayed on it.
When comparing the surfaces of the aluminium part, it is noticeable that the MM method (Fig. 11a) reproduces the surfaces less well than the 3DPM method (Fig. 11b) and the LA variant based on an STL (Fig. 11c). The reason for the comparatively better detection using the LA method is that the registered nominal object corresponds to the start contour and thus represents the supposedly best possible case. The comparison of the LA methods with STL (Fig. 11c) and 3DPM (Fig. 11b) shows similar results. In some places, such as the cylindrical structure, the surface is shown in a more recognisable way. Even in the visual comparison of the steel-air-aluminium transition the MM method shows insufficient results to the other methods (Fig. 12), because in the border area between two differently strong radiation absorbing components often no interface is detected, but a single material transition between the high absorbing and low absorbing materials was determined (Fig. 12a). The result is a good visibility of the highly absorbent material (Fig. reffig10a). The volume that lies between the steel and aluminium components is added to the weaker absorbing material and consequently the aluminium component surface is not recognised or is recognised less well. Because the application of the method LA with the start value of the ISO50 method is not reasonable and the method MM almost does not provide a correct surface determination for the aluminium component localised on the steel, and a correct STL as start contour for the method LA is usually not available (Fig. 12c), the use of the method 3DPM (Fig. 12b) is a reasonable alternative procedure to be used.

Structural resolution on the virtual boreholes
In the case of the bore specimen, bores up to 30 μ m are detected when determining the surface according to the LA method. Surfaces on smaller holes are not detected at all. With the method 3DPM, bores up to 30 μ m are also detected very well (Table 3). But even up to 20 μ m, holes are still partially separable. In the case of smaller holes, only the existence of inner structures is recognisable without being able to specify them more precisely.

Geometric deviations on the PTB calotte cube
The cube is a specimen whose spherical radii and spherical centre distances are calibrated. The cube is therefore suitable for the evaluative quantitative comparison of the surface finding methods. The determined radii of 33 spherical calottes were compared as well as 10 spherical calotte distances between 1.600 and 9.049 mm. With the exception of the 50 spot/voxel size, the deviations for the individual radii as well as for the distances were always in the single-millimetre range and increase with increasing voxel size. Figure 13 shows the average amount of deviation per radius measurement of the 33 calibrated radii depending on the surface determination method. For the LA method, the deviations increase almost linearly with the voxel size. In this case, the 3DPM method always Fig. 12 Comparison of surface contours on the basis of sectional images of the stacked universal multi-material specimen using surface determination: a MM, b 3DPM and c LA with STL starting contour Table 3 Borehole detection as a dependence of the borehole size (in μ m) on the bore specimen using the method 3DPM Fig. 13 Comparison of the average absolute value of deviation of 33 calibrated spherical calotte radii shows lower deviations in comparison. The centre-tocentre distances of the spherical calottes show the opposite picture (Fig. 14). The deviation totals for the 3DPM method are higher and thus comparatively slightly worse, whereby identical values were determined at a voxel size of 20 μm.

Conclusion
In this paper, a method for surface determination using a 3D edge detection approach was presented. It was shown on a multi-material test specimen that the Prewitt edge detection has less image artefacts compared to the Canny algorithm. Using different test specimens, it was possible to show that smaller surface structures as well as interface surface structures were better determined compared to conventional surface determination. Artefacts on the component surface and between two surfaces could be reduced. An improvement in the detectability of intermediate surfaces and surfaces of a weaker absorbing material can also be seen in the case of challenging multi-material components. When comparing measurements of calibrated radii of the same size, better measurement results were also achieved, even though there is minimal deterioration in the measurement of calibrated spherical calottes distances. Due to the advantages of the 3DPM method mentioned above and its relevance for joining connections, it is particularly suitable for the computed tomographic investigation of these joints created in versatile process chains.