2.1. CAD-Based Reconstruction and Fitting

Several CAD model–based fitting methods have been proposed for industrial robots.

Shah et al. ⁶ applied simulated annealing to fit multipart assembly models to point clouds and initialized the poses using oriented bounding boxes (OBB). However, real factory scans often contain nonmodeled elements, such as wire harnesses, bolts, and brackets, which distort the OBB and degrade the initialization accuracy. When the initial pose is inaccurate, the optimization may fail to converge to the global optimum.

Hu et al. ⁷ proposed a deep-learning-based framework using predefined CAD templates for object recognition and parameter estimation. Although effective under controlled conditions, articulated robots in automotive factories exhibit significant variations owing to their multiple links, subordinate components, closed-loop structures, and interchangeable end effectors. It is difficult to train a classifier that can accommodate all these variations.

Kawasaki et al. ⁸ introduced a fitting method for articulated robots; however, it required manual base-link initialization, increased the processing time, and increased the risk of human error. Their method also does not explicitly address nonmodeled components, resulting in unstable performance in real TLS scans. Furthermore, most existing fitting approaches implicitly assume single-chain link structures, limiting their applicability to the branched-chain mechanisms commonly found in conveyance systems.

2.2. ICP-Based Point-Cloud Fitting

The Iterative Closest Point (ICP) method ⁹ and its variants are among the most widely used methods for aligning two point clouds. Numerous extensions have been proposed; however, they share the common principle of iteratively estimating point correspondences and minimizing the geometric discrepancies between them. In its basic form, the ICP repeatedly assigns each point in the source cloud to the nearest point in the target cloud and computes the rigid transformation that best aligns the two sets.

Although effective for many rigid objects, the application of conventional ICP-based methods to articulated machines introduces several difficulties. First, ICP uniformly treats all points as equally informative, even though factory scans contain nonmodeled elements, such as wire harnesses, hoses, and fasteners, as well as occlusions and sensor noise. These factors generate geometric outliers that can significantly degrade the registration accuracy.

Second, ICP inherently assumes a single rigid body motion. By contrast, robots and conveyance systems consist of multiple links connected by joints, each of which requires its own motion. Independently estimating the pose of each link is challenging because naïvely applying ICP to individual parts can lead to inconsistencies between adjacent links. In articulated mechanisms, joint axes must remain coincident, and joint constraints must be satisfied; without incorporating such kinematic constraints, ICP solutions may violate the machine’s actual degrees of freedom or produce mechanically impossible configurations.

Finally, although many ICP variants attempt to mitigate noise and outliers, they still rely primarily on surface points that are heavily affected by nonmodeled attachments. In contrast, geometric edges, sharp intersections between links, or structural features, remain consistent between CAD models and real machines, making them more reliable cues for registration.

2.3. Edge Detection in Point Clouds

Point-cloud registration can be improved by exploiting geometric features rather than relying solely on the raw point positions. In this study, we focused on edge features that provide stable and distinctive geometric cues, even in cluttered factory environments. This section reviews existing methods for extracting edges from 3D point clouds and their applications to the registration process.

Extensive studies have been conducted on point-cloud edge detection. Hackel et al. ¹⁰ introduced a geometry-based contour detector to evaluate sharp changes in surface orientation; however, it relies on expensive 3D neighborhood searches and structured optimization. Pauly et al. ¹¹ proposed a multiscale PCA-based technique for extracting line-type features from point-sampled surfaces. Although this method performs well for high-density and relatively clean data, its susceptibility to noise limits its effectiveness for real-world factory scans. More recently, Ben-Shabat and Gould ¹² proposed DeepFit, which learns point-wise weights using a neural network, and applies weighted least squares for surface fitting. Although TLS is accurate and robust to density variations, the need to infer weights through a network introduces substantial computational costs, particularly for large-scale TLS point clouds.

In addition to feature extraction, several studies have leveraged linear features for point-cloud registration. Al-Durgham et al. ¹³ developed an automated random sample consensus (RANSAC)-based framework that exploits the invariant spatial separation and angular deviation between 3D line pairs to establish correspondences. Habib et al. ¹⁴ proposed a methodology that enforces the coincidence of conjugate linear features to integrate multisource data, including photogrammetric and LiDAR datasets, into a common reference frame. For large-scale urban environments, Stamos and Leordeanu ¹⁵ introduced a feature-based approach that organized segmented range data into a topological graph for scan placement. Furthermore, Date et al. ¹⁶ achieved efficient registration of bridge structures by employing a hash table-based RANSAC for an accelerated candidate search.

However, these existing approaches are primarily designed for static rigid objects and do not explicitly model the kinematic constraints of articulated or branched-chain industrial equipment. In contrast, our proposed method focuses on robust articulated model fitting by integrating explicit kinematic constraints with an edge-weighted registration framework, specifically addressing the challenge posed by nonmodeled components in real factory environments.

3.1. URDF of Articulated Robot

The URDF does not natively support closed-loop mechanisms. To represent articulated robots with such structures, we decompose the model into a main chain and one or more sub-chains, as shown in Fig. 3. This decomposition allows us to maintain a tree structure while modeling complex mechanisms.

We further modularized the model by separating the end effector (EEF) into its own URDF file. As EEFs vary in shape and function depending on the task, this separation allows for flexible reconfiguration. In factory environments, diverse robots can be constructed by combining a limited set of main-link structures with different EEFs. Reusing URDFs for the main chain, sub-chains, and EEFs simplifies the model construction, maintenance, and adaptation. Fig. 4(a) highlights this modular structure, in which the EEF is shown in red.

3.2. URDF of Conveyance System

These structured CAD models provide a practical basis for the development of branched-chain mechanisms with prismatic joints. The conveyance system has a simpler structure without closed loops. It consists solely of a main chain and can be represented using a single URDF. Fig. 4(b) shows the 3D assembly model of the conveyance system, demonstrating how the URDF encodes these branched-chain structures and their spatial relationships for 3D model fitting.

4.1. Normal Estimation

The efficiency of the proposed method stems from leveraging the data structure inherent in TLS. A TLS captures points using a polar coordinate system based on azimuth (\(\theta\)) and elevation (\(\phi\)) angles, allowing the entire point cloud to be treated as a 2D array, \({P[i][j]}\) as shown in Fig. 5(a). This structure permits highly efficient neighborhood searches by accessing neighborhood pixels (\(P{[i\pm\Delta\theta]}{[j\pm\Delta\phi]}\)), which avoids computationally expensive 3D k-d tree searches. For each point \(\boldsymbol{p}_{i,j}\), we define an adaptive neighborhood on the 2D grid ^21,22,23. The pixel radius (\(\Delta\theta\), \(\Delta\phi\)) scales inversely with the point’s distance \(d_{i,j}\) from the scanner. This ensures that the search window covers a consistent physical area and makes the estimation robust against density variations. Normal was computed in a two-stage process using weighted principal component analysis (WPCA) ^24,25. First, the initial normal \(\boldsymbol{n}_{i,j}^{(0)}\) is computed. This WPCA assigns a greater influence to points \(\boldsymbol{p}_k\) closer to the center \(\boldsymbol{p}_{i,j}\), thereby enhancing robustness against noise. This is achieved using the following Gaussian weighting function, \(w_{\textrm{pca}}\), where \({h_{i,j}}^2\) is derived from the neighborhood’s spatial extent:

The parameter \(h_{i,j}\) is adaptively determined for each point as half the maximum Euclidean distance to its neighboring points within the local search window. To preserve sharp edges, only neighbors \(\boldsymbol{p}_k\) whose initial normals \(\boldsymbol{n}_k^{(0)}\) are sufficiently aligned with the central point’s initial normal \(\boldsymbol{n}_{i,j}^{(0)}\) are retained.

The alignment threshold \(\tau_{\textrm{align}}\) is determined empirically according to the expected noise level of the laser scanner to distinguish between noise and geometric discontinuities. A second WPCA is then performed on this filtered set to compute the final normal \(\boldsymbol{n}_{i,j}\), which preserves the sharp edges and prevents smoothing across discontinuities.

By leveraging the 2D grid, applying WPCA, and performing this two-step refinement, our method efficiently computed high-quality normals for both smooth surfaces and sharp edges, which are crucial for subsequent edge detection.

4.2. Edge Candidate Scoring

With accurate normals computed, the potential edge points are identified. Instead of comparing only immediately adjacent neighbors, we evaluated the normal consistency within an annular region around each point \(\boldsymbol{p}_i\). Neighborhood \(\mathcal{N}_i\) is defined as all points \(\boldsymbol{p}_k\) whose distance to \(\boldsymbol{p}_i\) satisfies

4.3. Region Growing for Noise Removal

Edge candidates can include noisy points. To extract coherent structural edges, all edge candidates from all scanner positions are aggregated into a kd-tree. The algorithm iterates through candidates sorted by edge score \(S_i\), selecting the first unassigned point \(p_{i}\) as the seed. Neighboring points within a radius are added to the current region only if their normals satisfy

Regions smaller than a minimum threshold \(\tau_{\textrm{min}}\) are discarded as noise. This process continues until all candidates are assigned to either a valid region or are removed.

Finally, edge weight \(w_{\textit{edge},i}\) is assigned to each point \(\boldsymbol{p}_i\). Points belonging to valid regions that remain after the noise removal process (exceeding \(\tau_{\min}\)) are assigned a high weight \(w_{\textit{high}}\), while all other points are assigned a base weight of 1:

The value of \(w_{\textit{high}}\) was determined through sensitivity analysis by varying the parameter from 2 to 100. The minimum pose error is obtained at \(w_{\textit{high}}=10\). Smaller values reduce the influence of the edge constraints, whereas excessively large values increase the sensitivity to noise at the detected edges. Therefore, \(w_{\textit{high}}=10\) was adopted as the default value in this study, although users can adjust this parameter when applying the method to different datasets.

5.1. Base Link Fitting

The base link defines the origin and coordinate systems of the entire robot. Accurate 6-DoF alignment (XYZ translation and RPY rotation) is essential because errors propagate to all subsequent links. The fitting process was divided into rough and precise registration phases.

5.2. Main Link Fitting

Once the base link (\(L_1\)) is fixed, the pose of each subsequent main link (\(L_k\), for \(k\ge 2\)) is determined sequentially. This process is not a 6-DoF registration but rather an estimation of the single 1-DoF joint angle (\(\theta_k\)) around the kinematically constrained axis defined in the URDF file. The sequential fitting process is shown in Fig. 9(a). This estimation was also performed at the rough and precise stages.

5.2.2. Precise Registration (Rotational)

The angle \(\theta_{k}\) obtained from the rough registration (Section 5.2.1) is used as the initial pose for the final refinement. We perform an axis-constrained Weighted ICP that optimizes the transformation (\(\boldsymbol{R}\), \(\boldsymbol{t}\)) by minimizing the following objective:

\begin{align} E_{k} &= \sum_{i=1}^{N} w_{\textit{edge},i} \left\Vert\boldsymbol{R}\boldsymbol{p}_{i} + \boldsymbol{t} - \boldsymbol{q}_{i}\right\Vert^2 \notag \\ &\phantom{=~} + \lambda\left(\left\Vert\boldsymbol{Ra} + \boldsymbol{t} - \boldsymbol{a}_{0} \right\Vert^{2} + \left\Vert\boldsymbol{Rb} + \boldsymbol{t} - \boldsymbol{b}_{0}\right\Vert^{2}\right). \label{eq:10} \end{align}

The first term is the edge-weighted point-to-point error, where \(w_{\textit{edge},i}\) is the edge weight derived from Eq. \(\eqref{eq:6}\). The second term enforces the kinematic constraint by penalizing the deviations between the transformed axis points of the child link \((\boldsymbol{a},\boldsymbol{b})\) and the corresponding axis points of the parent link \((\boldsymbol{a}_0,\boldsymbol{b}_0)\).

Following the optimization, the resulting rotation matrix \(\boldsymbol{R}\) must be converted back into the single scalar joint parameter \(\theta_{k}\). The final angle \(\theta_k\) is determined by projecting the optimized rotation \(\boldsymbol{R}\) onto the kinematically defined joint axis. This projection ensures that the final \(\theta_k\) strictly adheres to the 1-DoF rotational constraint, eliminating any residual numerical errors in pitch, yaw, or translation that may have persisted from the weighted least-squares solution.

The EEF also consisted of multiple links, as shown in Fig. 10. Because the EEF is connected to the end of the main link, its position and orientation are determined in the same manner based on the axial relationship with the previous link.

5.3. Sub Link Fitting

The rotation angles of the sub-links can be uniquely determined from the rotation angles of the main links. The sub-link positions are computed as inverse kinematic solutions using the fixed main links.

6.1. Base Link Fitting

The conveyance system is equipped with multiple link mechanisms on top of the base link. To determine the position of the base link, feature points were extracted from both the TLS point clouds and the CAD model. Alignment was then performed by matching the corresponding feature points.

The conveyance system has upward-extending arms that hold the automobile parts in place. Therefore, the highest point of each arm is detected as a feature point.

To detect the feature points, a height map was generated by projecting the point cloud onto a horizontal grid, as shown in Fig. 11(a). Pixel value \(P_{(i,j)}\) represents the Z-coordinate (height) of the highest point in each grid cell. The feature points that were significantly higher than those of the surrounding pixels were selected.

Let the pixel value at \((i,j)\) be denoted as \(P_{(i,j)}\). The pixel values in the \(i\)-direction can be approximated by a second-order curve using the Taylor expansion, as shown in Eq. \(\eqref{eq:11}\).

Here, a larger value of \({\Delta^2 P}/{{\Delta i}^2}\) indicates a sharper convex shape, suggesting that the pixel corresponds to the vertex at the tip of the arm. The second-order difference between the pixels is computed as follows:

The same calculation was performed in the \(j\)-direction and along two diagonal directions. Using these results, the sharp convex vertices corresponding to the arm tip are detected as pixels that satisfy the following condition, where \(T\) is the negative threshold for convex detection:

The detected feature points are indicated in blue in Fig. 11(b). Next, the center axis of the conveyance system was determined based on these feature points. The central axis is computed to maximize the symmetry of the feature points. When a pair of symmetrical feature points was selected, the axis was set perpendicular to the line connecting them. To estimate the symmetry axis, the RANSAC method was used to determine the axis that maximized the number of symmetrical feature points. The detected center axis is shown in Fig. 11(b). Feature points that did not exhibit symmetry with respect to the central axis were excluded.

Feature points were extracted from the CAD model as vertices with large Z-coordinates. In the rough fitting process, the distances between pairs of symmetrical feature points were calculated, and the corresponding pairs were identified between the CAD model and point cloud. Once the corresponding pairs were determined, a rough fitting was performed by aligning them, as shown in Fig. 11(c).

Through rough fitting of the CAD model, the points around the base link can be extracted from the point cloud. Subsequently, an accurate fitting was achieved by applying a weighted ICP (Section 5.1.2) to align the CAD point clouds with the TLS point clouds. Here, the edge weight \(w_{\textit{edge},i}\) is the weight derived from the edge detection results in Eq. \(\eqref{eq:6}\).

6.2. Main Link Fitting

Once the base link was fixed, the pose of each subsequent main link was determined. The sequential fitting process and geometric constraints for the branched-chain mechanism are shown in Fig. 12. In conveyance systems, these links typically move along a 1-DoF prismatic (translational) axis, as defined in the URDF. As in the case of the articulated robot, the fitting process was performed at the rough and precise stages.

7.1. Qualitative Evaluation on Real Factory Point Clouds

The proposed method was applied to point clouds captured at an automotive factory. The input point clouds for the articulated robot and conveyance system contained approximately 5.5 million and 5.0 million points, respectively. The average point spacing (grid resolution), calculated as the mean distance to the nearest neighbor across all valid points, was approximately 1.66 mm for the articulated robot and 1.60 mm for the conveyance system. The OBB size was approximately \(\mbox{2.5~m} \times \mbox{2.0~m} \times \mbox{2.5~m}\) for the articulated robot, and \(\mbox{1.6~m} \times \mbox{1.7~m} \times \mbox{5.5~m}\) for the conveyance system. These datasets contain significant sensor noise, occlusions, and numerous nonmodeled accessories such as wire harnesses and covers, which act as geometric outliers. Despite these challenges, the method successfully aligned the CAD models of both articulated robots and conveyance systems, as shown in Fig. 13. By assigning high weights to stable geometric edges, the influence of nonmodeled components was effectively suppressed, allowing for accurate reconstruction of kinematically valid 3D assembly models.

To evaluate the effect of the edge-weighting strategy, we compared our approach with a conventional unweighted registration method using point clouds acquired from real robots (Fig. 14). In the unweighted case, small fitting errors in individual links accumulate along the kinematic chain, leading to a significant misalignment of the robot’s EEF. By contrast, the edge-weighted method anchors each registration step to reliable geometric features, effectively minimizing the cumulative error and ensuring precise alignment, even at the extremities of the machine.

7.2. Computational Efficiency

The efficiency of the proposed method was evaluated using a desktop PC equipped with an Intel Core i9-10900K CPU @ 3.70 GHz and 62.7 GiB RAM. The processing times for both target systems are summarized as follows:

• Articulated robot: Edge detection took 9.9 s, base-link registration required 296.7 s, and sequential fitting (main/sub-chains and EEF) was completed in 12.6 s.

• Conveyance system: Edge detection took 10.0 s, base-link registration required 143.0 s, and sequential fitting was finished in 37.3 s.

Base-link registration is the most time-consuming stage because it involves a RANSAC-based global search and feature matching to establish an initial coordinate system. However, the subsequent link-by-link pose estimation is significantly accelerated by leveraging the structured 2D grid of the TLS data for efficient neighbor searching.

These results demonstrated that the proposed method is computationally practical for processing large-scale point clouds in factory environments.

7.3. Quantitative Evaluation Setup

We assessed the robustness of the method using virtual point clouds generated from CAD models. The virtual datasets for the articulated robot and conveyance system consisted of approximately 4.7 million and 7.3 million points, respectively. The average point spacing was approximately 1.7 mm for the articulated robot and 2.36 mm for the conveyance system, respectively. This evaluation used three unique postures for both the articulated robot and the conveyance system. For each instance, we generated two datasets: one without and one with nonmodeled accessories (e.g., wire harnesses and covers). Fig. 15 shows the virtual CAD models. To simulate real acquisition artifacts, we nonuniformly placed virtual laser scanners to induce occlusions and applied Gaussian noise (\(\sigma=1\) mm) along the laser irradiation direction to the resulting point clouds. Fig. 16 illustrates the virtual data acquisition setup and the resulting noisy point cloud details.

7.4. Quantitative Results and Discussion

The errors in the rotational angles of the articulated robot and translational distances of the conveyance system are summarized in Tables 1 and 2, respectively. Under the “With Accessories” condition, the proposed edge-weighted precise fitting consistently achieved the lowest average error compared with the unweighted approach. These results demonstrate that assigning high weights to structural edges effectively suppresses interference from nonmodeled components, thereby significantly improving the robustness of posture estimation for both single-chain and branched-chain link mechanisms.

Although the observed accuracy improvements might appear marginal, they are critical for large-scale industrial systems where minor pose errors are amplified along extended kinematic chains. For the articulated robot in this study (overall height \(\approx 2.5\) m), the EEF extends approximately 1.0 m from the final joint to the welding tip. An angular error of merely 0.1° at the joint level thus results in a positional deviation of \(L\times \theta\approx 1.75\) mm at the tooltip. Such joint-level inaccuracies accumulate millimeter-scale endpoint errors, which are significant in precision-critical manufacturing. Similarly, in the 5.5 m long conveyance system, minor axis misalignments led to amplified deviations at the tips of the arm segments. Consequently, submillimeter pose accuracy is indispensable for high-fidelity virtual environments to ensure reliable collision detection and safe motion planning.