Automatic image alignment and stitching for
ultrasound-based robotic inspection of complex
geometry components

Ekaterina Iakovleva, David Roue, Philippe Bredif

► To cite this version:

Ekaterina Iakovleva, David Roue, Philippe Bredif. Automatic image alignment and stitching for ultrasound-based robotic inspection of complex geometry components. NDT & E International, 2024, 146, pp.103158. 10.1016/j.ndteint.2024.103158. cea-04982779

HAL Id: cea-04982779

https://cea.hal.science/cea-04982779v1

Submitted on 7 Mar 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic diagram of TFM imaging geometry.</figcaption></figure>

scan matching techniques [7,8]. The most common solution for scan matching in the robotics community is the ICP algorithm. This algorithm allows finding the rigid transformation between two overlapping point sets. It is based on an iterative process where it first establishes correspondences between two sets of points and then estimates a transformation that minimizes the distance between the corresponding points. This process is repeated with the new estimate until the convergence. The most often cited analysis of the algorithm comes from Besl and McKay [9]. The authors show that it achieves registration with good accuracy and terminates in a minimum number of iterations. Since the introduction of ICP, many variants have been introduced, and a brief overview of which can be found [10].

This paper presents a new stitching method that involves pairwise image registration. To mitigate the surface distortion in the stitched image caused by vibrations at high scanning speeds, two consecutive ultrasound images with overlapping areas are registered using the most popular ICP variant, also called point-to-point ICP [9]. This algorithm uses a point-to-point metric (Euclidean distance between corresponding points) and an exact closed-form solution for minimizing such metric. In this way, all image pairs are successively registered relative to each other to form a final image that provides an accurate surface approximation.

In this work, we focus on 2D TFM imaging using a 1D linear array, assuming that the probe moves parallel to the line formed by the array elements. However, we note that the proposed method can equally be applied to other imaging algorithms (delay and sum beamforming, plane wave imaging) and can be easily extended to 2D arrays and 3D imaging, since most versions of ICP are employed to match 2D and 3D point sets. We note also that our approach can be applied to form an image of the interior of the component using, for example, backwall echoes to perform registration.

This paper proceeds by reviewing TFM and ICP algorithms, and then introducing the new TFM stitching procedure using the ICP registration algorithm. Then our approach is compared to a reference model using synthetic data. Finally, the results obtained on experimental data are presented, demonstrating the advantages of the new stitching method.

## 2. Total Focusing Method

### 2.1. Principle

In this section, we recall the principle of conventional TFM, which will be used later to produce an image of a surface profile. TFM processing algorithm uses Full Matrix Capture (FMC) data, which is the set of signals collected by all transmitter-receiver pairs of an ultrasonic array to produce an image [11,12]. The intensity of the output image at pixel $p(x, z)$ is obtained by summing all received signals shifted in time:

$$I(p) = \left| \sum_{i,j=1}^{N_e} S_{ij}(\tau(x_i, x_j, p)) \right| \quad (1)$$

where $S_{ij}(\tau) = s_{ij}(\tau) + jH(s_{ij}(\tau))$ is the analytic representation of the signal $s_{ij}(\tau)$, received by element $j$ when emitting with element $i$, $H$ denotes the Hilbert transform, $N_e$ is the number of elements of the array and $\tau(x_i, x_j, p)$ is the corresponding time of flight given by:

$$\tau(x_i, x_j, p) = \frac{\sqrt{(x-x_i)^2 + z^2 + \sqrt{(x-x_j)^2 + z^2}}}{c} \quad (2)$$

where $x_i$ and $x_j$ are the lateral positions of transmitting and receiving elements, respectively, and $c$ denotes the sound velocity in the propagation medium. Fig. 1 schematically depicts the image formation process.

### 2.2. Grating lobe suppression in TFM imaging

In TFM imaging, the pitch of elements (center-to-center distance of elements) ideally should be less than one half of the operating wavelength [13–15]. When the sampling criterion is not satisfied, the TFM process leads to the formation of artifacts in the image, generally referred to as grating lobes.

In order to avoid their formation in the image, the method based on resampling the FMC data before TFM processing can be employed. This approach is widely used in imaging methods where imaging is performed in the frequency domain, in which the data with sufficient spatial sampling is reconstructed from the spatially undersampled initial data using zero-padding (or sinc interpolation) [16–18]. In the time domain, the resampling process can be performed with linear methods such as nearest neighbor, bilinear, bicubic, Lanczos or sinc interpolation [19]. In this case, the array aperture is synthesized with the desired pitch, $dx$, by subdivision of each interval $[x_i, x_{i+1}]$ into $\alpha > 1$ equal parts with respect to $dx/\alpha \le \lambda/2$, where $x_i$ denotes the position of the $i^{th}$ element along the x-axis. The signals corresponding to the new transmit/receive sub-elements $\left\{ \left(x_i + k \frac{dx}{\alpha}, x_j + l \frac{dx}{\alpha}\right) \right\}_{i,j=1,...,N_e-1}$ are approximated by a two-dimensional interpolation of the received signals $s_{ij}(\tau)$ at neighboring grid points $\{ (x_i, x_j) \}_{i,j=1,...,N_e}$. For example, in the case of bilinear interpolation, the interpolated signal corresponding to the new transmit/receive pair $\left(x_i + k \frac{dx}{\alpha}, x_j + l \frac{dx}{\alpha}\right)$ can then be written as:

$$s_{ijkl}(\tau) = \left(1 - \frac{k}{\alpha}\right) \left(1 - \frac{l}{\alpha}\right) s_{ij}(\tau) + \frac{k}{\alpha} \left(1 - \frac{l}{\alpha}\right) s_{i+1,j}(\tau) + \left(1 - \frac{k}{\alpha}\right) \frac{l}{\alpha} s_{i,j+1}(\tau) \\ + \frac{k l}{\alpha \alpha} s_{i+1,j+1}(\tau) \quad (3)$$

The time of flight for TFM imaging is then calculated by Eq. (2) on the new transmit/receive grid with virtual pitch $dx/\alpha$. Finally, the TFM image is obtained by combining the TFM images computed using the received and interpolated signals:

$$I_{\alpha}(p) = \left| \sum_{i,j=1}^{N_e} S_{ij}(\tau(x_i, x_j, p)) + \sum_{i,j=1}^{N_e-1} \sum_{k=1}^{a-1} S_{ijkl} \left( \tau \left( x_i + k \frac{dx}{\alpha}, x_j + l \frac{dx}{\alpha}, p \right) \right) \right| \quad (4)$$

where $S_{ijkl}(\tau)$ is the analytic representation of the interpolated signal $s_{ijkl}(\tau)$.

The similar idea applied to conventional delay-and-sum beamforming on receive using nearest neighbor approach can be found [20]. Although the resampling procedure increases the computational burden,
<figure><img src="image_3.png" /><figcaption>Fig. 2. The registration of ICP: (a) definition of correspondences between two point sets by searching closest points; (b) registration using rigid transformation.</figcaption></figure>

since it increases the volume of data to be processed, as we will see later, it allows us to work with larger separation $dx$, maintaining good contrast and resolution.

### 3. Iterative Closest Point algorithm

In this section, we introduce the classic ICP algorithm to perform 2D point-set registration. This method will be used in our stitching procedure to estimate the transformation parameters needed to align two overlapping TFM images into a common coordinate system.

Given two sets of points $Q = \{q_j\}_{j=1,...,M}$ and $P = \{p_i\}_{i=1,...,N}$ in $\mathbb{R}^2$, the objective is to estimate a 2D rigid transformation, represented by an orthonormal rotation matrix $R \in \mathbb{R}^{2 \times 2}$ and a translation vector $t \in \mathbb{R}^2$, with which $P$ is registered to be in the best alignment with the reference point set $Q$ (Fig. 2). The optimal transformation between two point sets can be found by minimizing the sum of squared errors:

$$ (R, t) = \arg \min_{R,t} \left( \sum_{i=1}^{N} \omega_i \|Rp_i + t - m_i\|^2 \right) \qquad (5) $$

where $m_i$ is the closest point in $Q$ to the point $p_i$, and $\omega_i$ is a weight for corresponding pair $(p_i, m_i)$.

Given an initial guess, $(R_0, t_0)$, the ICP algorithm solves the problem Eq. (5) using an iterative process. At each iteration $k \ge 0$, the following steps are performed:

1. Compute the coordinates of each point $p_i$ in $P$ in the reference frame, $Q$, using the current guess $(R_k, t_k)$:

$$ p_i^k = R_k p_i + t_k $$

2. Compute the set of correspondences between two point sets by:

$$ m_i = \arg \min_{q \in Q} (\|p_i^k - q\|) $$

3. Compute the weights for a point match using matching threshold, $d_{max}$:

<figure><img src="image_4.png" /><figcaption>Fig. 3. An example of TFM image stitched from two overlapping images of the same specimen profile using registration.</figcaption></figure>

$$ \omega_i = \begin{cases} 1, & \text{if } \|p_i^k - m_i\| \le d_{max} \\ 0, & \text{otherwise} \end{cases} $$

4. Compute the new transformation by minimizing squared distance Eq. (5) between the corresponding points:

$$ (\hat{R}, \hat{t}) = \arg \min_{R,t} \left( \sum_{i=1}^{N} \omega_i \|Rp_i^k + t - m_i\|^2 \right) $$

5. Update the complete transformation:

$$ R_{k+1} = \hat{R} R_k, t_{k+1} = \hat{R} t_k + \hat{t} $$

These steps are repeated until convergence (if the mean squared error Eq. (5) is within a given tolerance) or until the maximum number of iterations is reached. At each iteration, the fourth step of the ICP
<figure><img src="image_5.png" /><figcaption>Fig. 4. TFM imaging geometry used for simulation.</figcaption></figure>

algorithm can be solved by the Singular Value Decomposition (SVD) [21], where the calculation of the optimal transformation in Eq. (5) is based on the SVD of the cross-covariance matrix of two sets of points defined by:

$$H = \sum_{i=1}^{N} \omega_i (p_i - p_G)(m_i - m_G)^T$$

Here $p_G = \sum_{i=1}^{N} \omega_i p_i$ and $m_G = \sum_{i=1}^{N} \omega_i m_i$ represent the center of mass of both point sets. Therefore, the optimal rotation matrix is expressed as follows:

$$\hat{R} = VU^T \qquad (6)$$

where U are the left singular vectors and V are the right singular vectors of the cross-covariance matrix $H = U\Sigma V^T$. Moreover, the optimal translation vector is calculated by the following equation:

$$\hat{t} = m_G - \hat{R}p_G \qquad (7)$$

Note that the optimal transformation in Eq. (5) can also be calculated by quaternion-based methods [22,23] or by an algorithm using orthogonal matrices [24]. All these algorithms demonstrate the same performance and stability against noisy data [25].

The assignment of weights in the third step of the algorithm is related to the rejection of certain corresponding pairs in the least squares minimization Eq. (5). The pair $(p_i, m_i)$ is rejected if the distance between the point $p_i$ and its closest point $m_i$ is above the specified matching threshold, $d_{max}$. The purpose of this is usually to eliminate outliers and efficiently handle partially overlapping point sets, as classical ICP assumes a full overlap of the point sets being matched [26].

Finally, the second step of the algorithm, which is the most computationally expensive, especially for large datasets, can be solved using acceleration technique based on the k-d tree [27–29]. In this case, for each iteration, the two nearest neighbor points search can be performed in $\mathcal{C}(N \log(M))$ time instead of $\mathcal{C}(NM)$ time required for linear search.

## 4. TFM image stitching

In this section, we present the image stitching procedure designed to create a single 2D TFM image by combining multiple overlapping 2D TFM images (local images) obtained for different probe positions. As shown in Section 2, the local image is formed in the probe coordinate system. Therefore, a natural approach to stitch all images together is to choose the image corresponding to the first probe position as the reference and to then warp all of the other images into the reference coordinate system. An example of stitching two images is shown in Fig. 3.

We suppose that the probe position is defined in the xz-plane by a two-dimensional coordinate, $o \in \mathbb{R}^2$, and a rotation, $R \in \mathbb{R}^{2 \times 2}$, relative to

some coordinate system. Then, the transformation between two consecutive probe positions $k-1$ and $k$ can be calculated as follows:

$$q = R_{k-1,k}p + t_{k-1,k} \qquad (8)$$

where $R_{k-1,k} = R_{k-1}^{-1}R_k$ is the rotation matrix and $t_{k-1,k} = R_{k-1}^{-1}(o_k - o_{k-1})$ is the translation vector, $p$ is any point in $\mathbb{R}^2$ and $q$ is the transformed point.

Let us denote K images by $\mathcal{J} = \{I_k\}_{k=1,...,K}$ with the corresponding initial probe position estimates $\mathcal{T} = \{(R_k, t_k)\}_{k=1,...,K}$ provided by the robotic system.

Then the transformation $(R_{k-1,k}, t_{k-1,k})$ between each overlapping image pair $I_{k-1}$ and $I_k$ can be calculated by Eq. (8). This transformation allows to represent each point of the local image $I_k$ in the coordinate system of the local image $I_{k-1}$ and will be used as initial estimate in the iterative process described in the previous section.

Using these notations, our image stitching algorithm can be summarized in the following steps.

1. Detect keypoints (points of interest) for two consecutive images $I_{k-1}$ and $I_k$. In our application case, the keypoints are the surface points, which are detected by segmentation of TFM image using Otsu's thresholding method [30]. This method does not require input parameters. The optimal threshold value is determined by separating the image histogram into two classes (foreground and background) with maximum inter-class variance. The threshold is then applied to all image pixels in the z direction to retain the surface points with the highest intensity. At this stage, two sets of points, $Q_{k-1}$ and $P_k$, are created.

2. Compute the optimal transformation $(\hat{R}_{k-1,k}, \hat{t}_{k-1,k})$ such that $P_k$ is best aligned to $Q_{k-1}$ using ICP algorithm with an initial estimate provided by the relative displacement of the current probe given by Eq. (8).

3. Warp the current image $I_k$ into the reference coordinate system (first image, $I_1$) using the transformation obtained from Step 2 and all previously obtained motion estimates. The corresponding transformation is calculated by multiplying all pairwise transformations. Image warping can be performed using linear methods [19].

Repeat these steps for all images in $\mathcal{J}$, starting with the second image. The final composite image is then created by taking an average value at each pixel:

$$C(p) = \sum_{k=1}^{K} \omega_k(p) \tilde{I}_k(p) / \sum_{k=1}^{K} \omega_k(p) \qquad (9)$$

where $\tilde{I}_k(p)$ are the warped images and $\omega_k(p)$ is 1 at valid pixels and 0 elsewhere.
<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Conventional TFM; (b)–(f) TFM image $I_α(p)$ generated from FMC data oversampled by various factors $α$ using nearest neighbor interpolation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. (a) Conventional TFM; (b)–(f) TFM image $I_α(p)$ generated from FMC data oversampled by various factors $α$ using bilinear interpolation.</figcaption></figure>

## 5. Results with synthetic data

The stitching algorithm described above was tested on a set of synthetically-generated data. The motivation for using synthetic data in our study is that we know the exact transformation between two consecutive probe positions, and can evaluate the performance of the proposed algorithm relative to this alignment. For this purpose, we implemented two stitching methods, the proposed algorithm and the stitching procedure without registration step. In the latter case, the image is transformed into the reference coordinate system using initial probe positions Eq. (8). Both algorithms were implemented in Python and run on an Intel Core i9 2.4-GHz CPU with 32 GB of RAM. In our implementation the neighbors package of the sklearn module is used for two nearest neighbor points search (the second step of the ICP algorithm).

The FMC data were simulated with the CIVA software package ([31, 32]). The probe is a 1D array of 16 elements, the element pitch is 2.1 mm, the gap between two adjacent elements is 0.3 mm and the total aperture of the array is 33.3 mm. The elements are excited by a Gaussian signal with the central frequency of 3 MHz and a bandwidth of 50 % as -6 dB. The test specimen employed in the imaging is a steel block of 256 mm × 80 mm with relatively smooth irregular geometry.

Fig. 4 details schematic diagram of the cross-section of the test specimen. The distance between the array and the flat part of the surface

is 40 mm with water acting as a coupling medium. The sound velocity in the water is 1480 mm/s, such that the operating wavelength, $λ$, is around 0.5 mm. The local reconstruction was performed in the probe coordinate system. The image, centered on [0, 40] mm in the xz-plane (with respect to the probe centerline), has dimensions 10 mm × 35 mm, with a pixel size of 0.12 mm (≈$λ$/4 in water).

### 5.1. Local TFM imaging

Before discussing the results obtained with our stitching method, we first examine the local TFM reconstruction, since its quality directly affects the convergence to a correct alignment during the stitching process. Fig. 5(a) and Fig. 6(a) show the TFM reconstruction defined by Eq. (1) obtained for the probe positioned over a junction formed by two flat surfaces with different inclination angles (the region enclosed by the dashed line in Fig. 4). In this case, the image under test has dimensions 20 mm × 35 mm and the resolution of the imaging grid is 0.12 mm (≈$λ$/4). It can be seen that the TFM image exhibits high-intensity grating-lobe artifacts that make it difficult to extract the surface profile, as they can be attributed to true surface. One factor determining their appearance is the probe pitch, which in our experiment is 4λ. It is chosen with respect to the longitudinal wavelength of ultrasound in the test specimen material, which is 2 mm for the acoustic velocity of 5950 mm/s. To control the undesired lobes we use the resampling technique
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Stitched TFM images generated from the first dataset (with the addition of random noise to the initial probe positions within ±1 mm along the z-axis): (a) stitching along a linear probe trajectory; (b) stitching using the ICP algorithm; (c) stitching using a true probe trajectory (reference).</figcaption></figure>

described in Section 2.2. Here we present results of resampling with two simple and fast interpolation algorithms: nearest neighbor and bilinear interpolation. Fig. 5 (b)–(f) and Fig. 6(b)–(f) display the reconstructions $I_\alpha(p)$ defined by Eq. (4) for FMC data oversampled by a factor $\alpha = 2, 3, 4, 6$ and 8 with nearest neighbor and bilinear interpolation, respectively. As can be seen, both interpolation algorithms significantly improve TFM image quality. The bilinear interpolation is more complex than the nearest neighbor method but performs better due to its low-pass filtering properties. For bilinear interpolation function, the results show that already for $\lambda$ elements spacing (factor 4), the grating lobes are completely removed, so the surface profile can be accurately reconstructed using global thresholding method mentioned in Section 4, making the surface profile extraction procedure fully automatic. It can also be seen that the resampling method does not change the resolution in the TFM image, while it considerably improves the contrast. The main drawback of the proposed method is that it increases the computational burden, as the volume of data to be processed increases. To accelerate computation times, the surface imaging can be carried out in the frequency-wavenumber domain using spatial interpolation through zero-padding for image enhancement. However, these low-cost imaging techniques can only be applied to simple configurations. The approach

developed in the time domain can be used to image both the surface and the interior of the specimen with different scale factors.

In the rest of the study the resampling technique with bilinear interpolation and a scale factor 4 is employed for reducing grating lobes in the TFM image. The bilinear interpolation in resampling is used because it is easy to implement, efficient and, as the results show, has sufficiently good performance.

## 5.2. Stitching results

In this subsection, we examine the stitching results. Note that for the considered sample, which includes two flat sections at both ends (see Fig. 4), the ICP method cannot detect the horizontal movements of the probe over these sections. This drawback comes from the difficulty to find a correct match between two partially overlapping parallel flat surfaces with nearest neighbor search. Therefore, in this work, we have limited ourselves to measuring only the vertical movements and rotations of the probe.

To create a tiled TFM image of the whole specimen, the array was moved 45 times along the x-axis with a step of 5 mm, and therefore, the overlapping area between two consecutive local images is around 85 %

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Stitched TFM images generated from the second dataset (with the addition of random noise to the initial probe positions within ±1 mm along the z-axis and within ±1° in incidence): (a) stitching along a linear probe trajectory; (b) stitching using the ICP algorithm; (c) stitching using a true probe trajectory (reference).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Stitched TFM images generated from the third dataset (with the addition of random noise to the initial probe positions within ±1 mm along the z-axis and within ±3° in incidence): (a) stitching along a linear probe trajectory; (b) stitching using the ICP algorithm; (c) stitching using a true probe trajectory (reference).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Stitched TFM images generated from the fourth dataset (with the addition of random noise to the initial probe positions within ±1 mm along the z-axis and within ±4° in incidence): (a) stitching along a linear probe trajectory; (b) stitching using the ICP algorithm; (c) stitching using a true probe trajectory (reference).</figcaption></figure>

of the total array aperture. To simulate robot vibrations, random noise was added to the linear probe trajectory in order to disturb the vertical position and incidence angle of the probe in the xz-plane. Finally, the global image was stitched from all simulated local images in three ways.

a) The first image was stitched along the linear probe trajectory (nominal trajectory) without registration using composition process Eq. (9). In this case, the vertical probe centerline is considered parallel to the z-axis, so the incidence angle of the probe is 0°, and moves only along the x-axis with a regular step.

b) The second stitching was performed by the method described in the previous section. The nominal trajectory is used as an initial estimate in the registration step of the ICP algorithm.

c) The third image, which serves as the reference image here, was stitched without ICP registration along the perturbed probe trajectory used to simulate the robot vibrations (true trajectory). In this case, all individual images were transformed into the reference coordinate system by Eq. (8) using the perturbed probe positions and then composed using formula Eq. (9).

To measure the performance of the stitching algorithms, these three

stitched images were formed from each of the flowing datasets. Four TFM tiled scans were simulated with a uniformly generated error term added to the linear probe trajectory. For the first dataset (first tiled scan), the error term was set within ±1 mm only along the z-axis. For the next three datasets, we also added the error in incidence (the angle between the vertical probe centerline and the z-axis) within the range of ±1°, ±3° and ±4° for the second, third and fourth dataset, respectively.

For each dataset, the performances were measured by comparing the surface profiles extracted from the stitched images with the reference profile (true profile). The surface points were extracted using Otsu's thresholding method described in the previous section (the first step of the stitching algorithm) and then resampled to a regular grid with a sampling step of 1 mm (≈2λ) using linear interpolation. In our implementation, the ICP algorithm was limited to a maximum of 50 iterations and the value of the maximum matching threshold $d_{max}$ was fixed at 2 mm for the first two datasets and at 4 mm for the third and fourth datasets.

We now examine the stitching results. Fig. 7, Fig. 8, Figs. 9 and 10 show the results obtained for the first, second, third and fourth datasets, respectively, defined above. It can be easily seen that the quality of the
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Surface measurement and error estimation: (a) the reference profile (true profile) and four profiles extracted from stitched images obtained using the ICP algorithm; (b) the absolute error between the reference and measured profiles.</figcaption></figure>

stitched images degrades considerably when the exact positions of the probe are unknown, making them unsuitable for analysis (Figs. 7(a)–8 (a) and 9(a) and Fig. 10(a)). However, the results obtained with stitching using local image alignment show a significant improvement in surface image quality in terms of resolution (reduction of surface echoes splitting) and SNR (Signal-to-Noise Ratio) (Figs. 7(b)–8(b) and 9(b) and Fig. 10(b)).

Fig. 11(a) compares the reference profile and the profiles extracted from the stitched images formed using the ICP registration algorithm. The absolute errors found in comparison with the reference profile and measured profiles are presented in Fig. 11(b). The maximum absolute errors are 0.09 mm (≈λ/5), 0.4 mm (≈λ), 0.6 mm and 1.53 mm and the RMS errors are 0.02 mm, 0.14 mm, 0.28 mm and 0.73 mm for the first, second, third and fourth datasets, respectively.

It can be seen that, for the first three datasets, the extracted profiles are in good agreement with the reference profile. Despite the good image obtained for the fourth dataset, we can see that the stitching using pairwise registration accumulated registration errors resulting in erroneous alignment. The corresponding profile exhibits high reconstruction maximum (1.53 mm > 3 λ) and RMS (0.76 mm > λ) errors. This result shows that the resolving power of the ICP algorithm decreases as the deviation angle from the initial trajectory increases. Indeed, the good agreement between the profiles is obtained for small deviation angles, less than 1° (for the first two datasets). This is due to the fact that a large deviation angle reduces the overlapping area of two successive images. This increases the chances of incorrect correspondences and moves the final alignment away from the correct value. In general, the accuracy of the ICP method can be improved by setting an appropriate value of the maximum matching threshold $d_{max}$, which controls the correspondences, and/or by setting the initial guess close enough to the solution. In practice, it is quite difficult to find a common matching threshold and an initial guess suitable for all image pairs used in the stitching process. Another way to improve the convergence of ICP with a large deviation angle to the correct solution is to increase the overlapping area between two successive images by increasing the probe aperture and reducing the scanning step. However, this approach increases the acquisition time and the volume of data to be post-processed.

Note also that high reconstruction errors were observed on the convex part of the profile for all measured profiles. The ICP algorithm does not converge to the correct alignment in this region due to a poor initial dataset caused by poor imaging of the convex part.

Finally, in terms of computation time, stitching using the ICP registration algorithm is much slower than stitching along the nominal probe trajectory, which is explained by the iterative nature of the registration process. We found that the ICP algorithm converges quite fast, on average in 5–10 iterations. The stitched image computation time is 0.43

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Robot trajectory relative to the inspected component and array settings (a) and experimental setup (b).</figcaption></figure>

s and 0.07 s for the stitching using image alignment and for the stitching along the nominal probe trajectory, respectively, yielding a computational gain of 6 for the second algorithm. The indicated time does not take into account the time taken to form the local TFM images.

# 6. Results with experimental data

## 6.1. Stitching results

In this subsection, we evaluate our stitching method on a set of experimental data. The experimental setup, shown in Fig. 12, consisted of a matrix array with 256 (16x16) squared elements manufactured
<figure><img src="image_14.png" /><figcaption>Fig. 13. (a) FMC B-Scan image (scanning/time) for the first increment position; (b) Stitched TFM images generated for the first increment position along a linear probe trajectory; (c) Stitched TFM images generated for the first increment position using ICP algorithm.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. (a) FMC B-Scan image (scanning/time) for the second increment position; (b) Stitched TFM images generated for the second increment position along a linear probe trajectory; (c) Stitched TFM images generated for the second increment position using ICP algorithm.</figcaption></figure>

(Imasonic, France) with the same dimensions and ultrasonic properties as used to calculate the simulated data. That is the element pitch is 2.1 mm and the gap between two adjacent elements in two directions is 0.3 mm. To produce a 2D image, the FMC data were collected using only one row of the matrix. Experimental data are acquired using two 128-channel M2M Panther system (Eddyfi Technologies, Canada) and the robot UR10 (Universal Robots, Denmark), which offers six degrees of freedom and the ability to implement complex trajectories. The probe was moved at a scanning speed of 70 mm/s. To better understand the character of robot vibrations, the active row was scanned over a large area in two directions (scanning and increment axis), such that the row was parallel to the scanning axis. The probe was moved 45 times along the scanning axis with the 5 mm step and 26 times along the increment axis with the 3 mm step. Our method was tested on steel block of 256 mm × 80 mm × 80 mm comprising a planar extruded profile. The specimen was immersed in a water tank. The nominal height between the probe and the front of specimen is 40 mm on the flat surface.

As above, the local reconstruction was performed in the probe coordinate system with FMC bilinear resampling and a scale factor 4. The stitched image was formed for each increment position along the scanning axis with nominal trajectory (linear trajectory) provided by the robotic system. In total, 26 stitched TFM images were formed. Figs. 13 and 14 detail a selection of results obtained for two increment positions. The FMC B-Scan images collected during acquisition (Figs. 13(a) and Fig. 14(a)) exhibit vibration-induced noise resulting in surface echo splitting in the corresponding stitched TFM images obtained with the nominal robot trajectory (Figs. 13(b) and Fig. 14(b)). Compared to the stitched images formed by the new method (Figs. 13(c) and Fig. 14(c)), it can be seen that stitching with local image alignment significantly improves the surface image.

The results of profile reconstruction from stitched TFM images are presented in Fig. 15(a) and Fig. 16(a). As in the previous section, the surface points were extracted using Otsu's thresholding method and resampled to a regular grid with a sampling step of 1 mm using linear
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Surface measurement and error estimation for the first increment position: (a) reference profile (black line) and profiles extracted from stitched images obtained with ICP registration algorithm (green line) and without registration (blue line); (b) absolute error between the reference and measured profiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Surface measurement and error estimation for the second increment position: (a) reference profile (black line) and profiles extracted from stitched images obtained with ICP registration algorithm (green line) and without registration (blue line); (b) absolute error between the reference and measured profiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Probe motion estimates for two increment positions: (a) measured rotation angle; (b) measured vertical displacement.</figcaption></figure>

interpolation. The absolute errors found in comparison with the reference profile and measured profiles are presented in Figs. 15(b) and Fig. 16(b).

For profiles reconstructed with the ICP registration algorithm for two increment positions, the maximum absolute errors are 0.2 mm ($<$ λ/2) and 0.18 ($<$ λ/2) mm and the RMS errors are 0.08 mm ($<$ λ/6) and 0.07 mm ($<$ λ/6) for the first and second increment position, respectively. It

can be seen that, both ICP-derived profiles are similar and agree quite well with the reference profile.

The profiles extracted from the stitched images formed without image alignment shows high reconstruction maximum (0.58 mm and 0.63 mm > λ) and RMS (0.18 mm and 0.19 mm > λ/2) errors. As can be seen in Figs. 15(a) and Fig. 16(a), these profiles represent the concave part of the surface less well compared to the ICP-derived profiles.
<figure><img src="image_19.png" /><figcaption>Fig. 18. Interior TFM imaging geometry.</figcaption></figure>

Although they agree well visually with the reference profile, they contain a lot of noise, particularly at the junction of two flat sections. As we will see below, such profiles cannot be used for interior imaging without preliminary smoothing. Note also that, as for the simulated profiles discussed in the previous section, high reconstruction errors were observed in the region containing a convex part of the profile for all measured profiles. This was again caused by the poor imaging of this part.

The probe motion estimates obtained with ICP registration algorithm for the selected data are presented in Fig. 17. The result shows that the robot movements during inspection are affected by vertical vibrations, while the rotational vibrations are quite insignificant (maximum angle of -0.63° is observed for the second increment position). The vertical displacement within ±1.5 mm (>2λ) and the rotation angle within ±0.7° were observed from the results obtained for all collected data. At this stage, it is important to highlight that due to the good superposition of the measured profiles with the reference profile, it can be concluded that the horizontal vibrations of the robot used in the experiment are quite insignificant.

Note that if the mechanical movements of the robot result in only vertical vibrations, the ICP registration algorithm can be easily

simplified by setting the rotation angle and horizontal displacement to zero in the least squares minimization Eq. (5).

In this particular case, the algorithm converges in just one iteration to the vertical translation: $\hat{c}_z = m_{G,z} - p_{G,z}$, where $m_{G,z}$ and $p_{G,z}$ represent the vertical position of the center of mass of two point sets being matched. The correspondences between two point sets can then be computed by minimizing the scalar product: $m_i = \arg \min_{q \in Q} (\|\langle p_i - q, \vec{n} \rangle\|)$, where $\vec{n}$ is the unit vector in the horizontal direction. The overlapping region of two successive images can be calculated directly using the known displacement of the probe in the horizontal direction, so the maximum matching threshold $d_{max}$ can be disregarded in the algorithm. Note also that this simplified variant is much faster than the classical ICP, which is explained by its simplified probe motion estimation scheme.

## 6.2. Interior imaging

We now examine the quality of the measured profiles by comparing the interior TFM images of the sample in the region of strong vibrations of the robot during inspection. The configuration used for interior imaging is shown in Fig. 18. The probe was placed over the junction of two flat sections (the region enclosed by the dashed line in Fig. 4). A sample contains 5 SDH of 1.5 mm diameter at a depth of 18 mm with respect to the flat part parallel to the x-axis, spaced by 10 mm. The image has dimensions 20 mm × 60 mm with a pixel size of 0.1 mm. The FMC signals were simulated with CIVA using the true surface. The propagation time for TFM imaging were computed using beam propagation model implemented in CIVA [33], taking into account the measured profile.

The TFM results obtained with the true and measured profiles can be compared in Fig. 19. It can be seen, that the images corresponding to the reference profile (Fig. 19(a)) and the ICP-derived profiles for two increment positions (Fig. 19(c) and (d)) give similar results. However, a slight loss of reconstruction quality (in terms of location and SNR) can be observed for two SDHs located under the inclined part of the profile for both profiles. This was due to insufficiently accurate reconstruction of the profiles in this region, as can be observed in Figs. 15 and 16. The TFM image corresponding to the profile reconstructed without image alignment obtained for the second increment position is less accurate (Fig. 19(b)). We observe splitting of the SDH echoes and a significant reduction in contrast, which make it unsuitable for analysis.

<figure><img src="image_20.png" /><figcaption>Fig. 19. Interior TFM images using: (a) reference profile; (b) profile obtained without image alignment for the second increment position; (c) profile obtained with ICP for the first increment position; (d) profile obtained with ICP for the second increment position.</figcaption></figure>
# 7. Conclusions

In order to overcome the problem caused by the mechanical vibrations of a high-speed robot arm, we presented a new surface imaging method based on the classical point-to-point ICP registration algorithm. This method has been compared with standard surface imaging without using image alignment algorithm and has shown a significant improvement in surface image quality. The results of the surface reconstructions obtained on a synthetically-generated data also have shown that the resolving power of the ICP algorithm decreases for strong rotational vibrations, since it becomes quite difficult to efficiently handle partially overlapping images at large rotation angles. Employing robotic probe positioning, our method allowed us to understand the character of robot's vibrations and conclude that the mechanical movements of the robot used in inspection are affected by vertical vibrations. The results obtained are encouraging. However, although our method provides an accurate surface reconstruction, we have seen that it reduces computational power compared to standard stitching (without registration step) along the nominal probe trajectory.

Future work will include the application of the proposed concept to the case of more complex shaped surfaces and complex probe trajectories. A more advanced stitching method will also be investigated to improve the resolution of the ICP algorithm for strong rotational vibrations. Another investigation will consist of extending our stitching method for 3D surface imaging with a focus on optimizing its processing power.

## CRediT authorship contribution statement

**Ekaterina Iakovleva:** Conceptualization, Methodology, Software, Validation, Writing – original draft, Writing – review & editing. **David Roué:** Conceptualization, Validation, Writing – original draft, Writing – review & editing. **Philippe Brédif:** Conceptualization, Writing – original draft, Writing – review & editing, Methodology.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

The authors do not have permission to share data.

## References

[1] Hopkins D, Neau G, Le Ber L. Advanced phased-array technologies for ultrasonic inspection of complex composite parts. Proc., smart materials, structures, & NDT in aerospace 2011.

[2] Le Jeune L, Robert S, Dumas P, Membre A, Prada C. Adaptive ultrasonic imaging with the total focusing method for inspection of complex components immersed in water. In: AIP Conference proceedings, vol. 1650; March 2015. p. 1037-46.

[3] Malkin RE, Franklin AC, Bevan RL, Kikura H, Drinkwater BW. Surface reconstruction accuracy using ultrasonic arrays: application to non-destructive testing. NDT E Int 2018;96:26-34.

[4] McKee JG, Bevan RL, Wilcox PD, Malkin RE. Volumetric imaging through a doubly-curved surface using a 2D phased array. NDT E Int 2020;113:102260.

[5] Thomsen DK, Søe-Knudsen R, Balling O, Zhang X. Vibration control of industrial robot arms by multi-mode time-varying input shaping. Mech Mach Theor 2021; 155.

[6] Futami S, Kyura N, Hara S. Vibration Absorption control of industrial robots by acceleration feedback. In: IEEE Transactions on Industrial Electronics, Vols. %1 sur %2vol. IE-30, n° %13; 1983. p. 299-305.

[7] Lu F, Milios E. Robot pose estimation in unknown environments by matching 2d range scans. J Intell Rob Syst 1997;18:249-75.

[8] Minguez J, Lamiraux F, Montesano L. Metric-based scan matching algorithms for mobile robot displacement estimation. In: Proceedings of the 2005 IEEE International Conference on robotics and automation; 2005. p. 3557-63.

[9] Besl PJ, McKay ND. Method for registration of 3-D shapes. Sensor fusion IV 1992; 1611:586-606.

[10] Rusinkiewicz S, Levoy M. Efficient variants of the ICP algorithm. In: In Proceedings of the third International Conference on 3D Digital imaging and Modelling. Quebec City: Canada; May 2001. p. 145-52.

[11] Holmes C, Drinkwater BW, Wilcox PD. Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation. NDT&E Int Dec. 2005;38:701-11.

[12] Drinkwater BW, Wilcox PD. Ultrasonic arrays for nondestructive evaluation: a review. NDT&E Int. 2006;39:525-41.

[13] Johnson DH, Dudgeon DE. Array signal processing: Concepts and techniques. PTR Prentice Hall; 1993.

[14] Karaman M, O'Donnell M. Subaperture processing for ultrasonic imaging. IEEE Trans. UFFC. 1998;45:126-35.

[15] Martin CJ, Martinez O, Ullate LG, Octavio A, Godoy G. Reduction of grating lobes in SAFT images. In: 2008 IEEE ultrasonics symposium; November 2008. p. 721-4.

[16] Gu K, Wang G, Li J. Migration based sar imaging for ground penetrating radar systems. IEE Proc - Radar, Sonar Navig October 2004;151:317-25.

[17] Lévesque D, Blouin A, Néron C, Monchalin JP. Performance of laser-ultrasonic F-SAFT imaging. Ultrasonics 2002;40:1057-63.

[18] Cumming IG, Wong FH. Digital processing of synthetic aperture radar data, vol. 1. Norwood, Massachusetts: Artech house; 2005. p. 108-10.

[19] Getreuer P. Linear methods for image interpolation. Image Process Line 2011;1: 238-59.

[20] Shao-Yu P, Meng-Lin L. Linear array beamformation using virtual sub-wavelength receiving elements. In: IEEE International ultrasonics symposium IUS2015; 2015. 21-24 October.

[21] Arun KS, Huang TS, Blostein SD. Least square fitting of two 3-d point sets. In: IEEE Transactions on Pattern analysis and Machine Intelligence, vol. 9; 1987. p. 698-700.

[22] Horn BKP. Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society of America A April 1987;4:629-42.

[23] Walker MW, Shao L, Volz RA. Stimating 3-d location parameters using dual number quaternions. CVGIP Image Underst November 1991;54:358-67.

[24] Horn BKP, Hilden HM, Negahdaripour Sh. Closed-form solution of absolute orientation using orthonormal matrices. Journal of the Optical Society of America A July 1988;5:1127-35.

[25] Lorusso A, Eggert D, Fisher R. A comparison of four algorithms for estimating 3-D rigid transformations. In: Chez in Proceedings of the 4th British Machine vision Conference (BMVC '95). England: Birmingham; September 1995. p. 237-46.

[26] Zhang Z. Iterative point matching for registration of free-form curves and surfaces. Int J Comput Vis 1994;13:119-52.

[27] Bentley JL. K-d trees for semidynamic point sets. In: In Proceedings of the sixth annual symposium on Computational geometry; May 1990. p. 187-97.

[28] Friedman JH, Bentley JL, Finkel RA. An algorithm for finding best matches in logarithmic expected time. ACM Transaction on Mathematical Software September 1977;3:209-26.

[29] Simon D, Hebert M, Kanade T. Real-time 3-D pose estimation using a high-speed range sensor. In: Proceedings of IEEE International Conference on robotics and automation (ICRA '94). vol. 3; May 1994. p. 2235-41.

[30] Otsu N. A threshold selection method from gray-level histograms. IEEE transactions on systems, man, and cybernetics 1979;9:62-6.

[31] CIVA: Simulation software for Non-Destructive Testing,» [En ligne]. Available: https://www.extende.com/.

[32] Calmon P, Mahaut S, Chatillon S, Raillon R. CIVA: an expertise platform for simulation and processing NDT data. Ultrasonics 2006;44:e975-9.

[33] Gengembre N, Lhémery A. Pencil method in elastodynamics: application to ultrasonic field computation. Ultrasonics 2000;38:495-9.