PHOTOGRAMMETRY AND COMPUTED TOMOGRAPHY POINT CLOUD REGISTRATION USING VIRTUAL CONTROL POINTS

In this paper we propose a virtual control point based method for the registration of photogrammetry and computed tomography (CT) data. Because of the fundamentally different two data sources, conventional registration methods, such as manual control points registration or 3D local feature-based registration, are not suitable. The registration objective of our application is about 3D reconstructions of gyroscopes, which contain abundant geometric primitives to be fitted in the point clouds. In the first place, photogrammetry and CT scanning are applied, respectively, for 3D reconstructions. Secondly, our workflow implements a segmentation after obtaining the surface point cloud from the complete CT volumetric data. Then geometric primitives are fitted in this point cloud benefitting from the less complex cluster segments. In the next step, intersection operations of the parametrized primitives generates virtual points, which are utilized as control points for the transformation parameters estimation. A random sample consensus (RANSAC) method is applied to find the correspondences of both virtual control point sets using corresponding descriptors and calculates the transformation matrix as an initial alignment for further refining the registration. The workflow is invariant to pose, resolution, completeness and noise within our validation process.


INTRODUCTION
Sensor fusion is an important topic in many fields, because in real applications it is quite often difficult for a single sensor alone to provide the complete desired information. Therefore, the advantages of different sensors are integrated to strengthen the data characteristics. In the 3D digitization field of Tech Heritage (TH), the Gyrolog project (Fritsch et al., 2018;Fritsch et al., 2021) funded by the Federal Ministry of Education and Research of Germany, has innovatively introduced different methodologies, such as Photogrammetry, Computed Tomography (CT) as well as Endoscopy for the 3D digitization of the gyroscopic instrument collection of the University of Stuttgart. The combination of photogrammetry and CT has been discussed frequently in the medical field (Bolandzadeh et al., 2013), however, has not yet been applied often in TH digitization applications (Zhan et al., 2020). This paper mainly discusses a new registration method of photogrammetric and CT point clouds in such TH applications, where complete models are required. A point cloud is chosen as the common representation for photogrammetric surface data and CT volumetric data, due to the fact, that point cloud registration is an ongoing topic in the field of photogrammetry and computer vision, with various methods being put forward. The most frequently applied method is the iterative closest point (ICP) registration (Besl and McKay, 1992) or it's variants, which iteratively calculates the discrepancy of the overlap between two point clouds. Despite the wide application and continuous research with many algorithm variants, the requirement of sufficient initial registration makes it only suitable for the refinement of the registration. Methods based on automatic 3D features extraction such as the fast point feature histogram (FPFH) (Rusu et al., 2009) or nor- * Corresponding author mal aligned radial feature (NARF) (Rusu and Cousins, 2011) etc. have problems in the computation efficiency and even in detecting validate feature correspondences. These problems occur due to the original discrepancy between photogrammetry and CT data, such as density, edge and material characteristics, and the incomplete representation of CT scans. There are also works using geometric primitives such as (Alshawa, 2007) applies lines instead of points to implement an iterative process. But the application is intended to solve the registration of topographic terrestrial laser scanner data. (Yang et al., 2016) introduced a registration method based on semantic feature points, which also rely on the line feature and works for terrestrial laser scanning data. (Stamos and Leordeanu, 2003) calculates the intersection lines of neighboring planes and estimate the transformation matrix with at least two corresponding line pairs. (Theiler et al., 2012) put forward a terrestrial point clouds registration method by using virtual tie points from the intersection of plannar surfaces. However, the method deals with point clouds from the same source and a relatively easier situation for planes extraction. All the listed works are generally to align the point clouds based on primitive based information, they are more focusing on the same data type with same characteristics and lack the adaptability for multi-source data as well as complex structured 3D models. Other possibilities such as the artificial landmarks (Ayoub et al., 2007;Xin et al., 2013), however, have various drawbacks: (a) CT data contains no texture information, only geometry information could be used; (b) some of the surface information will be incomplete due to the characteristics of the CT scan; (c) unsharpness of corners for photogrammetry data which limits the accuracy of the manual control point picking process.
The proposed method below takes advantages of the characteristic of the artificial objects containing many regular geometry shapes such as cylinders, spheres or simply planes. For a best-fit process, point cloud segmentation is necessary to increase both the efficiency and precision. In our case, the region growing method (Vo et al., 2015) using just the normal and other attributes such as the color, delivers good results. Most frequent geometric primitives, such as planes, are fitted and used for control point determination by intersections. In addition, the primitives such as spheres or 3D circles could directly provide their centre points as control points. The correspondence of the extracted control point information from both datasets could be estimated via the RANSAC method (Buch et al., 2013). Afterwards, the registration is refined, based on the initial virtual control points selection.

METHOD
In the proposed workflow as shown in Figure 1, point clouds from CT and photogrammetry are generated from the multiple CT slices and multi-view imagery, respectively. After obtaining a CT surface point cloud, we implement the segmentation for the purpose of fitting geometry primitives for both datasets. The segmentation applies region growing methods based on different principles, such as color information, curvatures etc. In the next step, the virtual control points are calculated to estimate the transformation matrix with an refinement as a follow-up.

3D reconstruction
Generally 3D reconstruction is the technology of determining the 3D digital representation of a real object. The principle of photogrammetry is to use the overlapping information between multiview images to solve the 2D to 3D projection, and finally obtain the 3D surface model with texture. The quality of the final 3D model depends on factors such as the texture homogeneity of the object, the redundancy and the resolution of the photos, which reflect on the accuracy or even the completeness of the final 3D model. CT scanning records the attenuation information in each voxel space by X-rays passing through the object to reflect the internal structure and material properties of the object. The final 3D CT volumetric model consists of voxels with penetrating information. Though having the ability to recover the internal information in a non-destructive way, CT 3D reconstruction suffers from the noise and beam hardening artefacts and scattering. More details of 3D reconstruction as well as the CT surface extraction process could be referred to (Zhan et al., 2020).

Primitive fitting
The preliminary step of primitive fitting is point cloud segmentation, which has not been mentioned in most other point cloud registration methods using geometric primitives. However, for a watertight 3D model of a complex structured gyroscope object, the fitting process will be influenced by the spatial distribution of the points. Therefore, it is of vital importance to implement the point cloud segmentation into different clusters before fitting the primitives. Region growing is a conventional segmentation method, which is also validated in this work. A CT point cloud has high geometric accuracy regarding the voxel positions as well as the normals. Therefore, a region growing algorithm using the normals and curvatures works well in CT surface point cloud segmentation. As for photogrammetry, the primitives generally differ from each other in material and color.
Hence the color-based region growing segmentation is applied for the photogrammetric 3D model. For each segmented cluster, parameterized geometric primitives are used to fit to the discrete point cloud data. Like other fitting problems, various principles are available but mostly originates from RANSAC. Due to the diversity of the point clouds, a simple threshold definition is hardly to meet the requirements for a best primitive fitting problem. Maximum Likelihood Consensus (MLESAC), which improves the RANSAC by optimizing the inlier scoring rather than simply counting the inlier numbers, is adopted. The coefficients of the most frequent 3D primitives are listed as

Correspondance Estimation
With good fitting results, these parameterized primitives can be used to directly or indirectly obtain virtual control points. Because the fitting process is based on the RANSAC principle, the virtual control points are the representation of the overall characteristics of the point cloud clusters. Theoretically, they are more reliable than the manually selected control points or the feature points calculated using the neighboring information from the incomplete point cloud. After obtaining the virtual control points, their corresponding relationship needs to be solved. In general, after obtaining the key points, the corresponding descriptor will be calculated for the process of correspondence estimation. Global and local descriptors use global and neighborhood point information respectively, but these are not applicable to calculate virtual points, due to the fact that virtual points are outside the point cloud. The information of geometric primitives, which is used for virtual control point calculation, could be encoded. The virtual control points could be obtained from the center of a sphere, the intersection of the cylinder axis and the plane, or the intersection of three planes. (Theiler et al., 2012) applied mainly the plane angle for the condition of three planes intersection. If (a1, b1, c1) and (a2, b2, c2) stand for normal vectors of two intersecting planes, the angle between can be expressed as (1), α = arccos a1 · a2 + b1 · b2 + c1 · c2 a 2 1 + b 2 1 + c 2 1 · a 2 2 + b 2 2 + c 2 2 (1) We denote planes that are related to the intersection point as target planes, which can be ordered by the angles, and the rest non-target planes. Except the angles between target planes, the angles between each target plane and all non-target planes could also be calculated. All angle values are binned into a histogram. A similar to the idea of a point feature histogram (PFH), the binning process divides the angles into certain subdivisions, and counts the number of occurances in each subinterval. A similar process could also be applied to the cylinder axis and plane intersection situation by replacing the plane angle to line-plane angle. An example is shown in Figure 2.

Error analysis
After the correspondence estimation, the transformation matrix can be directly calculated to serve as an initial match. The next method is generally through ICP. ICP is mainly an iterative loop to calculate the nearest neighbors of two sets of point clouds, and usually good results can be obtained under good initial matching conditions. In the iterative process, the sum of the squared errors, as shown in (2) is minimized.
where Np is the number of points, xi and pi are corrresponding points, R and t are rotation and translation, respectively. Though ICP is widely used, it is a well-known fact, that most implementations of ICP do not consider precision measures as given in the field of photogrammetry. The alternative could be based on the least squares method. Here we use the wellknown Gauss-Helmert model to solve the overall 7 parameters transformation problem as shown in (3). In addition, for providng a best fusion model, a complete accuracy evaluation result using the law of error propagation can also be obtained with the assumption of precision values for both data sets, the target and the source data.
where X and x stand for target and source points, µ, X0 and R are scale, translation and rotation respectively. More detailed mathematical derivation could be referred to (Fritsch et al., 2021).

Object
The main experimental object for this study is a pneumatically driven directional gyro, manufactured by GM corp in US as shown in Figure 3. The object is chosen due to its rich geometric primitives design while few evident corner points within the rounded edges could be manually picked for CT and photogrammetry point cloud registration.

CT and photogrammetry 3D reconstruction
The sensors used for CT scanning and photogrammetry can be referred to (Zhan et al., 2020). As for photogrammetric 3D reconstructions, Figure 4 shows the camera poses represented by circular white rectangles via SfM (Structure-from-Motion) and the watertight 3D surface model from DIM (Dense Image Matching) and texturing procedure. CT is used to reconstruct internal structures of an object that is not visible from outside. Figure 5(a) shows the complete CT reconstructed 3D model visualized in VGStudio (VGStudio, n.d.). Figure 5(b) is the point cloud of CT data by representing each voxel by the center point with its intensity value.

CT 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. Thereafter a CT surface extraction is necessary to provide the similar information as given by the photogrammetric point cloud for registration. We apply a ball query method for each point to extract the points with less neighbouring points within a threshold range, which are interpreted as the points on the surface of the 3D CT model. An alternative method finding the convex hull could be found in (Zhan et al., 2020). Figure 6(a) shows the surface extraction result. In addition, by removing this shell of certain thickness will also avoid the visualization ambiguity of the integrated model on the surface.

Point cloud segmentation and primitive fitting
After obtaining a CT surface and photogrammetric point cloud, the implementation of a direct geometry primitive fitting may result in unsatisfying output, due to the highly complex structure of our objects, which are gyroscopes. Therefore, a point cloud segmentation as shown in Figure 6(b) using region growing methods based on different principles is necessary beforehand, considering the differences of the sensors as well as the object structure characteristics. A well-segmented point cloud contributes to an easier and more precise best-fit of geometric primitives based on RANSAC. The plane, as the most frequent geometric primitive in gyroscopes, is determined by the inliers, i.e. extracted co-planar points with a predefined threshold. Figure 7 displays examples of plane and cylinder fitting using the Point Cloud Library (PCL) (Rusu and Cousins, 2011). A control point could be determined by the intersection of every three non-parallel parameterized planes with precise coordinates. In addition, other best-fit primitives such as spheres or 3D circles could provide their center point directly as control points.

Transformation matrix estimation
In the next step, two sets of points are obtained as control points for the estimation of the transformation parameters as shown (a) plane fitting (b) cylinder fitting Figure 7. Primitives fitting on photogrammetry model in Figure 8. For the differences between two data-sets as well as the primitive fitting procedures, the control points may be disordered and different in quantity. To determine the transformation, the control points together with their descriptors as introduced in 2.3 are to be estimated for corresponding point pairs. The obtained result is used for the transformation matrix cal- Figure 8. Virtual control points of photogrammetric model culation as initial registration for the later refinement by ICP and/or the Gauss-Helmert Model. With two transformation matrices, the photogrammetry point cloud and the CT point cloud without the outer shell could be registered together. In our experiment both refinement methods delivered good registration results without too much differences. Cross sections of the photogrammetry and CT models are displayed in Figure 9. Additionally, more detailed precision measures could be obtained by the Gauss-Helmert Model estimation. The ground sampling distances (GSD) of the photogrammetry and CT for our application are 0.05-0.09mm and 0.06mm respectively. With 25 corresponding joint control points, the estimated standard deviation of unit weight σ = 1.03, the estimated precision of CT σCT = 0.07mm, and the estimated precision of photogrammetry is σ phot = 0.06mm. On the one hand, a high precision pose correspondence is realized between two data sources, and on the other hand the integrated model containing both a clean colored surface and internal structure information is derived.

Summary the findings
The proposed workflow takes into consideration the characteristics of the point clouds from different sensors for designing appropriate steps for a good registration. Due to the difference of CT volumetric data from the normal surface point cloud, the CT data is transformed into the point cloud format and the surface extraction is implemented to be fitted into the point cloud registration framework. As for the common incomplete representation of the surface of CT scanning data, virtual control points generated by the fitted primitives rather than the local information depending 3D features are applied for the transformation estimation.

Strength and weakness analysis
The proposed workflow has advantages in aspects (a) robustness against the noise of the point cloud; (b) robustness against incomplete CT data due the low penetration of some materials; (c) high efficiency on condition of good primitive fitting results. Except the proposed workflow, registrations with other principles such as 3D feature-based and manual control point-based registration are also under investigation. FPFH-based coarse registration suffers from the low efficiency. The relationship between the computing time and the number of points is plotted and the registration is shown in Figure 10 and 11. Though the selection of manual control points works for coarse registrations, it is very time consuming and delivers randomness in the process, with unsharp corners of the photogrammetry data. The proposed workflow involves several steps to reach the final registration. However, each single step could be improved, such as well-extracted CT surfaces, more robust point cloud segmentations, more adaptive primitive fittings regarding different objects as well as a finding more efficient descriptors for the virtual control points.

Outlook
With regard to the analysis in section 4.2, the automation of the workflow of each single step is necessary and can be improved.
In addition, a more detailed comparison regarding other registration strategies should be implemented to validate the proposed method, for general application scenarios.