Innovative Methodology for Multi-View Point Cloud Registration in Robotic 3D Object Scanning and Reconstruction

Abstract

The paper is concerned with the problem of multi-view three-dimensional (3D) point cloud registration. A novel global registration method is proposed to accurately register two series of scans into an object model underlying 3D imaging digitization by using the proposed oriented bounding box (OBB) regional area-based descriptor. A robot 3D scanning strategy is nowadays employed to generate the complete set of point cloud of physical objects by using 3D robot multi-view scanning and data registration. The automated operation has to successively digitize view-dependent area-scanned point cloud from complex-shaped objects by simultaneous determination of the next best robot pose and multi-view point cloud registration. To achieve this, the OBB regional area-based descriptor is employed to determine an initial transformation matrix and is then refined employing the iterative closest point (ICP) algorithm. The key technical breakthrough can resolve the commonly encountered difficulty in accurately merging two neighboring area-scanned images when no coordinate reference exists. To verify the effectiveness of the strategy, the developed method has been verified through some experimental tests for its registration accuracy. Experimental results have preliminarily demonstrated the feasibility and applicability of the developed method.

1. Introduction

Recently, automated three-dimensional (3D) object digitization, known under various terminologies such as 3D scanning, 3D digitizers or reconstruction, has been widely applied in many applications, such as 3D printing, reverse engineering, rapid prototyping and medical prosthetics. According to the sensing principle being employed, the current solutions can be generally classified into two main categories, namely hand-guided and automated scanning techniques. Hand-guided scanning allows for acquiring arbitrary shapes [1,2]. However, the effectiveness of this scanning method highly depends on the skills of the user and the scanning process is generally time-consuming. In contrast, the automated scanning solution using turntable-based 3D scanners is faster, less expensive, more automated, and easier to use. However, it is still not able to reconstruct objects with arbitrary or complex geometry [3,4]. To enhance the efficiency, the six-axis robot arm integrated with 3D imaging scanners has recently emerged as a technical developing trend for 3D surface scanning for objects with arbitrary or complex geometry [5,6,7]. Both Callieri [5] and Larsson [6] presented a system for automated 3D modeling consisting of a 3D laser scanner, an industrial six-axis robot and a turntable. The advantage of this system is that it automatically achieves the shape of the object from initial scanning and thus, if possible, will scan the object using an orientation of the scanner which gives the most accurate result. Meanwhile, an autonomous 3D modeling system, consisting of an industrial robot and a laser striper for efficient surface reconstruction of unknown objects was proposed by Simon Kriegel [7]. The system iteratively searches for possible scan paths in the local surface model and selects a next-best-scan (NBS) in its volumetric model. The data reconstruction of the robotic scanning systems reviewed above is mainly performed in two steps: data registration and data merge. First, an initial transformation matrix between two neighboring area scans can be traditionally determined by using the nominal scanning position of the robot arm or turntable as a starting reference, although it is not always accurate due to the unavoidable robot positioning errors. Thus, the matrix can serve as a reasonable initial estimate by using the Iterative Closest Point (ICP) algorithm [8] for further refinement of the transformation matrix. However, due to the lack of a global registration between successively scanned data, the existing scanning techniques still fail to robustly satisfy accurate object digitization with 100% surface coverage. Figure 1 shows two cases that are missing some parts of the object surfaces after performing automated object digitization using robot scanning.

The above issue is mainly caused by the fact that 3D sensor cannot measure point data from the scanned object surface that is in contact with its supporting ground. To resolve the issue, the object can be reoriented manually or automatically by using a robot grabber. Nevertheless, the changing orientation of the object produced by these two methods will definitively lose the initial pose estimate from before and after the orientation change. To estimate the pose transformation, some existing algorithms, such as local descriptor–based coarse alignment [9,10,11], use local-regional surface characteristics. As a general case, to compute the initial pose estimate, the geometric distances between corresponding 3D points existing in different views are minimized. The most common correspondences are points, curves and surfaces. The point signature [9] is a point descriptor for searching correspondences, which describes the structural neighborhood of a 3D point by the distances between neighboring points and their corresponding contour points. Johnson and Hebert proposed the concept of a spin image [10] to create an oriented point at a vertex in the surface, which is a two-dimensional (2D) histogram containing the number of points in a surrounding supporting volume. More recently, Tombari et al. [11] presented a 3D descriptor called the Signature of Histograms of OrienTations (SHOT) that encodes information about the surface within a spherical support structure. These methods generally work well for objects with unique local geometric characteristics. However, the methods can easily lead to less discrimination or sensitivity when dealing with object geometry having symmetrical or plural repeating surface contours. Therefore, modeling scanning completion with 100% surface coverage still remains one of the most challenging problems in current 3D object digitization.

Thus, accurately registering point cloud of two series of scanned point cloud into the reconstructed object model without an established coordinate reference has proven to be a difficult process [8,9,10,11,12]. To overcome the above difficulty, an automated robot scanning system integrated with a NTU-developed optical 3D measuring probe (NTU represents National Taiwan University) is developed on a global registration method to accurately register two series of scans into an object model underlying 3D imaging digitization by using the proposed oriented bounding box (OBB) regional area-based descriptor. The proposed descriptor is used to register two series of scans by matching the hardware and basic system configuration as well as data processing algorithms such as preprocessing, overlapping detection and robot path planning, which have been described in [12].

2. Global Registration Method Based on the OBB Regional Area-Based Descriptor

2.1. Principle Concept and Flow Chart Diagram of the Global Registration Method

To find the transformation between two series of scans that are acquired from multiple single scans, the OBB regional area-based descriptors of source point cloud are matched with the descriptors of the target point cloud. Figure 2 illustrates the concept of the proposed registration method based on the OBB regional area-based descriptor to register two series of scans which are represented by the red and green point cloud for the target and source 3D point cloud that are scanned and obtained from the optical probe. Each of the proposed descriptors includes two critical components. The first component contains the information about the OBB of the point cloud. Each OBB is represented by a datum corner, C(x_C, y_C, z_C), and three vectors, CC₁(x_max, y_max, z_max), CC₂(x_mid, y_mid, z_mid), CC₃(x_min, y_min, z_min), corresponding with the maximum, middle, and minimum dimensions of the OBB, respectively. The second component represents the spatial distribution of the surface area of the object in the OBB. The similarity between the regional area-based descriptor representing the source point cloud and the regional area-based descriptor representing the target one is then determined by using the normalized cross-correlation.

To make the proposed method clear in its operation procedure, the following flow chart diagram, shown in Figure 3, as well as Algorithm 1, is used to describe the proposed method in its main five steps. In step 1, the boundary points located near the scanning point missing zone can be detected through the geometric relationship with the support plane, using a set of connected edges that satisfy two conditions, in which all the boundary edges in the point missing zone have a similar orientation and each edge lacks an assigned triangle, either on the left or on the right side [7]. Then, the oriented bounding box (OBB), which is a rectangular bounding box that covers all the object point cloud, is determined in step 2. In the next step, the generation of the descriptors of the point cloud is performed. Then, the matching process between the descriptors is implemented and achieves the initial transformation that can be used for further refined registration using ICP.

Algorithm 1
Input: Two series of scans P and Q
Output: Transformation matrix
1	Estimate the boundary of the missing points and enclose the missing regions
2	Determine individual OBBs of the two series of scans
3	Generate the descriptors of P and Q
3.1	Calculate the OBB regional area-based descriptor of P
3.2	Translate the OBB of P to Q using the translation matrix defined by the centers of gravity of P with its estimated missing points and Q with its estimated missing points
3.3	Rotate Q around the three axes of the coordinate system defined by the center of gravity of Q as the origin and the three unit vectors considered as the three directions of the OBB of point cloud P with every increment Δθ on each axis, generating the OBB regional area-based descriptors of Q
4	Match the descriptors between P and every Q having an increment Δθ of rotation performed in Step 3.3
5	Determine the best machine and obtain the initial transformation matrix

2.2. Estimation of Missing Regions and Initial Translation

Assuming that two series of scans to be registered are represented as two sets of unorganized points P = {p₁, p₂,…, p_k} and Q = {q₁, q₂,…, q_k}, the initial transformation, which represents the coarse alignment between the two series, can be estimated as follows.

If P and Q represent the same free-form shape, the least squares method can be used to calculate the initial rotation and translation by minimizing the sum of the squares of alignment errors. Let E be the sum of the squares of alignment errors. Accordingly,

$$E={\displaystyle \sum _{i=1}^k{\Vert e_i\Vert }^2}={\displaystyle \sum _{i=1}^k{\Vert p_i-R\cdot q_i-t\Vert }^2}$$

(1)

where R is the rotation matrix, t is the translation vector, p_i and q_i is the corresponding point pair between the P and Q point clouds.

Let $\overline{p}$ be the centroid of the set of corresponding points in P {p₁, p₂, …, p_k} and $\overline{q}$ be the centroid of the set of corresponding points in Q {q₁, q₂, …, q_k}. Then $\overline{p}$ and $\overline{q}$ can be defined as:

$$\overline{p}=\frac{1}{k}{\displaystyle \sum _{i=1}^k{p}_i}$$

(2)

$$\overline{q}=\frac{1}{k}{\displaystyle \sum _{i=1}^k{q}_i}$$

(3)

The new coordinates of points are denoted by:

$$p{\text{'}}_i=p_i-\overline{p}$$

(4)

$$q{\text{'}}_i=q_i-\overline{q}$$

(5)

Note that:

$${\displaystyle \sum _{i=1}^{k}p{\text{'}}_i}={\displaystyle \sum _{i=1}^k\left(p_i-\overline{p}\right)}={\displaystyle \sum _{i=1}^k{p}_i}-k\overline{p}=0$$

(6)

$${\displaystyle \sum _{i=1}^{k}q{\text{'}}_i}={\displaystyle \sum _{i=1}^k\left(q_i-\overline{q}\right)}={\displaystyle \sum _{i=1}^k{q}_i}-k\overline{q}=0$$

(7)

The registration error term is then expressed and rewritten as:

$$e_i=p_i-R\cdot q_i-t=p{\text{'}}_i-R\cdot q{\text{'}}_i-t^{\prime }$$

(8)

where $t^{\prime }=t-\overline{p}+R\cdot \overline{q}$

The sum of the squares of the errors becomes:

$$\begin{array}{l}E={\displaystyle \sum _{i=1}^k{\Vert p{\text{'}}_i-R\cdot q{\text{'}}_i-t^{\prime }\Vert }^2}\\ ={\displaystyle \sum _{i=1}^k{\Vert p{\text{'}}_i-R\cdot q{\text{'}}_i\Vert }^2}-2t^{\prime }\cdot {\displaystyle \sum _{i=1}^k\left[p{\text{'}}_i-R\cdot q{\text{'}}_i\right]}+k{\Vert t^{\prime }\Vert }^2\\ ={\displaystyle \sum _{i=1}^k{\Vert p{\text{'}}_i-R\cdot q{\text{'}}_i\Vert }^2}+k{\Vert t^{\prime }\Vert }^2\end{array}$$

(9)

The first term in Equation (9) does not depend on $t^{\prime }$. Therefore, the total error is minimized if $t^{\prime }=0$, and then:

$$t=\overline{p}-R\cdot \overline{q}$$

(10)

Therefore, the initial translation matrix can be computed based on the difference between the centroid of P and the centroid of Q.

In addition, the weight loss at the missing parts of two the scan series may result in a significant bias between the translation matrix being estimated (by the difference between the centroid of P and the centroid of Q) and the real translation matrix. The shifting of two gravity centers (the red and green points shown in Figure 4) being created using a virtual camera from a model having the centroid (white point) in Figure 4a illustrates the influence of weight loss from the missing points to the centroid point position. Assume that the two point sets are separated as shown in Figure 4b. The accurate translation matrix to merge set Q (green point cloud) with set P (red point cloud) is defined through the difference between the current centroid of P and the original centroid of the green point cloud (green points). Therefore, it is necessary to estimate the original centroid of P on the object for translation.

The boundary points located near the missing part can be detected through the geometric relationship with the support plane and using a set of connected edges that satisfy two requirements, in which all edges have a similar orientation and each edge lacks an assigned triangle, either on the left or on the right side [7] (as shown in Figure 4c). Then, the Random Sample Consensus (RANSAC) algorithm [13] is used to fit a mathematical model to the boundary and enclose the missing region by adding extra points. Therefore, the initial translation matrix can be computed based on the difference between the centroid of enclosed P and the centroid of enclosed Q. From now on, P and Q are referred to as enclosed P and enclosed Q, respectively.

2.3. Determination of Oriented Bounding Box (OBB)

The oriented bounding box is a rectangular bounding box that covers whole object point cloud. Each OBB is represented by a corner, C(x_c, y_c, z_c), and three vectors, CC₁(x_max, y_max, z_max), CC₂(x_mid, y_mid, z_mid), CC₃(x_min, y_min, z_min), corresponding with the maximum, middle, and minimum dimensions of the OBB, respectively (as shown in Figure 5). The orientations of the OBB can be determined by using the covariance-based method, which is proposed by Gottschalk (1996) [14]. The algorithm is started by obtaining the centroid $\overline{p}\left(\overline{x},\overline{y},\overline{z}\right)$ of the object point cloud:

$$\overline{x}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{x}_i}\overline{y}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{y}_i}\overline{z}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{z}_i}$$

(11)

where n is the number of points in the object point cloud. The covariance matrix of the input point cloud COV is then computed as follows:

$$\begin{array}{cc}\hfill COV& =\frac{1}{n}{\displaystyle \sum _{i=1}^n\left(p_i-\overline{p}\right){\left(p_i-\overline{p}\right)}^T}\hfill \\ & =\frac{1}{n}\left[\begin{array}{lll}{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}\hfill & {\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}\hfill & {\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(z_i-\overline{z}\right)}\hfill \\ {\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}\hfill & {\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}\hfill & {\displaystyle \sum _{i=1}^n\left(y_i-\overline{y}\right)\left(z_i-\overline{z}\right)}\hfill \\ {\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(z_i-\overline{z}\right)}\hfill & {\displaystyle \sum _{i=1}^n\left(y_i-\overline{y}\right)\left(z_i-\overline{z}\right)}\hfill & {\displaystyle \sum _{i=1}^n{\left(z_i-\overline{z}\right)}^2}\hfill \end{array}\right]\hfill \end{array}$$

(12)

The above matrix has three real eigenvalues, λ₁ ≥ λ₂ ≥ λ₃, and three normalized eigenvectors, v₁(v₁₁, v₁₂, v₁₃), v₂(v₂₁, v₂₂, v₂₃), and v₃(v₃₁, v₃₂, v₃₃), respectively. These eigenvectors are considered as the three directions of the OBB. By projecting all points in the object point cloud onto the eigenvectors, three dimensions of the object can be determined as the distance between the nearest and farthest projected points in each eigenvector. Considering the coordinate system with the origin $\overline{p}$ and three unit vectors v₁, v₂, and v₃ (as shown in Figure 6), the coordinates of the point p_i(x_i, y_i, z_i) in this coordinate system are determined by following equation:

$$\left[\begin{array}{c}{x^{\prime }}_i\\ {y^{\prime }}_i\\ {z^{\prime }}_i\end{array}\right]=\left[\begin{array}{ccc}{v}_1{e}_1& v_1{e}_2& v_1{e}_3\\ v_2{e}_1& v_2{e}_2& v_2{e}_3\\ v_3{e}_1& v_3{e}_3& v_3{e}_3\end{array}\right]\left[\begin{array}{c}{x}_i-\overline{x}\\ y_i-\overline{y}\\ z_i-\overline{z}\end{array}\right]$$

(13)

$$\left[\begin{array}{c}{x^{\prime }}_i\\ {y^{\prime }}_i\\ {z^{\prime }}_i\end{array}\right]=\left[\begin{array}{ccc}{v}_{11}& v_{12}& v_{13}\\ v_{21}& v_{22}& v_{23}\\ v_{31}& v_{32}& v_{33}\end{array}\right]\left[\begin{array}{c}{x}_i-\overline{x}\\ y_i-\overline{y}\\ z_i-\overline{z}\end{array}\right]$$

(14)

$$\left[\begin{array}{c}{x^{\prime }}_i\\ {y^{\prime }}_i\\ {z^{\prime }}_i\end{array}\right]=T\left[\begin{array}{c}{x}_i-\overline{x}\\ y_i-\overline{y}\\ z_i-\overline{z}\end{array}\right]$$

(15)

where e₁, e₂, and e₃ are the three unit vectors of the original coordinate system.

Note that matrix T is the orthogonal matrix. It means that T⁻¹ = T^T.

Let:

$$\begin{gathered} {x^{\prime }}_{min}=\text{min}\left\{{x^{\prime }}_{\mathrm{i}},i=1,\mathrm{...},n\right\},{x^{\prime }}_{max}=\text{max}\left\{{x^{\prime }}_{\mathrm{i}},i=1,\mathrm{...},n\right\} \\ {y^{\prime }}_{min}=\text{min}\left\{{y^{\prime }}_{\mathrm{i}},i=1,\mathrm{...},n\right\},{y^{\prime }}_{max}=\text{max}\left\{{y^{\prime }}_{\mathrm{i}},i=1,\mathrm{...},n\right\} \\ {z^{\prime }}_{min}=\text{min}\left\{{z^{\prime }}_{\mathrm{i}},i=1,\mathrm{...},n\right\},{z^{\prime }}_{max}=\text{max}\left\{{z^{\prime }}_{\mathrm{i}},i=1,\mathrm{...},n\right\} \end{gathered}$$

The parameters of the OBB can be determined as follows:

$$\{\begin{array}{l}{x}_C=\overline{x}+v_{11}{x^{\prime }}_{min}+v_{21}{y^{\prime }}_{min}+v_{31}{z^{\prime }}_{min}\\ y_C=\overline{y}+v_{12}{x^{\prime }}_{min}+v_{22}{y^{\prime }}_{min}+v_{32}{z^{\prime }}_{min}\\ z_C=\overline{z}+v_{13}{x^{\prime }}_{min}+v_{23}{y^{\prime }}_{min}+v_{33}{z^{\prime }}_{min}\end{array}$$

(16)

$$\begin{array}{l}C{C}_1=\left({x^{\prime }}_{\max }-{x^{\prime }}_{\min }\right)v_1\\ C{C}_2=\left({y^{\prime }}_{\max }-{y^{\prime }}_{\min }\right)v_2\\ C{C}_3=\left({z^{\prime }}_{\max }-{z^{\prime }}_{\min }\right)v_3\end{array}$$

(17)

2.4. Generating Descriptors of Point Cloud

2.4.1. The OBB Regional Area-Based Descriptor

Each proposed feature descriptor includes two components [15]. The first component contains the information about the OBB of the point cloud. Each OBB is represented by a corner, C(x_C, y_C, z_C), and three vectors, CC₁(x_max, y_max, z_max), CC₂(x_mid, y_mid, z_mid), CC₃(x_min, y_min, z_min), corresponding with the maximum, middle, and minimum dimensions of the OBB, respectively. The second component represents the distribution of the surface area of the object in the OBB. We assume that the three directions of the OBB are divided into k₁, k₂, and k₃ segments (as shown in Figure 7a); the total number of subdivided boxes in the OBB is n = k₁ × k₂ × k₃. The surface area of the object in each subdivided box V_ijk (i = 0, ..., k₁ − 1; j = 0, ..., k₂ − 1; k = 0, ..., k₃ − 1) is denoted as f_v, where v = k(k₁k₂) + j(k₁) + I; f_v can be determined by subdividing the considered surface into the triangle mesh [16] (as shown in Figure 7c). The distribution of S_v is defined as the regional area-based descriptor, shown in Figure 7d. Thus, the feature descriptor FV can be described as follows:

$$FV=\{C,C{C}_1,C{C}_2,C{C}_3,f_0,\dots ,f_v,\dots ,f_{n_V}\}$$

(18)

where n_V = k₁k₂k₃ − 1.

2.4.2. Generating Database and Calculating Descriptors

Firstly, the regional area-based descriptor of point cloud P is computed as FV_P = {C_P; C_PC_1P, C_PC_2P, C_PC_3P, f_Pv, v = 0, …, n} in the OBB of P. Then we shift the OBB of P to Q using the translation matrix defined by the centers of gravity of P and Q. We consider the coordinate system O_Qx_qy_qz_q with the origin O_Q is the center of gravity of Q and three unit vectors ν₁, ν₂, and ν₃, which are considered as the three directions of the OBB of point cloud P. Rotating Q around the three axes of the coordinate system O_Qx_qy_qz_q is performed with every increment Δθ on each axis. For each rotation, the regional area-based descriptor Q in the OBB of P is FV_Qi = {C_Qi; C_QiC_1Qi, C_QiC_2Qi, C_PC_3Qi, f_Qv, v = 0, …, n} which can be then determined as described as Algorithm 2. To enhance the accuracy, the descriptors FV_P and FV_Qi can ignore the subdivided boxes in which the missing parts of point cloud P are located.

Algorithm 2. Generation of database

Input: Series of scans Q = {q1, q2,…, qn}, the OBB of P

Output: A set of descriptors of Q at different poses

for l = 0, 1, 2, ..., 2π/Δθ

for j = 0, 1, 2, ..., π/Δθ

for k = 0, 1, 2, ..., (π/2)/Δθ

θxl = l*Δθ, θyj = j*Δθ, θzk = k*Δθ

R = Rx (θxl) Ry (θxj) Rz (θzk)

$Q\stackrel{R}{\to }{Q}_i$

FVQi = {fQiv, v = 0, …, n}

end_for

end_for

end_for

2.5. Matching Descriptors

In the first step of the matching process, the OBB’s parameters of point cloud P are compared with the OBB’s parameters of point cloud Qi. If the OBB matching satisfies the given condition, the regional area-based descriptor of point cloud P is then compared with the regional area-based descriptor of Qi. These two OBBs should be satisfied with the following equation:

$$\{\begin{array}{l}\frac{C_P{C}_{P1}.C_Q{C}_{Q1}}{\Vert C_P{C}_{P1}\Vert .\Vert C_Q{C}_{Q1}\Vert }<{\alpha }_{thresh}\\ \frac{C_P{C}_{P2}.C_Q{C}_{Q2}}{\Vert C_P{C}_{P2}\Vert .\Vert C_Q{C}_{Q2}\Vert }<{\alpha }_{thresh}\\ \frac{C_P{C}_{P3}.C_Q{C}_{Q3}}{\Vert C_P{C}_{P3}\Vert .\Vert C_Q{C}_{Q3}\Vert }<{\alpha }_{thresh}\end{array}$$

(19)

where α_thresh is the given adequate threshold.

The similarity between the regional area-based descriptors that represent the object point clouds P and Q is determined by using the normalized cross-correlation. The normalized cross-correlation between F_P = {f_Pv, v = 0, …, n} and F_Q = {f_Qv, v = 0, …, n} is computed as follows:

$$C\left(F_P,F_Q\right)=\frac{{\displaystyle \sum _{v=0}^n\left(f_{Pv}-{\overline{f}}_P\right)\left(f_{Qv}-{\overline{f}}_Q\right)}}{\sqrt{{\displaystyle \sum _{v=0}^n{\left(f_{Pv}-{\overline{f}}_P\right)}^2}\cdot {\displaystyle \sum _{v=0}^{n_{}}{\left(f_{Qv}-{\overline{f}}_Q\right)}^2}}}$$

(20)

where ${\overline{f}}_P=\frac{1}{n+1}{\displaystyle \sum _{v=0}^n{f}_{Pv}}$ ${\overline{f}}_Q=\frac{1}{n+1}{\displaystyle \sum _{v=0}^n{f}_{Qv}}$

The matching is the best in Equation (20) in which the normalized cross-correlation $C\left(F_P,F_Q\right)$ reaches its peak.

2.6. Transformation Estimation and Refinement

After the best matching is defined, the correspondence feature vectors FV_P and FV_Q can be determined. Based on these feature vectors, the initial transformation matrix T_initial between two series of scans can be estimated by aligning the frame {C_Q, v_Q1, v_Q2, v_Q3} that represents the series of scans 1 to the frame {C_P, v_Qbest1, v_Qbest2, v_Qbest3} that represents the series of scans 2. Although the parameters of the transformation matrix can be obtained in the initial registration, the accuracy may not be satisfactory due to the error caused by the limited number of the rotating Q iterations. Fortunately, a refined registration such as the iterative closest point (ICP) algorithm can be performed to achieve precise registration between multi-view point cloud.

3. Results and Discussion

3.1. Case Study on Measured Data

To verify the feasibility of the developed methodology, some experiments were performed and evaluated. In the experiments, the point cloud of a hammer head were acquired from the NTU-developed optical 3D measuring probe being developed by using a random speckle pattern projection and triangulation measurement principle [12,15]. The dimensions of the hammer head are approximately 110 × 35.5 × 21 mm³. Firstly, the hammer head was placed on a fixed table, allowing for viewing the object from several positions around the object. However, the whole object surface could not be detected due to optical occlusion. Therefore, after the first series of automated scans to obtain point cloud P, the object was reoriented manually and the second series of automated scans was continuously carried out to acquire point cloud Q. Figure 8 illustrates the coarse registration process between the two series of scans. The width of the normalized cross-correlation peak depends on the number of iterations and Δθ, as shown in Figure 9. Meanwhile, the height of the normalized cross-correlation depends on the number of OBB segments k₁ × k₂ × k₃ and Δθ.

The performance of the registration can be estimated by using the distance from each overlapping point q_i in point cloud Q to the fitting control point p_i which is the projected point of q_i onto the triangle mesh of P. If d_i denotes the distance between a point p_i in P and its closest neighbor point q_i in Q, the mean distance μ and standard deviation σ, which are used to evaluate the performance of the object registration, are computed as follows:

$$\mathsf{\mu }=\frac{1}{n}{\displaystyle \sum _{i=1}^n{d}_i}$$

(21)

$$\sigma =\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(d_i-\mu \right)}^2}}{n-1}}$$

(22)

In this experiment, the mean deviation value is 0.0347 mm and the standard deviation value is 0.0235 mm. Figure 10 illustrates the other examples being achieved by the developed method.

3.2. Case Study on Synthetic Data

In this section we provide a discussion on the influence of parameters on the performance of the proposed method. In addition, a comparison against several already published local feature descriptor methods is also introduced. The considered approaches are: Spin Images (SI), Point Signatures (PS), Point Feature Histograms (PFH), and the Signature of Histograms of OrienTations (SHOT). All methods were implemented in C++ using the open source Point Cloud Library (PCL). To estimate the computation cost of the registration process, the experiments are processed on a computer with an Intel core i7 processor with 3.40 GHz and 8 GB RAM.

In real application, it is very difficult to distinguish different error sources such as shape measurement error (noise, surface sampling, etc.), correspondence error (occlusion, outliers, etc.), and registration error. In order to evaluate the coarse registration errors, we have generated synthetic data, which is provided by the Industrial Technology Research Institute (ITRI) (as shown in Figure 11). The dimensions of the socket model, connector model, cylinder model, and Brazo Control model are 45 × 25 × 25 mm³, 41 × 33 × 25 mm³, 35 × 35 × 35 mm³, and 86 × 53 × 30 mm³, respectively. Then, the point cloud corresponding with different views of the models are extracted by the NTU-developed virtual camera software (Precision Metrology (PM) lab, National Taiwan University, Taipei, Taiwan). The motion between green and red point cloud is a displacement of 50 mm in each axis and a rotation of π/4 around a fixed axis. Therefore, we have precise knowledge of the motion parameter rotation matrix R and translation vector t to serve as a reference for error estimation and validation of the coarse registration methods. The measures used to determine the accuracy include the rotation error and translation error. In order to evaluate the accuracy, the estimated transformation is compared to an expected one.

Results obtained show that sampling does not considerably affect the accuracy of the registration results. Thus, low resolution surfaces can be used in the proposed approach to reduce the computation time. Meanwhile, as is shown in

Table 1. Experimental results using the connector model obtained by the proposed coarse registration method with different samplings and increments Δθ.

Points	Increment Δθ (rad)	Translation Error (mm)	Rotation Error (rad)	Time (s)
1000	0.05	1.371	0.036	1.063
	0.02	1.371	0.017	15.328
	0.01	1.371	0.009	126.439
5000	0.05	1.352	0.039	7.244
	0.02	1.352	0.013	110.212
	0.01	1.352	0.012	904.237
10,000	0.05	1.368	0.035	18.332
	0.02	1.368	0.015	203.795
	0.01	1.368	0.009	1762.824

, errors in the proposed registration algorithms extremely depend on the rotation increment Δθ (rad). To enhance the accuracy, the rotation increment should become smaller. However, the increment has a significant impact on the runtime of the proposed coarse registration process. As an example, the experiment with 1000 points in Table 1, with computation time in the case Δθ = 0.01 (rad), at 126.439 (s), was much higher than that in the case Δθ = 0.05 (rad), at only 1.063 (s).

Another important characteristic of the developed method is the total number of subdivided boxes in the OBB, n = k₁ × k₂ × k₃. From

Table 2. Experimental results using the connector model obtained by the proposed coarse registration method with different numbers of the oriented bounding box (OBB) segments k₁ × k₂ × k₃.

Points	OBB Segments: k1 × k2 × k3	Translation Error (mm)	Rotation Error (rad)	Time (s)
1000	4 × 5 × 6	1.371	0.057	4.052
	7 × 8 × 9	1.371	0.017	15.328
	10 × 11 × 12	1.371	0.017	50.953
5000	4 × 5 × 6	1.352	0.073	43.570
	7 × 8 × 9	1.352	0.013	110.212
	10 × 11 × 12	1.352	0.013	329.641
10,000	4 × 5 × 6	1.368	0.035	67.327
	7 × 8 × 9	1.368	0.015	203.795
	10 × 11 × 12	1.368	0.015	936.841

, it can be seen that if the OBB contains many subdivided boxes, it is more expensive to find the best rotation. For a fast registration procedure, fewer subdivided boxes are preferable. However, it is also required that the number of subdivided boxes is great enough to achieve the expected rotation. For instance, in the experiment using a connector model of 1000 points, the rotation error rose to 0.057 (rad) for the OBB of n = k₁ × k₂ × k₃ = 4 × 5 × 6, while it reached just 0.017 (rad) if the OBB of n = k₁ × k₂ × k₃ = 7 × 8 × 9.

Besides, the ratio of the overlapped area is especially crucial. As shown in

Table 3. Experimental results using the connector model obtained by the proposed coarse registration method with different ratios of overlapped area.

Models	Ratio of Overlapped Area	Translation Error (mm)	Rotation Error (rad)	Time (s)
Connector	90%	1.066	0.013	15.328
	80%	1.371	0.017	13.543
	60%	2.452	0.042	10.476
	50%	5.796	0.137	9.548
	30%	12.347	0.532	7.365
Brazo	90%	1.023	0.012	16.462
	80%	1.455	0.015	13.443
	60%	2.674	0.044	11.348
	50%	5.357	0.159	8.348
	30%	15.642	0.758	6.769

, translation and rotation errors in the proposed approach grow in direct proportion to the percentage of the non-overlapping region. This change is especially significant with over 40% of the outliers. This is because it is difficult to estimate translation exactly without a high level of overlapping. The excellent results are obtained for over 80% of shape the overlapping. In this case, the translation and rotation for tested models were at under 1.5 (mm) and 0.02 (rad), respectively. As described in the Introduction, the proposed registration method is developed to register two series of surface scans which exist with a certain degree of optical occlusion from each individual viewpoint of the probe scanning. To accurately register two series of scans under such a circumstance, it is important to have them overlap as much as possible. However, realistically, the overlapped ratio really depends on the complexity of the scanned object’s geometry and how the object is placed on the rotation table. In general, the higher the overlapped ratio between two series of scans, the more accurate the registration that can be achieved.

The result of

Table 4. Experimental results using synthetic data obtained by the proposed method and local feature descriptor methods.

Models	Method	Translation Error (mm)	Rotation Error (rad)	Time (s)
Brazo	Spin Images (SI)	5.893	0.074	6.762
	Point Signatures (PS)	3.642	0.056	103.433
	PFH	4.543	0.047	43.232
	SHOT	2.632	0.053	5.432
	Proposed	1.371	0.017	13.543
Connector	Spin Images (SI)	7.443	0.093	8.342
	Point Signatures (PS)	5.378	0.053	123.533
	PFH	6.313	0.077	48.436
	SHOT	3.472	0.023	7.922
	Proposed	1.174	0.014	17.468
Cylinder	Spin Images (SI)	6.421	1.235	9.672
	Point Signatures (PS)	8.432	2.348	103.573
	PFH	3.445	1.764	34.562
	SHOT	5.472	0.973	4.358
	Proposed	1.421	0.013	12.667
Socket	Spin Images (SI)	9.243	2.645	9.457
	Point Signatures (PS)	9.831	2.569	103.376
	PFH	5.219	1.323	52.436
	SHOT	4.642	1.873	5.782
	Proposed	1.684	0.015	13.578

reveals the performances of the proposed approach and 3D key point descriptors. It shows that the developed method in this paper outperforms in registration accuracy with reasonable operation efficiency. As we know, one significant common problem using feature descriptors is that they usually fail to deal with object geometry that has symmetrical or plural repeated surface contours. We demonstrate the developed approach by using four test cases comprising two regular shapes, a Brazo and a connector, and two featureless models, a cylinder and a socket with 80% shape similarity which was provided by the Industrial Technology Research Institute (ITRI) (as shown in Figure 11). Given the reported results, it is clear that the proposed approach performs more accurately than the other methods on all the tested data. The computation efficiency of the proposed method is approximately the half of the fastest method, the Signature of Histograms of OrienTations SHOT. However, its time efficiency still ranks within a reasonable range.

To test the noise influence, the synthetic experiments were carried out with respect to three different levels of Gaussian noises: σ = 0.01, 0.02 and 0.03 (standard deviation). The coarse registration results on the connector point-sets are shown in

Table 5. Experimental results using the connector model obtained by the proposed coarse registration method with different levels of Gaussian noise.

Points	Gaussian Noise (σ)	Translation Error (mm)	Rotation Error (rad)	Time (s)
6000	0.01	1.339	0.018	20.862
	0.02	1.481	0.019	20.277
	0.03	1.585	0.021	20.839

. We found that the proposed approach is able to achieve acceptable results with a certain level of noise.

4. Conclusions

In this study, a global registration method is developed to overcome one of the most challenging problems remaining in current 3D scanning systems for accurate image registration, even with no reference coordinate existing between the neighboring scanned surface patches. The key innovation in this work is the strategy in determining a robust registration transformation between two neighboring scans. The experimental results demonstrate the validity and applicability of the proposed approach. The registration error can be controlled within a few tens of micrometers in one standard deviation for an object size reaching 100 millimeters.