COMPUTED TOMOGRAPHY DATA COLOURING BASED ON PHOTOGRAMMETRIC IMAGES

Nowadays various methods and sensors are available for 3D reconstruction tasks; however, it is still necessary to integrate advantages of different technologies for optimizing the quality 3D models. Computed tomography (CT) is an imaging technique which takes a large number of radiographic measurements from different angles, in order to generate slices of the object, however, without colour information. The aim of this study is to put forward a framework to extract colour information from photogrammetric images for corresponding Computed Tomography (CT) surface data with high precision. The 3D models of the same object from CT and photogrammetry methods are generated respectively, and a transformation matrix is determined to align the extracted CT surface to the photogrammetric point cloud through a coarse-to-fine registration process. The estimated pose information of images to the photogrammetric point clouds, which can be obtained from the standard image alignment procedure, also applies to the aligned CT surface data. For each camera pose, a depth image of CT data is calculated by projecting all the CT points to the image plane. The depth image is in principle should agree with the corresponding photogrammetric image. The points, which cannot be seen from the pose, but are also projected on the depth image, are excluded from the colouring process. This is realized by comparing the range values of neighbouring pixels and finding the corresponding 3D points with larger range values. The same procedure is implemented for all the image poses to obtain the coloured CT surface. Thus, by using photogrammetric images, we achieve a coloured CT dataset with high precision, which combines the advantages from both methods. Rather than simply stitching different data, we deep-dive into the photogrammetric 3D reconstruction process and optimize the CT data with colour information. This process can also provide an initial route and more options for other data fusion processes.


INTRODUCTION
3D reconstruction is the procedure of generating 3D appearance and structure of real objects in close range. Based on different principles, various methods can be utilized, such as frequently used monocular cameras, binocular cameras, video cameras and laser scanners (Ma, Liu, 2018). In addition to the surface reconstruction of the real object, 3D imaging methods originating from the medical field such as CT or MRI are also used to obtain the internal structure of the object. Different methods have their unique advantages under certain application conditions. However, there is still no single feasible way to bring all advantages together. Therefore, data fusion is the way to provide optimal quality information through combining the sensitivity and specificity of different sensors.
The object involved in this research and presented below is a classical gyroscope with a detailed and complex structure. Overall, it is part of a German pilot project to digitize an existing collection of gyroscopic instruments in 3D Virtual Reality models (Fritsch et al., 2018). Therefore, the digitization task relies on methods with high precision options. Photogrammetric 3D reconstruction with appropriately collected high resolution images is able to reach high quality textured surface models. However, sometimes the results might not be sufficient if the viewing angles of the cameras are limited, or homogeneous surfaces without enough features are being processed, or the illumination situations are unsuitable. In our application, * Corresponding author endoscopes are introduced as a supplementary method for local reconstructions. Besides the photogrammetric application with various sensors, CT is another method to complement the data for the internal and invisible part of the object.
CT is used to reconstruct internal structures of an object that is not visible from outside. It is an imaging technique taking a large number of radiographic measurements from different angles, in order to generate slices of the object. All reconstructed slices are then integrated into a uniform 3D coordinate system to construct the complete 3D volumetric representation. CT has been widely used not only in the medical fields, but also in industrial (Kastner, Heinzl, 2015;Shen et al., 2004;Carmignato, 2012) and archaeological domains (Yoshikawa et al., 2003). However, the CT volumetric data lacks coloured texture.
Besides CT, photogrammetry with its main task of 3D reconstructions have been developed for decades. This mature technology utilizes multi-view images to generate a 3D surface model with texture information, which is preferable for human perception and interpretation. Despite some limitations, such as regions with less features or changing illumination (Gienko, Terry, 2014) conditions could lead to unsatisfying results in 3D reconstructions, the photorealistic texture could provide a perfect supplement for CT data.
Early research related to CT colour texturing are in the medical field for diagnostics purposes, which focus on the data segmentation. Kerr et al. (2000) created a CT colour lookup table  for visualization, and Park et al. (2005) proposed a semiautomatic way of colouring CT data using Adobe Photoshop ® . There is also work on mapping the texture from digitized colour portraits to the 3D CT model (Xia et al., 2000). Those approaches either assign artificial textures or are mapping the texture from 2D images to 3D CT data under strict experimental conditions. More recent work draws attention on photogrammetry, which aims to fuse CT 3D volumetric data and photogrammetric point clouds. Their main concern is the registration accuracy (Maal et al., 2008), (Khambay et al., 2002). In principle, the data fusion can provide more information for diagnostics in the medical field, but the data itself and its attribution do not have accuracy improvement at all, which is the main contribution of this work.
In the case presented below, the CT data has a voxel size of 64.5μm and present a high precision and completeness in comparison with the photogrammetric point clouds, but have no colour information. A direct way to assign colour values to CT is to find the correspondences with the photogrammetric data. However, poor quality point clouds from photogrammetry may not be able to provide sufficient and correct texture information for colouring CT data. Therefore, it is advisable to turn to the colour sourceimages for 3D reconstructions. The colour information from the image pixel could be assigned to CT voxels with the transformation relationships between CT data and photogrammetric images (Fritsch et al., 2018). As it is shown in the workflow chart ( Fig. 1), the first step is the generation of CT surface point clouds and photogrammetric point clouds of the object's surface. Secondly, we take the surface of CT data for accomplishing the alignment with the photogrammetric 2.5D surface model. Thirdly, the typical procedure for alignment between CT surface data and photogrammetric data is to select control points from both data sets to obtain a coarse-registered model and then an ICP (Iterative Closest Point) algorithm or its variant could be applied for refinement.

THEORY AND METHODOLOGY
Fourthly, after the alignment, the transformations between CT data and images for the 3D reconstruction are determined. According to each image pose, the CT surface data can generate a depth image (see Section 2.5), which corresponds with the original photogrammetric image. Fifthly, the exclusion of the ambiguous pixels on the depth image due to the projection of 3D points should be removed from the view in reality by making using of range values. And finally, all points from CT surface data are processed in the same way in order to obtain a full coloured CT surface representation.

Photogrammetric Point Cloud Generation
The photogrammetric 3D reconstruction utilizes multi-view images. After the data acquisition, the pose of the images is estimated by the principle of SfM (Structure-from-Motion). The result will be applied in dense image matching process for a dense point cloud of the object surface. The meshing and texturing step will make a more realistic model by unwrapping and calculating texture images colour. The workflow can be summarized in Figure 2.

Threshold Determination of CT Scan
The reconstructed CT volume represents the object with attenuation values carried by each voxel. The overall histogram of the final volumetric data illustrates the different peaks with respect to the various materials used in object, since each material attenuates in different way when the X-ray passes through. Theoretically, the overall histogram of the final CT volume calculated with the FBP (filtered backprojection) algorithm carries multiple peaks representing different materials, Fig. 3. However, the X-ray images have information not only about the accumulated attenuation on the detector but also the noise which can comes from the nature of monochromatic X-ray source causing beam hardening artefacts and scattering (Krumm et al., 2008). Based on the severity of the noise which can distort the shape of each peak one can separate materials with the simple thresholding (Otsu,1979). Additionally, some parts of scanned object can be simply not visible during the visualization as a result of high X-ray energy spectrum: Materials such as paper or plastic will have low attenuation in comparison with the metal and dense materials with high atomic number.

CT Point Cloud Surface Extraction
The CT volume can be converted to the point cloud according to the spatial coordinates of each voxel. The coordinate of each voxel is calculated with the use of the grid origin of the regular CT volume and grid spacing, which describes the initial position of the volume and space between two voxels respectively. However, in the case of the CT volume the colour information is missing due to the absence of an appropriate energy spectrum on detector side. Alternatively, the voxel values can be normalized and mapped to [0.255] range and stored as grayscale for a colour field of the point cloud.
Having the whole point cloud of the CT volume, the surface is necessary to be extracted for the colouring process. The solution is to define the threshold to the volume data to separate the background and the material and then characterize the problem as finding the convex hull which is the set of those points which are enclosed in the convex polygon. By projecting imaginary parallel rays to the surface of the volume, the outer point will be touched by the ray. Due to the possibility of noise and voids, which can happen in the thresholding selection process, a certain number of voxels after the first intersected voxel need to be checked to justify the integrity of the volume.

Photogrammetric and CT Data Registration
The registration of pair-wise point clouds is achieved normally by the iterative closest point (ICP) algorithm or some variants. But first of all, a coarse registration is necessary to provide an initial position for the source and the target point cloud. The coarse registration here is realized by picking manually control points from both data sets. In the next step, a 7-parameter transformation is used, which is expressed as, where and are coordinate vectors of control points from the source and target point clouds respectively, contains the translation parameter ( 0 , 0 , 0 ) , is the scale and is the (3 × 3) rotation matrix. This transformation can be reformulated as a linear model through Taylor expansion, = ( , , , , , , ) where denotes differential changes of translations, rotations and the scale.
When applying the well-known least squares Gauss-Helmert model for the determination of the unknown parameters x, its first order can be formulated as The second order or dispersion model D is as follows: with v the vector of inconsistencies, P the weight matrix of the observations and the variance of the unit weight.

Depth Image Generation and Application
A depth image is an image channel containing the distance information from a viewpoint to the scene objects. It can be obtained from multi-view images for the purpose of 3D point calculation. However, a depth image can also be generated from a 3D point cloud with adequate camera parameters for other applications. The process is a projection from the 3D space to the 2D image frame, where pixels carry not colour but range values to the corresponding 3D point.
If the camera parameters from the pose estimation and calibration are adopted, the depth image should also be overlapping with the original photogrammetric images in the projection areas. To get the conceived result, the intrinsic parameters, including the calibrated focal length and the principle point coordinates, and extrinsic parameters are to be provided to define the camera. In addition, the width and height of the depth image are needed to be set.

Colour Information Assignment
Because the range image generation process projects all points onto the image plane with the camera intrinsic parameters, some points on the back side, which should not be projected according to the law of light propagation, will fall on the image plane. This will cause ambiguities when finding the corresponding pixels between the depth image and the undistorted image. The ambiguity may happen as indicted in Figure 5. Another possibility for neighbouring pixels can also be that the light ray travels through two surfaces and have two intersection points. Therefore, the way to solve the problem is by taking advantages of the range values of neighbouring pixels. The procedure to eliminate the wrong pixels is as follows:  Create the depth image according to the camera pose.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition)  List the pixel coordinates and colour information according to the depth image and corresponding photogrammetric image for each 3D point.  Take a certain size window of the pixel, and compare the range values between the target pixel and its neighbouring pixels. The pixel which holds much bigger value ranges will be excluded from the colouring process.  The qualified pixel will get the colour information from the corresponding pixel of the undistorted photogrammetric images.
By applying this processing pipeline, the corresponding pixels of the back side of the object are eliminated. However, a single point is still visible from multiple images, which means one point from CT surface might get more than one colour values. Therefore, a further range comparison is necessary to determine the final colour value.

Experimental Object
This research is implemented under the project Gyrolog, which digitizes a huge collection of gyroscopes, maintained by the University of Stuttgart, Germany. The goal is to preserve the instruments virtually in 2D and 3D. Thus, for the 3D VR digitization we combine different 3D imaging technologies. The structures of the gyroscopic instruments are of high complexity and therefore challenging the 3D digitization tasks.
The experiment is conducted with the Honeywell Golden Gnat sensor, one of the objects from the Gyrolog project. It was a miniature rate gyro and mostly part of guidance and control systems during the cold-war era. It is a floated, single-axis gyroscope and weights circa 150-160gr. The Golden Gnat is driven by a hysteresis motor that turns the rotor (the gyro part of the gyroscope) with 24,000 rpm (rotations per minute) and senses angular velocities from around 0.01°/sec. to 50°/sec. Figure 6. The Golden Gnat of Gyrolog project

3D Imaging Sensors
Our 3D reconstruction is implemented with a full frame D-SLR camera, Sony α 7R III with Zeiss Loxia 2.4/25mm lens. The resolution is 7953 * 5304 pixels. Figure 7. Sony α 7RIII camera, Zeiss Loxia 25mm lens The object to be digitized is placed on a turntable under an appropriate lighting configuration. The camera is mounted on a tripod to take pictures. Normally the imaging distance is around 40-50 cm. By rotating the turn-table 360 degrees, the pictures from different views can be obtained. It's also necessary to adjust the height of the camera to make two or more circles of images. The point density is calculated to be about 90 μm theoretically from the intrinsic parameters of the camera and an average imaging distance. In reality, inappropriate illumination, textureless surface or poor image configurations will influence the precision of the object coordinates (Remondino et al., 2017). The setup can be seen in the Figure 8.

Figure 8. Photogrammetry setup
With an industrial CT scanner, electrons are accelerated in the X-Ray tube in the direction of a target. Due to the deceleration of the particles, a beam is generated within the X-Ray band. This beam can be used to penetrate an (also metallic) object and generate a projection on the X-Ray detector. Taking projections from different angles of the object under investigation the machine provides data to reconstruct a 3D model of the object. In principle, all the internal information of the object is stored as the intensity of the voxels. However, the pre-processing and extraction of the information can be challenging. The threshold extraction between high penetrability and low penetrability is easily done by visualizing the intensity histogram of the whole objects, but threshold extraction between the noise of the air and highly penetrating objects can be difficult. Therefore, our experiment is mainly concentrating on the metal surface rather than the plastics parts. Further research on how to properly extract the threshold of high penetrability material is undergoing.

Data acquisition process of the Golden Gnat:
With the setup introduced in section 3.2, the Golden Gnat is imaged with Lego bricks fixing its base. One of the original images is shown in Figure 10. The photogrammetric 3D reconstruction is implemented via RealityCapture software. To better focus on the object part for the experiment, the point cloud is filtered to be displayed in Figure 11. Regarding CT data, an example of X-ray image is displayed in Figure 12. Figure 12. An example of an X-ray image of the Golden Gnat According to Figure 3, the peaks of the histogram are very close and the contrast between materials is low. Picking up the threshold by applying the thresholding algorithm (Otsu, 1979) give unsatisfying result, less clear than expected as shown in Figure 13.

Figure 13. Threshold selection by the Otsu algorithm
Applying the visual inspection, the threshold value can be selected in order to improve the quality by balancing between the noise and object integrity as shown in Figure 14. Figure 13 and 14 are made in VgStudio (Volume Graphics, 2020).

Point Cloud Registration of the Golden Gnat:
The registration of CT and photogrammetric data of the Golden Gnat is divided into two steps. First of all, the control points from both datasets are manually selected. The control points should distribute evenly on both surface point clouds and should avoid The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) coplanar situations. Secondly, with an initial registration, an ICP algorithm is applied with PCL under C++, to calculate the closest distance between the overlap of the two point clouds. The aligned CT point clouds and the photogrammetric source point clouds are displayed in Figure 16.

Depth image generation of the Golden Gnat:
After the transformed CT surface data is acquired, the depth image of the CT point cloud needs to be calculated. The depth image generation process is realized under PCL from aligned CT point clouds. In this research, the needed camera intrinsic and extrinsic information is exported from the 3D reconstruction of the Reality Capture software. However, the definition of the parameters, especially the camera poses, which have multiple representation possibilities, are crucial for the implementation of the procedure.
The depth image corresponds to the pose of the original image of Figure 10 and is displayed in Figure 18. It should be noted here, that the result does not correspond to the original photogrammetric image but to the undistorted image, which needs to be calculated through the original image and applying camera distortion parameters. Figure 19 indicates the overlap between the depth image and the undistorted photogrammetric image. Figure 19.Overlap display of the depth image and the corresponding photogrammetric image 3.3.5 Colour information assignment for the Golden Gnat CT surface point cloud: By applying the method stated in section 2.6, a window size of 5*5pixels is chosen to exclude the ambiguous pixels for the colour assignment process. Figure 20 displays the final results of the coloured CT point clouds. It can be seen that except the wired part all is coloured with the information of its wrapping plastics. Thus, the CT point cloud is coloured with satisfying quality from photogrammetric images.

Figure 20. CT coloured point cloud of the Golden Gnat
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition)

Comparison of two data
A comparison of a zoomed local display of both the original photogrammetric point cloud and the coloured CT point cloud is shown in Figure 21 and 22. From Figure 21 and 22, it can be concluded that the high precision geometry of the CT data and texture advantage of photogrammetry data are integrated in the coloured CT point cloud, which is the aim of our research.

The Limitation of the Research
The experiment is implemented with some assumptions, which can also be regarded as limitations of the methods. First of all, the threshold extraction of the CT data with all possible materials is not completely solved, which limits the surface overlap between photogrammetry data and CT data. Secondly, surface extraction algorithms also need to be more robust, due to the fact that the data of CT can be noisier for some occasions. Only a qualified CT surface data is meaningful for the colouring process. Otherwise, it will cause wrong colouring operations.

Future work and conclusion
Photogrammetry is a classical method for 3D reconstruction, while the quality is affected by the viewing configurations and the illumination. CT data preserves high precision and completeness regardless of light and viewing conditions, however, lacks colour information.
In this research, the proposed method provides a new perspective of integrating advantages of different 3D imaging technologies together rather than simply stitching the data. The method can be extended to different sensors and in different ways. The coloured CT data has not only the high precision and regular geometry, but carries now also the colour information.
The method can still be improved with regard to CT surface extraction results to get rid of the application limitations. Moreover, the ambiguous point exclusion can also be avoided if the ray trace geometry is solved in the depth image generation process.