# Review of 3D topography stitching and registration algorithms: procedure, error evaluation and mathematical modelling

Robin Guibert, Frederic Robache, Julie Lemesle, Saïd Hanafi, Christopher Brown, Maxence Bigerelle

► To cite this version:

Robin Guibert, Frederic Robache, Julie Lemesle, Saïd Hanafi, Christopher Brown, et al. Review of 3D topography stitching and registration algorithms: procedure, error evaluation and mathematical modelling. Surface Topography: Metrology and Properties, 2025, 13 (1), pp.013001. 10.1088/2051-672X/adabfd. hal-05231882

HAL Id: hal-05231882

https://uphf.hal.science/hal-05231882v1

Submitted on 1 Sep 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
# Review of 3D topography stitching and registration algorithms: procedure, error evaluation and mathematical modelling

Robin Guibert<fnref n="1" />,<fnref n="2" /> ID, Frederic Robache<fnref n="1" /> ID, Julie Lemesle<fnref n="3" />,<fnref n="4" /> ID, Saïd Hanafi<fnref n="1" />, Christopher A Brown<fnref n="5" /> ID and Maxence Bigerelle<fnref n="1" />,<fnref n="6" /> ID

E-mail: robinguibert.contact@gmail.com

Keywords: stitching, topography, roughness, registration

## Abstract

Surfaces are the privileged places of interaction between physical phenomena and objects. Roughness studies, especially when performing multiscale analysis, are tools of choice to understand physical phenomena and their scales of application. However, profilometers, especially optical systems, must compromise between field of measurement and resolutions. Stitching is an assembling technique aiming to solve this compromise by combining elementary maps, such as images or topographies. Stitching generates high resolution over a large field of measurement maps, which increases the measurable scale range and facilitates the correct identification of physical phenomena at their scales of application. This article proposes a review of 3D topography stitching algorithms. After explanations on the use cases of 3D topography stitching, the stitching procedure from elementary maps acquisition to the obtention of the stitched map is described step-by-step. Secondly, errors in measurement and stitching are presented with the sources of errors and the error evaluation methods. Lastly, the mathematical modelling of 3D topography is detailed to better understand the optimization process used in the in-plane and out-of-plane registration steps of the stitching algorithms. Comparison of algorithms involved in stitching are proposed so that researchers might find the most suitable algorithm to their needs. Overall, this work aims at introducing researchers and metrologists to important multidisciplinary notions for the use and design of 3D topography stitching algorithms and offers a tutorial-based approach.

## Abbreviations

<table><tr><td>SPM</td><td>Scanning Probe Microscopy</td></tr><tr><td>AFM</td><td>Atomic Force Microscope</td></tr><tr><td>FFT</td><td>Fast Fourier Transform</td></tr><tr><td>RMS</td><td>Root Mean Square</td></tr><tr><td>CNC</td><td>Computerized Numerical Control</td></tr><tr><td>SEM</td><td>Scanning Electron Microscope</td></tr><tr><td>FIB</td><td>Focused Ion Beam</td></tr><tr><td>AO</td><td>Alternative Optimization</td></tr></table>

<table><tr><td>EBL</td><td>Electron Beam Lithography</td></tr><tr><td>CT</td><td>Computed Tomography</td></tr><tr><td>SIFT</td><td>Scale Invariant Feature Transform</td></tr><tr><td>SURF</td><td>Speeded-Up Robust Features</td></tr><tr><td>BRISK</td><td>Binary Robust Invariant Scalable Keypoints</td></tr><tr><td>MSER</td><td>Maximally Stable Extremal Regions</td></tr><tr><td>ORB</td><td>Oriented FAST and Rotated Brief</td></tr></table>

<footnote n="1">¹ Univ. Polytechnique Hauts-de-France, CNRS, UMR 8201 - LAMIH - Laboratoire d'Automatique de Mécanique et d'Informatique Industrielles et Humaines, F-59313 Valenciennes, France</footnote>

<footnote n="2">² Anscale, 76600 Le Havre, France</footnote>

<footnote n="3">³ Valutec, Univ. Polytechnique Hauts-de-France, CEDEX 9, 59314 Valenciennes, France</footnote>

<footnote n="4">⁴ U.R Concept, 59300 Valenciennes, France</footnote>

<footnote n="5">⁵ Surface Metrology Laboratory, Worcester Polytechnic Institute, Worcester, MA 01609, United States of America</footnote>

<footnote n="6">⁶ INSA Hauts-de-France, F-59313 Valenciennes, France</footnote>
<table><tr><td>RANSAC</td><td>RANdom SAmple<br>Consensus</td></tr><tr><td>MSAC</td><td>M-estimator Sample<br>Consensus</td></tr><tr><td>PCA</td><td>Principal Component<br>Analysis</td></tr><tr><td>SiC</td><td>Silicon Carbide</td></tr><tr><td>SEE</td><td>Stitching Error Estimator</td></tr><tr><td>MLE</td><td>Maximum Likelihood<br>Estimation</td></tr><tr><td>DOF</td><td>Degree of Freedom</td></tr></table>

# 1. Introduction

Physical phenomena and their couplings set the world in motion. The expression of such physical phenomena occurs on the interface between objects and their environment. In other words, surfaces are the privileged place where physical phenomena are expressed. Roughness studies use surface topography to provide a better understanding of such mechanisms. Also, scale, i.e., a collection of wavelengths in surface morphologies, is an important factor, as each physical phenomenon occurs at a given scale. Therefore, the use of multiscale decomposition to separate physical phenomena by their scale of applications facilitates the study of such physical phenomena. It is then necessary to acquire as much scales as possible to increase the probability of detection of the morphological expression of physical phenomena. The limitation of scale range is in fact inherent to surface metrology: a compromise exists between field-of-view and resolution. A higher resolution is accompanied by a smaller field-of-view. Using an unsuitable resolution would impede the detection of morphological structures linked to the studied physical phenomena. On the other hand, unsuitable field-of-view would interfere with statistical analysis if not enough morphological structures were measured, especially in the case of large structures, resulting in a kind of cherry-picking bias. Setting the correct balance between resolution and field-of-view is challenging, as both expected scales of application of the physical phenomena, technical limitations of measurement techniques and available parts for the measurement systems should be considered. In fact, the selection of a measurement technique is directly linked to the expected scales to observe. Atomic Force Microscope (AFM) produces high-resolution maps at the cost of a limited field-of-measurement, generally a few microns. For a larger map size, interferometer can be used for sub-micrometer resolution over hundreds of microns maps, followed by confocal microscope and focus variation microscope. Tactile profilometer can be useful for larger map, but resolution is then limited by the acquisition time of the topography. In optical system, resolution and field-of-view depend on the objective lens: high magnification is accompanied with smaller field-of-

measurement. For Scanning Probe Microscopes (SPM), probes with various tip sizes influence lateral resolution: a smaller tip increase resolution. However, field-of-measurement is not fixed for SPM and can be modified by metrologist. The drawback is then the acquisition time, as each point in captured one at a time, contrary to optical systems that measure several points at a time.

Stitching is an assembly technique of elementary maps, such as topographies, images or any type of maps, which produces a seamless high resolution over a large field of measurement stitched map. If enough similarities are present in the overlapping regions of elementary maps, stitching can be of interest, including for visual roughness measurements that use images directly to estimate a specimen roughness. Also, lateral and vertical resolutions remain the same as elementary maps. The stitched map then contains a larger scale range, which is essential to multiscale roughness studies. In topography, stitching is widely used, especially with 3D optical profilometers. Figure 1 compares an elementary map measured over 140 × 105 µm (640 × 480) and a stitch map measured over 2001 × 2001 µm (9127 × 9127 points) from an abraded polymer disk measured with a white light interferometer Zygo NewView™ 7300 (Zygo, Middlefield, Connecticut, USA), and a 50x objective lens. Lateral resolutions on both topographies are 0.219 µm but the field of measurement is more than 14 times bigger on the X-axis and more than 18 times on the Y-axis on the stitched map. The global view of the stitched topography reveals a groove network, where each groove is intersected by other grooves, i.e., a fractal structure.

However, stitching algorithms are often complex, potentially computation-intensive, and misunderstood. Stitching by essence is a multidisciplinary field, combining mechanics, computer science and mathematics, which creates a steep learning curve for any researcher new to stitching. This article proposes a review of stitching algorithms, with a focus on 3D topography stitching. A second objective of this review is to provide a guided tutorial to better understand each step involved in 3D topography stitching. The use cases of 3D topography stitching, the stitching procedure, its error evaluation and the mathematical modelling of the stitching problem are presented with references to the literature. Comparison of algorithms and techniques, along with sum-up tables, provides elements for the selection of a suitable stitching algorithm and an overview of the notions useful in stitching.

# 2. Use cases of 3D topography stitching

Stitching of 3D topographies can be used for a large variety of cases, making it a versatile tool for surface metrology and scale analysis. First case is the measurement of high-resolution topographies over a large field
<figure><img src="image_2.png" /><figcaption>**Figure 1.** An elementary topography and a stitched topography. Both maps have the same lateral resolutions, but the stitched map contains more points over a larger field of measurement.</figcaption></figure>

of measurement. By combining elementary topographic maps, stitching increases the field of measurement. It is then possible to use a 50x objective lens for high resolution, and to stitch many elementary topographies to obtain a field of measurement equivalent to a 5x objective lens, or much larger. Stitching does not increase resolution directly but allows the use of objective lens with better resolution, which field-of-measurement is increased by stitching. The gain in resolution is then quantified as the ratio between the resolution of the higher magnification objective lens over the resolution of the lower magnification objective lens. Stitching allows the capture of the larger points-of-interest of the surface morphology, while keeping the smaller details observable (Yu and Peng [1]). It is possible to use stitching to capture impossible field of measurement, much larger than what objective lens can offer. Large aspheric mirrors with diameter over one meter were measured with stitching by Yan et al [2] as precision over a large surface was required for the correct measurement of new large telescopes (several tens of meters) and also to extend the vertical dynamic range that incorporates the overall form of the mirror. Using stitching also increases the scale range, which is essential to measure as much as possible scales for multiscale analysis (Brown et al [3]). Scale can be defined as an ensemble of wavelengths (or frequencies) as explained in Brown et al [3]. By capturing more scales with stitching, less assumptions about the relevant scales are needed, which dictate the selection of measurement apparatus, objective lens for optical system and probe for Scanning Probe Microscopy (SPM). The use of multiscale analysis can then identify the relevant physical phenomena at the relevant scales. The number of points in a stitched map is also higher, which brings stability to 3D roughness parameters (ISO 25 178-2 [4], ASME B46.1 [5]) as most are based on statistics.

The second stitching case is shown in figure 2: stitching with a change of objective lens to get better lateral resolutions on the same field of measurement. The specimen is a sandblasted TA6V rod. Topographies were measured with a Bruker ContourGT (Bruker, Billerica, Massachusetts, USA) equipped with 5×, 20× and 50× objective lenses, five times in each configuration. The field of measurement is 1258 × 943 µm. Lateral resolutions are 1.97 µm for the 5× (640 × 480 points), 0.494 µm for the 20× (2552 × 1915 points) and 0.198 µm for the 50× (637 × 4782 points). Details on the stitched maps are more observable and the higher number of points allows a better statistical stability when computing roughness parameters.

Roughness parameters can then be computed on these maps. The coefficient of variation measures how roughness parameter values vary relatively to its mean value (table 1). Using stitching for 3D topographies can increase the statistical stability and relevance of roughness parameters by allowing the use of objective lens with higher magnification. Also, note that some roughness parameters are highly influenced by the sampling interval, such as Sds, Sdq and Sdr, and so, their coefficient of variation increases more compared to other parameters when using stitching without resampling.

It is also easier to measure a statistically representative topography for a given morphology, as maps can be larger with 3D topography stitching. Moreover, this facilitates the study of heterogeneous surfaces, including their interface, with a large enough field of measurement to capture all topographic sub-regions. Finally, more and more researchers use stitching with topography measurement system to capture the form of surfaces. Wang et al [6] evaluate the material removal distribution on free-form polished surfaces. This application requires large field of measurement with the constraint to follow a free-form surface and
<figure><img src="image_3.png" /><figcaption>Objective lens 20x – 5 x 5 stitching – Lateral resolutions: 0.494 µm</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Objective lens 50x – 10 x 13 stitching – Lateral resolutions: 0.198 µm</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Sandblasted topographies, with and without stitching. Field of measurement is the same on all topographies, but different objective lenses were used (5×, 20× and 50×) and lateral resolutions are higher.</figcaption></figure>

stitching allows the extension of the vertical scanning range by merging elementary maps with shifted vertical scan range. Stitching was used for the measurements of spheres by Laubach *et al* [7] with a line-scanning interferometer that follows the specimen curvature thanks to a five-axis system while the

specimen is fixed on a rotation table. The stitching algorithm corrected displacements generated by the rotation of the five-axis system around the specimen curvature. Similarly, the deformation of spherical surfaces with high numerical aperture was measured by Pan *et al* [8] with stitching where the specimen was
<table><caption>**Table 1.** Coefficients of variation from various roughness and morphological parameters influenced by stitching and the use of different objective lenses. Some parameters vary 100 times more than others.</caption><thead><tr><th>Ranking</th><th>Parameter</th><th>Coefficient of variation [%]</th><th>Ranking</th><th>Parameter</th><th>Coefficient of variation [%]</th></tr></thead><tbody><tr><td>1</td><td>Z center</td><td>216.6949</td><td>36</td><td>Sku</td><td>13.3222</td></tr><tr><td>2</td><td>Sdv</td><td>142.1763</td><td>37</td><td>Sbi</td><td>12.9695</td></tr><tr><td>3</td><td>Ssk</td><td>134.3876</td><td>38</td><td>Z value</td><td>12.4920</td></tr><tr><td>4</td><td>Sds</td><td>133.442</td><td>39</td><td>Svk</td><td>12.3637</td></tr><tr><td>5</td><td>Sda</td><td>123.638</td><td>40</td><td>Sq</td><td>10.8734</td></tr><tr><td>6</td><td>Spc</td><td>115.905</td><td>41</td><td>Sv</td><td>10.8057</td></tr><tr><td>7</td><td>Ssc</td><td>105.0488</td><td>42</td><td>S10z</td><td>10.5282</td></tr><tr><td>8</td><td>Motif sphere radius</td><td>102.3888</td><td>43</td><td>Mean depth of furrows</td><td>9.9654</td></tr><tr><td>9</td><td>Density of furrows</td><td>95.7312</td><td>44</td><td>Vvc</td><td>9.9128</td></tr><tr><td>10</td><td>Smr</td><td>88.1655</td><td>45</td><td>Vv</td><td>9.8837</td></tr><tr><td>11</td><td>Sdr</td><td>82.8385</td><td>46</td><td>Vvv</td><td>9.6600</td></tr><tr><td>12</td><td>Volume</td><td>79.7255</td><td>47</td><td>Sal</td><td>9.4495</td></tr><tr><td>13</td><td>Sha</td><td>77.2525</td><td>48</td><td>Sa</td><td>9.4274</td></tr><tr><td>14</td><td>Shv</td><td>76.6252</td><td>49</td><td>Smc</td><td>9.3992</td></tr><tr><td>15</td><td>Area</td><td>72.4082</td><td>50</td><td>Motif height</td><td>9.1112</td></tr><tr><td>16</td><td>Number of motifs</td><td>65.0246</td><td>51</td><td>Sxp</td><td>8.3162</td></tr><tr><td>17</td><td>Coflatness</td><td>63.0420</td><td>52</td><td>Sp</td><td>8.0957</td></tr><tr><td>18</td><td>Spd</td><td>59.9483</td><td>53</td><td>Vmc</td><td>7.7084</td></tr><tr><td>19</td><td>Sdq</td><td>56.0107</td><td>54</td><td>Sr1</td><td>7.1591</td></tr><tr><td>20</td><td>Motif min pitch</td><td>50.8688</td><td>55</td><td>Sk</td><td>7.0264</td></tr><tr><td>21</td><td>Minimum diameter</td><td>46.1627</td><td>56</td><td>Sz</td><td>6.5810</td></tr><tr><td>22</td><td>Mean diameter</td><td>44.3771</td><td>57</td><td>Aspect ratio</td><td>5.1919</td></tr><tr><td>23</td><td>Equivalent diameter</td><td>44.1157</td><td>58</td><td>Y value</td><td>3.4495</td></tr><tr><td>24</td><td>Maximum diameter</td><td>42.8822</td><td>59</td><td>Y center</td><td>3.2994</td></tr><tr><td>25</td><td>Motif pitch</td><td>39.3796</td><td>60</td><td>Motif roundness</td><td>2.5821</td></tr><tr><td>26</td><td>Motif max pitch</td><td>33.0263</td><td>61</td><td>Str</td><td>2.3782</td></tr><tr><td>27</td><td>Perimeter</td><td>32.7246</td><td>62</td><td>Motif orientation</td><td>2.2056</td></tr><tr><td>28</td><td>Maximum depth of furrows</td><td>20.2932</td><td>63</td><td>Svi</td><td>1.9955</td></tr><tr><td>29</td><td>Vm</td><td>20.2932</td><td>64</td><td>X value</td><td>1.7205</td></tr><tr><td>30</td><td>Vmp</td><td>20.2932</td><td>65</td><td>X center</td><td>1.6123</td></tr><tr><td>31</td><td>Spk</td><td>19.3794</td><td>66</td><td>Compactness</td><td>1.3706</td></tr><tr><td>32</td><td>S5p</td><td>17.1050</td><td>67</td><td>Sci</td><td>1.1667</td></tr><tr><td>33</td><td>S5v</td><td>16.0201</td><td>68</td><td>Number of neighboring motifs</td><td>1.0394</td></tr><tr><td>34</td><td>Std</td><td>15.1617</td><td>69</td><td>Sfd</td><td>0.8979</td></tr><tr><td>35</td><td>Motif form factor</td><td>14.7095</td><td>70</td><td>Sr2</td><td>0.1734</td></tr></tbody></table>

rotated between each position. A comparison of stitching algorithms for the measurements of cylindrical surfaces with an interferometer was proposed by Chen *et al* [9], highlighting the need for wider dynamic range and larger field of measurement in interferometry. Kong and Chen [10] used stitching for large aperture acylindrical surfaces with mild aspherical departure in order to limit the complexity of large optics testing as well as its cost. Chen *et al* [11] indicated that large steep convex spheres with small ratio of the radius of curvature to the aperture (R/number) could not be captured in a single measurement and stitching overcame this limitation. Chen *et al* [12] also used stitching to increase lateral resolution and the vertical range of measurement in interferometry for large optical surfaces.

The relevance of stitching for some SPM measurement systems is lower than for optical systems, as SPM can easily modify their field of measurement. Some

configurations of SPM, such as large field-of-view AFM, uses high-end displacement table. However, such system is expensive compared to using stitching with a classic system, as limited range displacement table can focus on robustness and precision and can extend their field of measurement by stitching algorithms for a lower price. Also, vertical scanning range is limited on large field of view AFM, requiring very flat surfaces and increases the risk of undesired contact during displacement that could damage the measurement system or the surface. Stitching then outcomes this large field-of-view AFM configuration.

## 3. The stitching procedure

### 3.1. Step-by-step review of stitching algorithms

Stitching of elementary topographic maps is complex, and many steps are involved in the stitching
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Stitching procedure of 3D topographies. This stitching procedure highlights the importance of acquisition preparation as the quality of elementary maps needs to be high so that the stitching algorithm determines the correct registration before map merging. Reproduced from [14]. CC BY 4.0.</figcaption></figure>

procedure, combining advanced algorithms. The overall goal of stitching algorithms is to compute the geometric transformations, such as translations or rotations, to apply on each elementary map, so that overlapping regions are aligned between neighboring maps, i.e., all elementary maps are registered. Transformed elementary maps are then merged with algorithms preventing visible discontinuities and fusion artefacts in the stitched map. The stitching procedure is summed up by Marinello et al [13] and by Lemesle et al [14]. Figure 3 presents the stitching procedure. Even though stitching algorithms might differ, especially in the registration steps, this stitching procedure is always followed.

## 3.2. Data acquisition

The data acquisition step is the most important one for stitching, as the quality of topographic elementary maps and the accurate and robust positioning between the measurement of two elementary maps are linked to the final quality of the stitched map and to the difficulty to align each elementary map. Stitching might even be impossible if the quality of elementary maps is poor. Missing values, predominant noise, illumination heterogeneity are common quality problems in topographic measurement system. The positioning of each elementary map can either be free, which is called unstructured stitching (Huang et al [15]), either predetermined using a grid, rectangular (Pan et al [8], Wang et al [6]) or spherical (O'Donohue et al [16]) depending on the measurement system, which is called structured stitching. In either case, a common part in the elementary maps must be present to be correctly aligned during stitching. This common part is called an overlap and the commonly used overlap percentage is comprised between 10% and 30% in both X and Y axes. Chen et al [17] suggest at least 10% overlap between elementary maps. On the other hand, Seo et al [18] use 30% overlap for its stitching of a freeform surface. Several works use 20% overlap (Guibert et al [19], Wang et al [6], Lu et al [20], Deltombe et al [21] and Deltombe et al [22]). No influence of overlap over stitching precision in 3D topography stitching has been found by Marinello et al [13] for 50%, 60%, 70%, 80% and 90%, which are very

high overlap values. A larger overlap percentage offers more points for the correct alignment of elementary maps, but also increases the number of elementary maps to measure, and so, increases the total acquisition time and the stitching time. The choice of the overlap percentage is a compromise between a large enough overlapping region for the map comparison in the registration steps and a small enough overlapping region to limit the number of elementary map acquisitions.

Some paper described the use of multiaxis system for the stitching of freeform surfaces (Seo et al [18], Spichtinger et al [23] or Lu et al [20] in the case of measurement of spherical Fresnel lens). The interest is to keep the normal alignment between the surface and the objective lens for each elementary map. In this case, the quality of each elementary map is higher, and stitching is easier as there is less noise on elementary maps. However, more degrees-of-freedom should be considered in the stitching algorithm compared to a 3-axis translation system. In the case of surface topography, common profilometers are interferometers, focus variation microscopes, confocal microscopes, tactile profilometers and Atomic Force Microscopes (AFM) according to Leach [24]. The selection of the measurement system depends on the availability of the system, its speed, its precision, the specimen surface integrity and morphology, and finally the specimen size.

Stitching was used with interferometer in optics by Murphy et al [25] and Cimma et al [26], in polymer abrasion by Guibert et al [19], in artistic painting by Mironova et al [27] and with a focus on fractal dimension by Bigerelle et al [28]. Focus variation microscope was also used in the study of Mironova et al [27]. Stitching of 2D profiles obtained by laser confocal microscopy and tactile profilometer was performed by Fu et al [29]. Other researchers focused on using stitching with AFM, such as Marinello et al [13], Montgomery et al [30] and Zhao and Yang [31]. Note that while stitching can be used either in 2D or 3D measurement system, it is more relevant for field optical systems, which measure a fixed field of measurement, than for Scanning Probe Microscopes (SPM) that have a variable field of measurement.
### 3.3. Preprocessing

The second step of the stitching procedure is the preprocessing step. Once elementary maps are acquired, it can be useful to apply preprocessing to facilitate stitching. Preprocessing can be data cleaning, by using Fast Fourier Transform (FFT)-based filtering (Pan et al [8]), error removal using gaussian process (Liu et al [32]), thresholding with the help of Otsu's thresholding (Otsu [33]) or height signal and noise separation (Lei et al [34]). The determination of reference plane, by removing form using polynomials (generally, a 1st order polynomial for stitching, but 2nd order can be used for symmetrical surfaces and 3rd order for asymmetrical surfaces) can help the stitching algorithm by limiting the vertical drift between elementary maps. Also, some transformations might be applied to facilitate the stitching procedure, such as using directional filters (Sobel X and Sobel Y, as explained by Kanopoulos et al [35], Gao et al [36] and Zhang et al [37]) to highlight the contours of morphological structures. For non-directional highlighting of morphological structures, gradient-based operators are useful for stitching (gradient, Laplacian, or composite filters such as Laplacian of gaussian filters or differences of gaussian filters). Note that noise can be amplified if the elementary map quality is poor or if multiple operators are applied.

### 3.4. In-plane registration

Then, elementary maps need to be aligned using the overlapping sub-regions of the maps. This alignment process is called the registration step. In the case of 3D topography, registration can be separated in two sub-steps: in-plane registration for XY-translations and Z-rotation, and out-of-plane registration for Z-translation and XY-rotations. Both in-plane and out-of-plane registrations are required for 3D topography stitching and are performed either successively or at the same time. Outside of stitching, registration can also be useful for the comparison of surfaces before and after a process to evaluate the morphological transfer function of the process [38] or for colocalization.

Five categories exist for in-plane registration: manual stitching, geometric algorithm, fiducialization, global optimization and computer vision-based algorithms. Figure 4 presents a classification of in-plane registration algorithms. First category is manual stitching, which consists in the manual identification of displacements to apply to register two elementary maps. Metrologists can then use transparency while overlapping two elementary maps, until a visual matching is found. This technique depends on the metrologist's visual ability and experience and no quantification is possible, making this technique hard to evaluate. Alternatively, the metrologist can rely on colocalization: moving the elementary maps is still manual but a numerical value can be computed, giving

more indications to help registration. Typically, cross-correlation and root mean square (RMS) difference can be used to respectively find maximum similarities or minimum differences between the overlapping elementary maps. This technique is more robust than the transparency technique, and also provides quantification. Another technique is semi-automatic: the region pointing technique. A metrologist indicates small sub-regions on each elementary map and applies an automatic stitching algorithm between these sub-regions. By doing so, the solution space for the registration is reduced and the algorithm run time is lower. The region pointing technique is more precise than the two previous techniques, as finer displacements to apply on elementary maps during stitching are detected with an automatic stitching algorithm. Manual stitching technique are simple but suffer from three drawbacks. First, the metrologist's role is significant, which adds arbitrary choice. Thus, repeatability is poor with manual stitching techniques. Also, nothing prevents false inputs in the process, as a metrologist can make mistakes in the determination of displacements to apply or in the selection of relevant sub-regions, especially in medium and large size stitching. Manual registration is prone to errors and lacks the robustness of automatic stitching algorithm. Lastly, the speed of the stitching procedure using manual registration is lacking, as a metrologist must interact at each local stitching and cannot be as fast as a computer for such repetitive task. This also limits the maximum number of elementary maps, as stitching 10 000 elementary maps would take too much time. Montgomery et al [30] indicated that a manual stitching could reach 30 min per image pair, in contrast to automatic stitching which only takes a few seconds.

Geometric stitching is the second kind of stitching algorithm. This is the most basic automatic stitching algorithm, and most manufacturers use it in their profilmometers, as it always gives a stitched map, although not the best registered. This technique is based on the metadata of the elementary maps. Geometric stitching is fast and robust but depends on the accuracy of the measurement system, especially of the displacement table. It is even more useful in structured stitching. Absolute positions in XYZ-axes can be used, as done by Dong et al [39] who used the absolute position of each elementary topography of concave spherical mirrors in a circular grid. The relative position of each elementary map in the stitching grid is also used in geometric stitching. Also, hardware-assisted stitching, using extra information such as accelerometers, which can measure the position and rotation of each map during the acquisition. King and Bibby [40] used an interferometer, combined with a 5-axes Computerized Numerical Control (CNC) table and a 7-axis partholder allowing rotary tilting for the topographic measurement of optics parts. However, geometric stitching lacks in accuracy, despite its high speed. It is nonetheless interesting as it could be used as
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Categories of in-plane registration algorithms, summing-up techniques used in manual and automatic stitching algorithms.</figcaption></figure>

prepositioning, before a finer registration with a more accurate stitching technique, as explained by Paul *et al* [41].

The third category is fiducialization, which consists in adding markers, i.e., points-of-interest, over the surface. Stitching then focuses on these markers to perform in-plane registration. Adding markers is interesting as distinctive known shapes are easily detected for in-plane registration, as the markers can be large enough to make the signal noise negligible in comparison of such large structures. King [42] used fiducialization with an alignment of diamond markers. Kammers and Daly [43] performed Scanning Electron Microscope (SEM) image stitching by adding spherical markers in the case of small scale measurements. Markers can use unique shapes, such as circle, diamond or square, or all being similar with a special pattern allowing the unicity of the solution of in-plane registration. Focused Ion Beam (FIB) can be used for the addition of markers on the sample surface (Rust and Todd [44]), as well as powders or droplets (needling drying before measurement). Also, native markers such as in microfluidic devices can already exist on the sample as shown in Xu *et al* [45]. The main drawback of fiducialization is the non-conservation of the original object surface and surface integrity cannot be preserved as markers are added. Postprocessing is then required to avoid perturbation during the calculation of roughness parameters. The surface marking can also be difficult for small scale surfaces, and a not

enough accurate marking process can pollute the smaller scale.

In-plane registration can also be based on global optimization techniques to determine the correct in-plane displacements to apply on elementary maps. Global optimization techniques in stitching usually aim to minimize differences between overlapping regions of the elementary maps in a least square sense. A stitching algorithm based on global optimization was used by Chen *et al* [9] to compute the correct displacements, Δx, Δy and Δz for the registration steps. Wang *et al* [6] modeled the error minimization problem as a linear least-square problem and use QR factorization to solve it. Chen *et al* [17] also modelled their stitching with two subproblems: first, a configuration optimization problem, i.e., a linear least-square problem solved by using blockwise QR decomposition combined with the sparse technique, and secondly, a non-linear optimization problem solved with iterative Alternative Optimization (AO) to minimize displacements errors in a 6-axis measurement system. This multiaxis measurement system offered the advantage to measure each elementary map normally, despite the curvature of the measured spherical lens. Note that the stitching algorithm optimizes all elementary maps displacements at once and globally, in contrast to sequential stitching. These techniques are more accurate than geometric in-plane registration, but also more noise-sensitive, which makes the stitching of flat or small-scale surfaces difficult. To speed-up the global optimization process and to increase
accuracy, pre-positioning can be used before a finer
registration based on global optimization, which reduc-
ces the number of solutions to explore during the
optimization process. Note that global optimizations
are computation-intensive, and the algorithm run
time increases with higher degrees of freedom and the
number of points on each elementary map. Selecting a
suitable global optimization algorithm and a suitable
objective function also contributes to the run time of
the stitching procedure and the registration precision.

Lastly, computer vision can be used for in-plane
registration. Computer vision is a subfield of comp-
uter science focused on the detection of points-of-
interest and objects in a given dataset, such as an
image. Two categories of computer vision-based algo-
rithms are used in stitching: direct methods and fea-
ture detection methods (Pravenaa [46]). Direct
methods explore the whole elementary map to find
similarities between overlapping regions, such as
cross-correlation-based algorithms. Correct displace-
ments are determined when cross-correlation is max-
imized (Marinello et al [13]). As cross-correlation is
based on convolution, such method is computation-
intensive and might not be suitable for large stitching
involving a high number of elementary maps due to a
high run time. Cross-correlation-based stitching was
used by Paul et al [41] on images of an Electron Beam
Lithography (EBL) and by Lemesle et al [14] in a com-
parison with proposed novel 3D topography stitching
algorithms. This category of computer vision-based
stitching algorithm is still noise-sensitive, but less than
global optimization techniques.

On the other hand, feature detection methods are
two-step processes, where points-of-interest, called
features, are first detected, before applying in-plane
registration only on the detected features. This strategy
is especially interesting for large elementary maps,
where direct methods can be much slower. Feature
detection algorithms were first developed for images
but can also be applied on topographic maps (Lu et al
[20]) and point clouds. A review of 3D point cloud
stitching was proposed by Agier [47], alongside with a
mathematical model of the optimization process in
stitching, and stitching was applied to 3D Computed
Tomography (CT) scans. Numerous algorithms exist
for feature detection and are called descriptors. First
major descriptor is the Scale Invariant Feature Trans-
form (SIFT) developed by Lowe [48]. Improvement
and new use cases of SIFT were also proposed in the
following years by Brown and Lowe [49], and Brown
and Lowe [50], before using SIFT for automatic image
stitching (Brown and Lowe [51]). SIFT is now widely
used in image recognition (Moretti et al [38]), stitching
and machine learning. An open-source detector called
Speeded-Up Robust Features (SURF) was proposed by
Bay et al [52] and was adapted for 3D topography
stitching by Li et al [53]. Derived descriptors were
developed by Valgren and Lilienthal [54], such as
U-SURF, which does not compute out-of-plane

rotations and is faster than SURF, or SURF-128, which
is faster and more accurate than SURF. Feature
descriptors can either detect features focusing on cor-
ners, such as the Harris Corner descriptor (Harris and
Stephens [55]), or on lines, such as Binary Robust
Invariant Scalable Keypoints (BRISK, Leutenegger et al
[56]), FAST (Rosten and Drummond [57]), Maximally
Stable Extremal Regions (MSER, Matas et al [58] and
Obdržálek et al [59]), Min Eigen (Shi and Tomasi
[60]), Oriented FAST and Rotated BRIEF (ORB,
Rublee et al [61]) and KAZE (Alcantarilla et al [62]).
The main differences between these feature descrip-
tors are the type of detected features (corner or edge),
the description of feature orientation, the number of
neighboring pixels considered during feature extrac-
tion and the speed of the feature detection algorithm.
These features descriptors were applied on over-
lapping regions of two topographic maps in figure 5.
KAZE descriptor was the most suitable on these topo-
graphic maps (this can change depending on the sur-
face morphology and scale), as more features with high
stability are found. On the contrary, FAST and MSER
descriptors poorly performed. These detection pro-
blems can be due to signal noise, false positive errors
or similarities (such as periodicity). To avoid influence
from outlier features, feature detection stitching algo-
rithms use algorithms to select stronger features, such
as RANdom Sample Consensus (RANSAC) algorithm
(Moretti et al [38] and Juan and Oubong [63]) or a
MATLAB equivalent algorithm called M-estimator
SAmple Consensus (MSAC). While features selected
by these algorithms have a higher stability, a random
component is used, and the repeatability of the stitch-
ing process is drastically reduced. Multiple different
stitching outputs can be obtained when using feature
detection algorithms.

To compare in-plane registration algorithms, a
3 × 3 stitching grid was measured on a mirror-
polished copper rod with indents from a macro-
indenter by a focus variation microscope Alicona Infi-
niteFocus G5 (Alicona, Raaba, Austria) and a 10x
objective lens. Each elementary topography is linked
to a focus-stacked image, as, in focus variation micro-
scopy, topography is computed from a stack of images
where height is determined when a point is focused.
Stitching was performed with eight different in-plane
registration algorithms. Figure 6 shows the center part
of the stitched topography, including overlapping
regions. Geometric stitching, global optimization
stitching with genetic algorithm and simulated
annealing, direct methods with cross-correlation on
image and the rising sun algorithm variants (naive
reflectance stitching and two-step reflectance stitch-
ing), and feature detection methods with SURF and
KAZE features produce good quality stitching. How-
ever, the left indent in deformed for geometric stitch-
ing and global optimization stitching, as it is located in
overlapping regions of elementary maps, and in-plane
registrations presented small alignment errors. This
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of feature descriptors for two overlapping regions of elementary topographic maps. Red circles and green crosses represent the locations of the features detected by each feature descriptor respectively in the left and right topographies to register. Similar features are linked together with yellow lines. The number of detected features and feature matching vary with the feature descriptor.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison of the central part of a 3 × 3 stitched map. Geometric stitching (a), direct methods (b, e, f), feature detection methods (c, d) and global optimization stitching (g, h) succeeded in stitching, but the left indent is slightly deformed in (a), (g) and (h) due to small registration errors.</figcaption></figure>
<table><caption>**Table 2.** Comparison of in-plane registration methods. Each stitching algorithm category compromises between registration accuracy and speed. This table allows metrologists to select a suitable algorithm depending on their project constraints.</caption><thead><tr><th>Method</th><th>Accuracy</th><th>Robustness</th><th>Speed</th><th>Surface integrity</th></tr></thead><tbody><tr><td>Manual</td><td>Operator dependent</td><td>—</td><td>—</td><td>+++</td></tr><tr><td>Geometric</td><td>Apparatus dependent</td><td>+++</td><td>+++</td><td>+++</td></tr><tr><td>Fiducialization</td><td>+++</td><td>++</td><td>++</td><td>—</td></tr><tr><td>Global optimization</td><td>++</td><td>++</td><td>—</td><td>+++</td></tr><tr><td>Direct</td><td>+++</td><td>++</td><td>—</td><td>+++</td></tr><tr><td>Feature detection</td><td>+++</td><td>—</td><td>++</td><td>+++</td></tr></tbody></table>

deformation can be easily detected as indent has a known standard shape. However, such defect can remain undetected in surface with complex and unknown morphology.

A comparative table of the five in-plane registration algorithm categories is provided in table 2. Manual stitching can be a first approach method but its poor robustness, speed and being prone to human errors make it unsuitable for large stitching. Automatic methods are then more interesting as human error is removed, although measurement and stitching errors still exist. Fiducialization can be used as a last resort method for nearly pathological cases of 3D topography stitching, by adding markers to the surface that can be easily detected. However, as surface integrity is not conserved, other automatic stitching in-plane registration methods should be tested beforehand.

The selection of a suitable in-plane registration algorithm also depends on surface morphology. For elementary maps with distinctive features, cross-correlation-based and feature detection algorithms are of interest, especially at medium and larger scales. For medium scale and homogeneous surfaces, global optimization of differences can be used. The key for correct registration is the ratio between surface features and uncertainty amplitude. As this ratio is smaller, such as when features are at small scales, more complex in-plane registration algorithms should be used, such as multimap-based algorithm. Lemesle et al [14] provided a benchmark of some in-plane registration algorithms on four test cases with various surface type, whose data are gathered in table 3. This benchmark is interesting, as except for the tested in-plane registration part, all other parts of the algorithms are the same. Therefore, differences in performance are only due to differences in the in-plane registration algorithms. Performance of each tested algorithm is evaluated on run time, memory and accuracy using mean repositioning error and Stitching Error Estimator (SEE). A clear trade-off between run time and accuracy is shown: more accurate registration algorithms require more time. Multimap-based in-plane registration algorithms (naïve reflectance, exhaustive enumeration reflectance and two-step reflectance, variants of the rising sun stitching algorithm) outperform other techniques in terms of success and accuracy, but are the slowest. Note that exhaustive

enumeration reflectance and two-step reflectance algorithms have the same accuracy, but the two-step reflectance is optimized regarding speed based on heuristics for the determination of a correct reflectance angle. Nonetheless geometric stitching is a good fallback technique, especially for pathological cases, such as perfectly plane surface or periodic surface with no noise, provided that the measurement system accuracy and robustness are good enough. If the measurement system is not accurate and robust, geometric stitching should not be used, and fiducialization should be preferred.

### 3.5. Out-of-plane registration

As topographies are 3D maps, stitching algorithms must also handle out-of-plane registration, i.e., correction of the translation along Z-axis and the rotation along X and Y axes to compensate rigid body errors (figure 7). Out-of-plane registration is necessary to avoid vibration effect and vertical drifts.

A first approach method for out-of-plane registration computes the difference of mean planes between overlapping regions of elementary maps. Registration is then performed by removing this difference to one of the elementary maps. This method is fast and efficient for low noise surfaces or flat surfaces, but it is not suitable for more complex surfaces, such as curved or plateau surface. Vertical drift might not be fully compensated, and no out-of-plane rotation is compensated. Other techniques are the use of Principal Component Analysis (PCA) or first order polynomial to remove the first order form of the elementary maps. While out-of-plane rotations are applied with both techniques, complex morphology is still a challenge to register. A more precise out-of-plane registration algorithm is based on an optimization process, aiming to minimize height differences between overlapping regions of the elementary maps to register. While more precise, this method is also more computation-intensive and therefore slower than the previous methods. Optimization processes are categorized as local optimization, which is fast but may stop optimization in a local extremum, and global optimization, which is slower but more accurate as the global extremum can be reached without interference from local extremum.

Simulated annealing is a global optimization algorithm that can be used for stitching. Van Laarhoven and Aarts [65] explain this algorithm with an analogy
<table><caption>**Table 3.** Benchmark of some in-plane registration algorithms for various surface types. For all metrics, lower is better. Data are gathered from Lemesle *et al* [14].</caption><thead><tr><th>Test case</th><th>In-plane registration algorithm</th><th>Run time [s]</th><th>Memory usage [Mb]</th><th>Mean repositioning error [Pixel]</th><th>Stitching error estimator (SEE) [-]</th></tr></thead><tbody><tr><td rowspan="5">Polymer abrasion with homogeneous morphology</td><td>Geometric</td><td>0.47</td><td>15.25</td><td>NA</td><td>NA</td></tr><tr><td>Cross-correlation based</td><td>0.41</td><td>15.25</td><td>1.4</td><td>-0.47</td></tr><tr><td>Global optimization of differences</td><td>0.51</td><td>15.31</td><td>8.01</td><td>-0.36</td></tr><tr><td>Naive reflectance</td><td>0.5</td><td>20.1</td><td>1.37</td><td>-0.47</td></tr><tr><td>Exhaustive enumeration reflectance</td><td>2.64</td><td>25.57</td><td>1.37</td><td>-0.47</td></tr><tr><td rowspan="2">Independent Multiphysics mirror-polished copper rod with heterogeneous morphology and scale</td><td>Two-step reflectance</td><td>0.76</td><td>17.42</td><td>1.37</td><td>-0.47</td></tr><tr><td>Geometric</td><td>18.37</td><td>400.06</td><td>NA</td><td>NA</td></tr><tr><td rowspan="7">Arm skin replica with large and distinctive features</td><td>Cross-correlation based</td><td>19.31</td><td>523.78</td><td>579.25</td><td>-0.54</td></tr><tr><td>Global optimization of differences</td><td>9.9</td><td>400.88</td><td>16.5</td><td>-0.58</td></tr><tr><td>Naive reflectance</td><td>20.73</td><td>477.32</td><td>10.4</td><td>-0.82</td></tr><tr><td>Exhaustive enumeration reflectance</td><td>80.29</td><td>407.92</td><td>10.85</td><td>-0.82</td></tr><tr><td>Two-step reflectance</td><td>24.12</td><td>401.59</td><td>10.85</td><td>-0.82</td></tr><tr><td>Geometric</td><td>28.7</td><td>681.46</td><td>NA</td><td>NA</td></tr><tr><td>Cross-correlation based</td><td>14.62</td><td>714.05</td><td>44.13</td><td>-0.61</td></tr><tr><td rowspan="5">Titanium alloy abrasion with homogeneous surfaces and more than one hundred elementary maps</td><td>Global optimization of differences</td><td>38.99</td><td>686.99</td><td>33.67</td><td>-0.88</td></tr><tr><td>Naive reflectance</td><td>44.76</td><td>741.18</td><td>12.71</td><td>-1.17</td></tr><tr><td>Exhaustive enumeration reflectance</td><td>176.18</td><td>686.9</td><td>12.74</td><td>-1.17</td></tr><tr><td>Two-step reflectance</td><td>51.82</td><td>686.9</td><td>12.74</td><td>-1.17</td></tr><tr><td>Geometric</td><td>24.19</td><td>223.11</td><td>NA</td><td>NA</td></tr><tr><td rowspan="5">Titanium alloy abrasion with homogeneous surfaces and more than one hundred elementary maps</td><td>Cross-correlation based</td><td>20.98</td><td>222.72</td><td>21.55</td><td>-0.96</td></tr><tr><td>Global optimization of differences</td><td>33.38</td><td>216.69</td><td>63.66</td><td>-0.82</td></tr><tr><td>Naive reflectance</td><td>32.46</td><td>219.44</td><td>2.92</td><td>-1.56</td></tr><tr><td>Exhaustive enumeration reflectance</td><td>99.07</td><td>214.8</td><td>2.98</td><td>-1.56</td></tr><tr><td>Two-step reflectance</td><td>38.89</td><td>214.99</td><td>2.98</td><td>-1.56</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Illustration of piston, tilt and tip between adjacent subaperture measurements. Stitching of 3D topographies aims at registering, i.e., aligning adjacent elementary maps, by finding the correct geometric transformations to apply on one map so that overlapping regions present less differences between each other. Reproduced with permission from [64]. © (2008) COPYRIGHT Society of Photo-Optical Instrumentation Engineers (SPIE).</figcaption></figure>

to the annealing processing of solids. A faster version of this algorithm was proposed by Szu and Hartley [66]. Numerous combinatorial problems were solved with simulated annealing by Kirkpatrick et al [67]. A silhouette-based registration and stitching algorithm based on simulated annealing was proposed by Lensch et al [68]. Another global optimization algorithm is the genetic algorithm, which is designed after natural selection of genes, starting from an initial population of solutions that will be iteratively optimized by keeping elite solutions (best solutions for an iteration), merging and mutating the parameters to optimize, as explained by Holland [69], Forrest [70], Davis [71] and Gendreau and Potvin [72]. Mühlenbein [73] detailed local search based on genetic algorithms. Image stitching and registration using genetic algorithms were used by Xu and Sun [74], Yetis et al [75], Wang and Xu [76] and Prejmerean [77]. Particle swarm optimization can also be of interest and was described by Kennedy and Eberhart [78], Shi and Eberhart [79] and Bratton and Kennedy [80]. Poli et al [81] wrote an overview about the particle swarm algorithm. A variant of the particle swarm algorithm was proposed by Anantathanavit and Munlin [82], called the radius particle swarm optimization. Zhu et al [83] stitched 3D topographies of cylindrical surfaces using a particle swarm optimization-based algorithm. Optimization-based out-of-plane registration still faces difficulties with surface periodicity as multiple best solutions might exist, such as in sinus surfaces, waffers and measurement marbles.

### 3.6. Blending

Next step of the stitching procedure is the blending, which aims at merging elementary maps into a stitched map and to equalize values in the overlapping regions to avoid the appearance of seams, as explained by Ma et al [84]. Binary and fuzzy logic masks are applied to

each element of the overlapping regions of the elementary maps to blend. Blending is necessary as false structures might be created by stitching due to height differences in overlapping regions, which are always present due to noise and measurement uncertainties. Therefore, it is impossible to measure exactly the same topography twice.

Blendless stitching consists in assembling each elementary map without using a blending and likely results in the appearance of seams on the stitched map. It is fast and might be sufficient if morphological structures are much larger than noise. However, blendless stitching is not recommended as the seams might interfere with the computation of segmentation, such as watershed and Wolf pruning segmentation, or might impede on topography interpretation.

Mean blending can alternatively be used: each element in the overlapping regions of the stitched map is computed by the means of values of the elements in each elementary map (Mooney and Stahl [85]). Although slower than blendless stitching, seams do not appear if the vertical drift has been fully compensated. However, if the vertical drift is not appropriately corrected, a mean blending generates three-level stairs: one on the first elementary map, one on the second elementary map and a mean level in the overlapping regions. Mean blending is the most used 3D topography stitching blending in measurement apparatus.

Another approach is to use weights, between 0 and 1, on each element in the overlapping regions depending on the measurement quality of the element. Generally, weights are higher for element closer to the center of the elementary map and lower at the border in optical topography measurement system. Weight values are then dependent on the spatial position of each element inside the elementary map (Pravenaa [46]). This type of blending is called alpha blending. Alpha blending, under the name of fade-in fade-out
algorithm, was used in 3D topography stitching by Wang et al [6] to avoid seam appearance.

A variant of the alpha blending is feathering, which applies a gaussian blur on the weight mask before blending. Feathering is interesting to avoid ghosting, i.e., object or structure duplication due to unexpected movements during acquisition. Feathering is more used in image stitching than in topography stitching, as surface stability is searched in topography measurement in laboratories. The drawback of feathering is the decreasing in resolution for high frequency signal.

Another type of blending, called biband blending, separates low frequency and high frequency signals before blending, so that frequency bands are merged separately before adding resulting merged maps together. A generalized variant of biband blending is multiband blending, which uses the same concept but with more frequency bands.

Poisson blending, also called Poisson editing, uses a different approach compared to the previously described blending algorithms. Poisson editing determines the smoother assembling of each structure in the overlapping regions by resolving the Poisson equation (equation (1), from Szeliski et al [86]) as explained by Dessein et al [87].

$$E_1 = \sum_{i,j} s_{i,j}^x [f_{i+1,j} - f_{i,j} - g_{i,j}^x]^2 \\ + s_{i,j}^y [f_{i,j+1} - f_{i,j} - g_{i,j}^y]^2 \quad (1)$$

Where $f_{i,j}$ is the result image after Poisson editing, $g_{i,j}^x$ and $g_{i,j}^y$ are the target gradient values, $s_{i,j}^x$ and $s_{i,j}^y$ are the gradient constraint weights representing smoothness weights, based on the notation from Szeliski [88].

Poisson blending was used in an image blending algorithm developed by Szeliski et al [86] to limit seam appearance. The proposed algorithm is a faster version of the Poisson blending algorithm, using less memory and less computation intensive.

Lastly, graph cut blending searches for the best frontier between overlapping regions of elementary maps so that a more coherent merging is obtained. This method involves once again binary masks, and the blending quality is linked to the ease of segmentation of the overlapping regions of the elementary maps.

Table 4 groups topographies blending with various algorithms. Vertical drift was not fully compensated purposely, so that the effect of blending on the overlapping regions are more visible.

A summary table comparing the previously described blending algorithms is provided in table 5, based on the blending power, i.e., the ability to merge elementary maps without seams, the error sensitivity, the speed and the conservation of surface fidelity after blending.

Alpha blending is then a good choice for 3D topography stitching as it is both fast, with a good blending power and a good surface fidelity after blending. Graph cut is even better regarding surface fidelity but

it is slower. Most blending algorithm variants, such as feathering or multiband blending, are designed to solve specific defects, such as ghosting. Finally, Poisson blending is interesting in its approach and is suitable for image stitching, but its bad surface fidelity makes it not suitable for 3D topography stitching, which requires a good metrology and fidelity.

### 3.7. Stitched topography

After all steps of the stitching procedure are complete, a stitched map is obtained. A 3D stitched topography of a polymer rod abraded with a grade 80 SiC paper (FEPA) is shown in figure 8. Elementary maps were measured with a white light interferometer Zygo NewView™ 7300 equipped with a 50x objective lens. A larger field of measurement is provided by stitching, while keeping a high resolution, allowing the acquisition of large morphological structures and smaller details with the same measurement gauge.

### 3.8. Data classes in 3D topography stitching

Stitching of topographic elementary maps can use various data classes inside the 3D topography stitching algorithm. As different profilometry techniques exist, the data classes available for stitching may vary. First, metadata (i.e., data about the elementary dataset) of the elementary maps can be useful in stitching, such as absolute position or relative position inside the stitching grid of each elementary maps. Metadata can be sufficient if good enough. This data class is also useful for prepositioning before applying more precise stitching algorithm. Metadata-based stitching algorithm precision depends on the robustness, the accuracy, and the vertical and horizontal resolutions of the profilometer. Markers can also be of interest when using fiducialization, and marker-based registration can be efficient, provided that markers or marker patterns are unique. Alternatively, 3D topographic data can use raw topographic data where overlapping regions of elementary maps are compared to perform in-plane and out-of-plane registrations. It is a fast method, although not as fast as metadata-based stitching. But it is also possible to facilitate the stitching procedure by preprocessing elementary maps. In this case, data class is then transformed topographic data, with operations such as filtering, segmentation or thresholding. It is also possible to compute other map types from topographic such as reflectance (Lemesle et al [89]). The last data class is multimap data, i.e., multiple maps of different types, such as image and topography, aligned element by element. In-plane registration can then use one type of data while out-of-plane registration uses topographies. Note that multimap data classes can be used at the same time depending on the stitching algorithm. Lemesle et al [14] proposed three variants of a reflectance-based multimap stitching algorithm, called the rising sun stitching algorithm. Some variants of this algorithm
<figure><img src="image_11.png" /><figcaption>**Table 5.** Comparison of blending algorithms. Selecting a blending algorithm is a compromise between blending power and error sensitivity on one hand and speed on the other hand.</figcaption></figure>

<table><thead><tr><th>Method</th><th>Blending power</th><th>Error sensitivity</th><th>Speed</th><th>Surface fidelity</th></tr></thead><tbody><tr><td>Blendless</td><td>—</td><td>—</td><td>+++</td><td>-</td></tr><tr><td>Mean</td><td>+</td><td>+</td><td>+++</td><td>+</td></tr><tr><td>Alpha</td><td>+++</td><td>++</td><td>++</td><td>+++</td></tr><tr><td>Feathering</td><td>+++</td><td>++</td><td>+</td><td>++</td></tr><tr><td>Biband</td><td>++</td><td>++</td><td>-</td><td>++</td></tr><tr><td>Multiband</td><td>+++</td><td>-</td><td>-</td><td>++</td></tr><tr><td>Poisson</td><td>+++</td><td>+++</td><td>-</td><td>+</td></tr><tr><td>Graph cut</td><td>+++</td><td>++</td><td>+</td><td>+++</td></tr></tbody></table>

add contrast optimization of reflectance maps before applying in-plane registration algorithm. Out-of-plane registration is then performed using original topographic maps, from which reflectance maps were computed. Table 6 shows a summary of the different data class that can be used in stitching algorithm.

A compromise between stitching registration precision and speed of the stitching algorithm exists as

more accurate stitching algorithms require more operations or more complex data classes. It is also common to use metadata for prepositioning before using another data class for a finer registration. Memory usage also increases when more complex data or multiple maps are used.

# 4. Evaluating errors in 3D topography measurements and stitching

## 4.1. Errors and sources

Stitching of 3D topographies is prone to errors at multiple steps: either due to single measurement error sources, either due to stitching procedure error sources as explained by Wyant and Schmit [90]. Single measurement error includes optical errors, such as optical distortions inherent to optical systems, which complexity increases with the use of more lenses inside the optical system. Optical distortions are said to be in pincushion, if the distortion points towards the center of the image, or in barrel if the distortion is towards the external part of the image. As each lens has its own
<figure><img src="image_12.png" /><figcaption>**Table 6.** Comparison of data classes for stitching. More complex data classes can improve stitching accuracy at the cost of speed and memory consumption. Also, it is possible to combine multiple data classes to get higher accuracy with compatible stitching algorithms.</figcaption></figure>

<table><thead><tr><th>Data class</th><th>Speed</th><th>Accuracy</th><th>Memory size</th></tr></thead><tbody><tr><td>Metadata</td><td>+++</td><td>Apparatus dependent</td><td>+++</td></tr><tr><td>Markers</td><td>++</td><td>+++</td><td>++</td></tr><tr><td>Raw topographic</td><td>++</td><td>++</td><td>++</td></tr><tr><td>Transformed topographic</td><td>+</td><td>+++</td><td>++</td></tr><tr><td>Multimap</td><td>+</td><td>+++</td><td>-</td></tr></tbody></table>

distortion, due to lens characteristics and imperfection, a complex combination of pincushion and barrel distortion might appear, which sometimes is called a moustache distortion. These distortions are mathematically defined by the Brown-Conrady distortion model (equations (2) and (3)) based on the works of Conrady [91] and Brown [92]:

$$ 
\begin{aligned}
x_u = & x_d + (x_d - x_c)(K_1 r^2 + K_2 r^4 + \dots) \\ 
& + (P_1(r^2 + 2(x_d - x_c)^2) + 2P_2) \\ 
& \times (x_d - x_c)(y_d - y_c)(1 + P_3 r^2 + P_4 r^4 + \dots)
\end{aligned} 
\quad (2) 
$$

$$ 
\begin{aligned}
y_u = & y_d + (y_d - y_c)(K_1 r^2 + K_2 r^4 + \dots) \\ 
& + (2P_1(x_d - x_c)(y_d - y_c) + P_2) \\ 
& \times (r^2 + 2(y_d - y_c)^2)(1 + P_3 r^2 + P_4 r^4 + \dots)
\end{aligned} 
\quad (3) 
$$

Where $(x_d, y_d)$ is the point of the distorted image as projected on image plane using a given lens, $(x_u, y_u)$ is the point of the undistorted image as projected by an ideal pinhole camera, $(x_c, y_c)$ is the center of the distortion, $K_n$ is the n-th radial distortion coefficient, $P_n$ is the n-th tangential distortion coefficient and $r = \sqrt{(x_d - x_c)^2 + (y_d - y_c)^2}$ is the Euclidean distance between the point in the distorted image point and the distortion center. In pincushion distortion, $K_1$ is a positive value, while in barrel distortion, $K_1$ is a negative value. King [42] explained how to observe and quantify distortion in topography with the help of diamond markers. As the amplitude of distortion is higher in the external part of the image, each image can be cropped to reduce distortion.

Vignetting is also an optical phenomenon which occurs when the light passes through a lens. As reflection modifies the light path, a part of the reflected light beams does not reach the sensor as explained by Pruss et al [93]. The light intensity then decreases on the external part of the image, as shown in the work of Ma et al [84]. Vignetting can cause topographic differences compared to the
true surface and can decrease the signal-to-noise ratio. The vignetting correction is either based on a combination of algorithms, either on cropping of the region-of-interest of the sensor or on a complex combination of lenses in the objective lens. Note that optical distortion and vignetting can be limited with the use of high-quality objective lenses, but this comes at a high price.

Time-dependent errors are also present as both the single measurement acquisition and the stitching acquisition process of multiple images, including table displacements, are time-dependent process. Thermal deviation is a time-dependent error source in surface topography and microscopy. Small scales can be modified by the material thermal expansion. Thermal vibration can also appear if thermal stability is not reached. This can be observed in the interferometry fringes. Thermal vibrations are generally low frequencies. The differences of height in overlapping regions due to thermal deviation and vibration can infer with the stitching of elementary maps. A controlled room limits the thermal influence during acquisition. Thermal deviation was observed during the image acquisition of an electron beam lithography, which offsets the positions of elementary maps by 10 nm (Paul et al [41]). Surface contamination is also a time-dependent error source, as dust or humidity droplets can contaminate the surface during the stitching acquisition process, especially for large stitching that can last several hours or days. Once again, a clean controlled room in temperature and humidity reduces this error source. More generally, the time-stability of the surface to be measured is essential for stitching, especially when measuring the smaller scales.

Positioning errors are linked to the precision and robustness of the moving axes of the measurement system. At the beginning of a stitching acquisition process, command positions of each elementary maps are calculated for X, Y and Z axes. The measurement system moves from position to position, and acquires an elementary map at each position. However, positioning errors can occur. If this is not corrected by the stitching procedure with a suitable algorithm, in-plane positioning errors can result in ghosting, i.e., a structure duplication which does not exist on the true surface. In the case of out-of-plane positioning errors, a vertical drift is produced, creating a vertical discontinuity in the stitched map. Once again, a suitable stitching algorithm is necessary to avoid vertical drift. Note that stitching algorithms by design aims at compensating positioning errors during the registration steps. However, more complex and time-costly stitching algorithms are needed when positioning errors become high, and in some extreme cases, stitching might not possible.

The second category of error sources is linked to the stitching procedure. First, as stitching requires a computer, analogic data acquired by the sensor are turned into digitized data, which precision is not infinite. This results in height value rounding, which impedes on correct registrations of elementary maps. Although rounding errors mainly interfere with the stitching of very small

scales, these errors are systematically present in each measured dataset. A second problem is the creation of seams on the stitched map (Ma et al [84]). Seams are visible lines caused by value differences in the overlapping regions. In a measured dataset, measurement uncertainties make each overlapping region slightly different and stitching algorithm should take these differences into account during registration. Also, blending is used to merge all elementary maps with each other's without creating seams in the process. Seams should not be present in the stitched map as surface fidelity is decreased and false morphological structures might be generated. Then, roughness parameters, especially segmentation-based parameters might be modified by the presence of seams. Last kind of stitching error source is absent data in the stitched dataset (Ma et al [84]), resulting in bad registration. In the case of topography, these absent data are non-measured points and are common. A 3D topography stitching algorithm needs to be compatible with non-measured points to be of interest to metrologists.

Uncertainty propagation in stitching algorithm is limited. First, most preprocessing operations and registration steps are applied on proxy maps of the original elementary maps, only used to facilitate computation of registration. Thus, these operations do not propagate uncertainty. However, the stitched map is assembled in the blending step from the original elementary maps. The blending operation is the source of uncertainty propagation in stitching. Mathematically, blending operation is described in equation (4).

$$a_{ij} = \sum_{n=1}^{N} w_{n,ij} a_{n,ij} \quad \text{where} \quad \sum_{n=1}^{N} w_{n,ij} = 1 \\ \forall w_n \in [0; 1], \forall i, j \in \mathbb{N}^{+2}, \forall N \in \mathbb{N}^{+} \quad (4)$$

Where $a_{ij}$ is the value of the stitched map A at the i-th row and the j-th column, $a_{n,ij}$ is the value of the elementary maps $A_n$ at the i-th row and the j-th column, $w_{n,ij}$ is the weight to apply on each value of $a_{n,ij}$ at the i-th row and the j-th column. The values of $w_{n,ij}$ are deterministic and depend on the selected blending algorithm. In non-overlapping regions, $w_{n,ij} = 1$ for only one elementary map and other elementary maps weight are $w_{n,ij} = 0$. The uncertainty propagation for each element $a_{ij}$ is then (equation (5)):

$$\sigma_{ij} = \sum_{n=1}^{N} w_{n,ij}^2 \sigma_{n,ij}^2 \\ + \sum_{n=1}^{N} \sum_{m(m \neq n)}^{N} w_{n,ij} \\ w_{m,ij} \rho_{ij} \sigma_{n,ij} \sigma_{m,ij} \\ \text{where} \quad \sum_{n=1}^{N} w_{n,ij} = 1 \quad \forall w_n \in [0; 1], \forall i, \\ j \in \mathbb{N}^{+2}, \forall N \in \mathbb{N}^{+} \quad (5)$$

Where $\sigma_{ij}$ is the standard deviation of element $a_{ij}$ of the stitched map A, $\sigma_{n,ij}$ is the standard deviation of element $a_{n,ij}$ of the elementary maps $A_n$ and $\rho_{ij}$ is the correlation coefficient between two $\sigma_{n,ij}$.
## 4.2. Quantification of stitching error

Stitching error quantification is difficult for two reasons. First, the optimization processes in the stitching registration steps aim at minimizing differences (or maximizing similarities) between overlapping regions, and secondly, blending hides defects after elementary map assembling. However, some attempts to evaluate stitching errors exist, either at the local level, either at the global level. Root mean square (RMS) differences of overlapping regions is a first approach to local error quantifications when computed for each elementary stitching. This approach is questionable to evaluate errors as RMS differences are often used as an objective function in registration. Therefore, it is not a good value for error quantification as it is minimized during the stitching procedure. Also, RMS does not detect all stitching effects and defects due to averaging and blending. Furthermore, fine textures and non-periodic features cannot be well described by RMS as this metrics is not local enough.

A quality predictor was built by Tian et al [94] based on support vector regression and radial basis function that uses quality vectors determined from local features and global features extracted with the SIFT descriptor. The values in the quality vectors from this paper are the distance of ghosting objects, the out-of-plane misalignment measure, the overall structural distortion and the geometric error (defined as the phase angle between line pairs in this work).

A third technique for the local stitching error evaluation is the use of cross-correlation between overlapping maps to detect the number of pixels to move in X and Y translations, respectively $\Delta_{x,offset}$ and $\Delta_{y,offset}$, in order to align surface in each local stitching (Lemesle et al [14]). In-plane repositioning error is then computed with the norm of these values, as in equation (6). This technique is mainly for evaluating in-plane repositioning errors.

$$ Error = \sqrt{\Delta_{x,offset}^2 + \Delta_{y,offset}^2} \qquad (6) $$

If the hypothesis of linearity of the cumulated errors in stitching is assumed, since measurement apparatuses have bounded uncertainties in repositioning errors, i.e., tolerance, cumulated errors should be linear statistically. A linear regression is then used to evaluate how linear the accumulation of errors is. The Stitching Error Estimator (SEE) is then calculated with equation (7), where $p$ is the p-value of the linear regression and $N$ is the number of elementary maps to stitch (Lemesle et al [14]). The last two error quantification values present the advantage to avoid the use of a reference surface to assess the stitching quality.

$$ SEE = \frac{\log_{10}(p)}{N} \qquad (7) $$

It is also possible to compute summary statistics on local errors, such as means or standard deviation, to find if errors are almost constant or if errors increase, such as in stitching chain reaction failure. Note that

the stitching procedure can only be as accurate as the measurement apparatus. Uncertainties in confocal microscopy and apparatus geometric errors are evaluated by Wang et al [95]. Similarly for focus variation microscopy, Pérez et al [96] proposed leads to correct translation errors on X and Y axes.

More methods for stitching error evaluation are proposed at the global level and are mostly based on the use of a known reference surface that is split in elementary maps and stitched together. The stitched map is then compared to the reference surface so that differences and errors are detected. Chen et al [9] used simulated surfaces for this purpose. Peak-to-valley height differences and root mean square of height differences are often used for the comparison, as explained by Yan et al [2]. Gao et al [97] simulated interferometry stitching and computed residual errors, i.e., differences between the stitched topography and the ideal topography in an element-wise fashion. The root mean square difference is then computed on residual errors for stitching error quantification. However, problems arise with this approach. First, the pixels of the stitched surface and the reference surface might not be aligned, and local peaks might be mismatched, and so, differences are not relevant for error evaluation. Secondly, these values are often used as objective function in registration and therefore are not suitable to correctly assess stitching errors as these values are minimized. Also, the root mean square of height differences can sometimes be smaller than the vertical resolution and can be due to uncertainties for small scales and signal noise, as explained by Huang et al [15]. Reference surface can also be a physical object such as mirrors in optics. Rommeveaux and Barrett [98] compared the measured radius of curvature of a mirror with interferometry and stitching, to the nominal radius of curvature of the mirror. Marinello et al [13] calibrated stitched topographies by introducing scaling factors along X and Y axes, called $C_x$ and $C_y$, as tools for stitching performance evaluation (Bariani [99]). Scaling factor $C_x$ is computed with equation (8), where $l_x$ is the side length of a repeated structure in a reference surface, such as a quadratic cell, along X-axis, and $l_{sx}$ is the side length of a repeated structure in a stitched surface. Respectively, the scaling factor $C_y$ is computed for Y-axis. If $C_x = C_y = 1$, registration is considered exact.

$$ C_x = \frac{l_x}{l_{sx}} \qquad (8) $$

Uncertainty in stitched topographies can also be evaluated with more complex statistical analysis method. Chen et al [100] proposed the use of Maximum Likelihood Estimation (MLE) to evaluate stitching error in cylindrical surfaces measured by interferometry.

Lastly, comparison of roughness parameters (such as in ISO 25178-2 standard [4]) between a non-
<table><caption>**Table 7.** Degrees of freedoms of body models for 3D topography stitching. More degrees of freedom make registration more difficult, and sometimes, unsolvable in a reasonable time.</caption><thead><tr><th>Model</th><th>Number of Degrees of Freedom (DoF)</th></tr></thead><tbody><tr><td>Deformable body</td><td>DoF<sub>unique</sub> * N * n * m</td></tr><tr><td>Affine deformable body</td><td>DoF<sub>unique</sub> * N * n * m</td></tr><tr><td>Rigid body</td><td>DoF<sub>unique</sub> * N</td></tr></tbody></table>

stitched and a stitched surface with a higher magnification objective is proposed by Pawlus *et al* [101]. While parameters such as arithmetical mean height Sa, root mean square height Sq, autocorrelation length Sal, texture aspect ratio Str, are stable between both topographies, stitching has a higher influence over other roughness parameters, with higher deviations on latter roughness parameter values.

## 5. Mathematical modelling of 3D stitching

The stitching problem in 3D topography is decomposed in the in-plane and out-of-plane registrations, i.e., the determination of the transformations to apply to each elementary map to obtain a stitched map. Once these transformations are determined, it is interesting to understand how they are applied to each elementary map.

Three main models exist for 3D topography stitching: deformable body, affine deformable body and rigid body. The more general model is the deformable body, where each elementary map can be deformed pixel by pixel independently to fit the overlapping region of another map. The deformable body model is valid as the optical distortion in an optical system and the probe thermal dilation due to friction in contact profilometers generate local deformations in the measured map compared to the true surface. This model is however hard to solve as a high number of degrees of freedom (DoF) are involved (see table 7), where the total number of degrees of freedom *DoF* depends on the number of allowed degrees of freedom per map element *DoF<sub>unique</sub>* , the number of elementary maps to stitch *N*, and the number of rows *n* and columns *m* in each map. As an example, stitching four elementary maps composed of 640 rows and 480 columns for 6 degrees of freedom (XYZ-translations and XYZ-rotations) involved a total number of 7 372 800 degrees of freedom.

A second model is the affine deformable body, since deformations are linked to physical phenomena and therefore are not random but follow laws (similarly to physical laws). Interactions between these complex physical phenomena increase the complexity of these laws. The objective is then to determine these laws and to correct the map distortion. In the hypothesis of affine deformable body, the laws are affine, i.e., the distance between points might vary, and a common transfer function for deformation exists for all

points in the map. Each pixel can move, but all pixels follow a similar law, such as in homothety. The number of degrees of freedom is the same as in the deformable body model, but the correction is simpler to determine and to apply on each elementary map. Affine deformable body model is more justified in image stitching as multiple images can be captured with a camera in a rotational movement. Image homography might require correction before and during the stitching procedure to compensate perspective. While in a laboratory setting, 3D topography stitching often relies on a normal acquisition between a specimen and the measurement system, and the displacement tables more often propose only transition movements for the displacement over the stitching grid, and sometimes add rotational axes. Perspective influence is then smaller than in image. Angle of steep slopes might change with perspective. Also, rectangular features might change into trapezoid features. In-plane orientation of morphological structures and form factor are modified by perspective. Optical distortions can also be included in the affine deformable body. Straight lines in topography become curved with distortions and motif circularity gets lower. Note that deformation correction should be preferably applied before stitching if some deformation laws and parameters are known, such as optical distortion correction depending on optics parameter. Otherwise, stitching should consider the affine deformable body model if such deformations are not negligible. If surface curvature is non-uniform, the non-linear nature of distortion over an elementary map and perspective make the true surface curvature difficult to retrieve.

Last simplified model is the rigid body model, as described by Kong and Chen [10] and Lee *et al* [102]. In the rigid body model, each point has a constant in-plane distance between each other's, and all points are geometrically transformed with the same correction for all degrees of freedom. The total number of degrees of freedom is then drastically reduced. While the general rigid body model uses 6 degrees-of-freedom (XYZ-translations and XYZ-rotations), such as in Chen *et al* [17], most 3D topography stitching algorithms focus on 3 degrees-of-freedom (XYZ-translations) as explained by Lu *et al* [20]. The following sections assume the rigid body model for the stitching problem.

A surface dataset can be represented under two forms. First, a surface data can be represented as a 2D matrix (i.e., a map) where each element of the matrix *A* is a height value *z*<sub>ij</sub> along *n* rows and *m* columns (equation (9)). The coordinates along X and Y axes are implicit with the position of each element inside the matrix and the distance between the points is obtained by multiplying by the lateral resolution. The second form is a *p*<sup>*n*</sup> 3 matrix, where *p* is the total number of pixels. In the case of a conversion from a 2D matrix, *p* = *n*\* *m* pixels. Each column represents an axis, such as X, Y and Z axes (equation (10)). The values *x*<sub>*i*</sub>, *y*<sub>*i*</sub>
<table><caption>**Table 8.** Elementary operations for in-plane and out-of-plane registrations in 3D topography stitching. Note that in-plane operations are based on index shifting and interpolation while out-of-plane operations modify values in the elementary maps.</caption><thead><tr><th>Transformation</th><th>Transformation type</th><th>Elementary operation</th></tr></thead><tbody><tr><td>X-translation</td><td>In-plane</td><td><b>B</b>(i + Δ<sub>x</sub>, j + Δ<sub>y</sub>) = <b>A</b>(i, j) ∀ i, j ∈ ℝ<sup>+</sup></td></tr><tr><td>Y-translation</td><td>In-plane</td><td></td></tr><tr><td>Z-translation</td><td>Out-of-plane</td><td>f(<b>A</b>, γ) = <b>A</b> + γ</td></tr><tr><td>X-rotation</td><td>Out-of-plane</td><td>g(<b>A</b>, α, β) = <b>A</b> + α<b>X</b> + β<b>Y</b></td></tr><tr><td>Y-rotation</td><td>Out-of-plane</td><td></td></tr><tr><td>Z-rotation</td><td>In-plane</td><td>h(<b>A</b>, θ) = rot<sub>θ</sub>(<b>A</b>)</td></tr></tbody></table>

and $z_i$ as respectively the values of the i-th pixel of map A on the X, Y and Z axes. This representation is called a point cloud. The length corresponds to the number of points in the map. The point cloud allows multiple points to share the same XY positions with different height values but needs more memory in computer and is slower to compute compared to the 2D matrix form. However, it is useful to go from one form to another for some transformation operations.

$$A = \begin{pmatrix} z_{1,1} & \cdots & z_{1,m} \\ \vdots & \vdots & \vdots \\ z_{n,1} & \cdots & z_{n,m} \end{pmatrix} \qquad (9)$$

$$A = \begin{pmatrix} x_1 & \cdots & y_1 & \cdots & z_1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_i & \cdots & y_i & \cdots & z_i \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_p & \cdots & y_p & \cdots & z_p \end{pmatrix} \qquad (10)$$

Applying transformations to elementary maps is equivalent to applying elementary operations which depend on the transformation to apply. Related to in-plane transformation, the elementary operation is an index shifting in 2D matrix for the XY-translations, and the Z-rotation consists in applying a rotation matrix to a point cloud, after conversion from 2D matrix. In the case of out-of-plane transformations, the elementary operations are based on the addition of a constant value to the height values for Z-translation, and the addition of a plane, under the form of a 2D matrix, for XY-rotations. Table 8 sum ups the elementary operations depending on the transformation to apply. Seo *et al* [18] detailed the elementary operations for the X and Y rotations (tilt and tip) as the determination of adjacent planes to correct.

Where **A** and **B** are overlapping datasets to register, *i* and *j* are map indices, *α* and *β* are plane

During registrations, optimization processes are used to determine the best alignment in overlapping regions of elementary datasets. An objective function is then selected for the optimization processes. Two approaches exist. First, the goal of registration can be to minimize the differences between two overlapping regions. The objective function is then based on least-square minimization (equation (11)).

$$f(\mathbf{A}, \mathbf{B}) = \sum_{i=0}^{m-1} \sum_{j=0}^{n-1} (a_{ij}^* - b_{ij})^2 \quad \forall i, j \in \mathbb{N}^{+2} \qquad (11)$$

Where *f* is the objective function, **A** and **B** are two overlapping datasets, $a_{ij}^*$ are the elements of the geometrically transform dataset to align the overlapping regions of **A** and **B**, *i* and *j* are respectively element indices in the 2D matrix for the X and Y axes.

On the other hand, the aim of registration can be to maximize similarities of overlapping regions. Cross-correlation can then be used as an objective function to determine the best alignment of elementary maps (equation (12), from Rankov *et al* [103]). Cross-correlation is used by Rankov *et al* [103] for image registration. Cross-correlation-based registration algorithms rely on an optimization process to find similarities between overlapping regions. The displacement to apply is determined by using a moving window, extracted from the first map, and a fixed window, extracted from the second map. The moving window is moved along the fixed window, and the cross-correlation is computed for each position of the moving window. When the maximum of cross-correlation is reached, maximum similarity between overlapping regions is found. The correct displacement is then the displacement applied to the moving window so that cross-correlation is maximized.

$$CC = \frac{\sum_{x=0}^{L-1} \sum_{y=0}^{K-1} (w(x, y) - \bar{w})(f(x + i, y + j) - \bar{f}(i, j))}{\sqrt{\sum_{x=0}^{L-1} \sum_{y=0}^{K-1} (w(x, y) - \bar{w})^2} \sqrt{\sum_{x=0}^{L-1} \sum_{y=0}^{K-1} (f(x + i, y + j) - \bar{f}(i, j))^2}} \qquad (12)$$

coefficients respectively for the X and Y axes, $γ$ is the piston coefficient for the Z axis and $rot_θ$ is the operator for the rotation of an angle $θ$.

Where $w(x, y)$ is the value of the pixel at the element indexes *x* and *y*, and $\bar{w}$ is the mean value of the pixels inside the moving window. $f(x + i, y + j)$ refers to the
<table><caption>**Table 9.** Operations to apply on an elementary map depending on the stitching model. As the stitching model uses more dimensions, more elementary operations are combined, and therefore, stitching algorithm complexity is higher.</caption><thead><tr><th>Stitching model</th><th>Transformation operation</th></tr></thead><tbody><tr><td>1D</td><td>b(i + Δ<sub>x</sub>) = a(i) ∀ i ∈ N<sup>+</sup></td></tr><tr><td>2D</td><td>b(i + Δ<sub>x</sub>) = (f ◦ g<sub>x</sub>)(a(i)) ∀ i ∈ N<sup>+</sup><br>= a(i) + αx(i) + γ ∀ i ∈ N<sup>+</sup></td></tr><tr><td>3D</td><td>B(i + Δ<sub>x</sub>, j + Δ<sub>y</sub>)<br>= (h ◦ g ◦ f)(A(i, j)) ∀ i, j ∈ N<sup>+</sup><br>= rot<sub>0</sub>(A(i, j) + αX(i, j) + βY(i, j) + γ) ∀ i, j ∈ N<sup>+</sup></td></tr></tbody></table>

value of the composite image at a given pixel inside the moving window, and $f(i, j)$ is the mean value of all composite image pixels inside the moving window. The window dimensions in pixels are represented by $K$ and $L$.

Assuming the rigid body hypothesis, the stitching problem can be modelled depending on the elementary operations described previously for each transformation. Three stitching models are described in table 9: 1D model (profile), 2D model (profile with tilt and tip) and 3D model (surface topography).

These mathematical models of the stitching problem assume that registration can find an exact match between the overlapping regions of the elementary maps. This case is generally found in simulated surfaces, but never happens with measured maps as uncertainties exist during the acquisition of elementary maps and the stitching procedure. Obtaining two identical overlapping regions in elementary maps with successive measurements is unlikely (Lemesle et al [104] and Lemesle et al [105]). Therefore, the stitching problem is about finding geometrical transformations that minimize differences (or maximize similarities) between overlapping regions of the elementary maps to stitch. Also, non-measured points need to be considered. During the stitching procedure, non-measured points should be considered as wildcard values, that can match any other value in the overlapping regions. Especially, the challenge with non-measured points is not their presence, but their spatial distribution on maps which can be different in each elementary map. In the case of contact profilometers, capped height values can be obtained if heights are out of the Z scanning range. These values are misleading in stitching and should be replaced by non-measured points before registration. Finally, as non-measured points should be considered as wildcards, the number of potential solutions in the optimization processes of registration steps increases, making the stitching algorithm run time longer.

## 6. Conclusion and further works

Topography measurements can rely on stitching, an assembly technique producing a seamless high resolution over a large field of measurement map. The scale range is therefore increased, and both large structures

and small details can be captured on the same map. More scales can be explored with multiscale analysis when stitching is used, which helps the detection of physical phenomena at their scale of application. This article reviews existing stitching algorithms, especially 3D topography stitching algorithms. First, the stitching procedure is described: the elementary map acquisition, the preprocessing, the in-plane and out-of-plane registrations, the blending and the obtention of the stitched map. Especially, the in-plane registration is categorized according to six classes: manual, geometric, fiducialization, global optimization-based algorithm, computer-vision-based direct and feature detection methods. Blending algorithms are also presented. Algorithms are compared according to their precision and their run time. Secondly, measurement and stitching errors, their sources and their evaluation methods are detailed. Evaluating stitching errors is difficult, and only a few articles proposed some methods, mainly with flaws. This part is a lead to investigate in future research. Lastly, the mathematical modelling of the stitching problem is shown, from the description of elementary operations in rigid body models to the optimization process and objective functions used in the registration steps of the stitching procedure. The consideration of non-measured points in 3D topography stitching is also explained. Overall, 3D topography stitching is fundamentally a multidisciplinary field, combining knowledge in topography, metrology, optimization, computer vision and algorithmics. This article is both a review and a tutorial for any metrologist wanting to exploit high resolution over large field of measurement topography, by detailing important notions and algorithms useful for 3D topography stitching and roughness studies.

This review of the literature on 3D topography stitching has presented the multidisciplinary nature of the algorithms and the following points have been categorized depending on the domain:

- **Mechanical analysis:** Scales are of primary importance in roughness studies and the use of stitching allows a wider scale range to be explored during the multiscale analysis of physical phenomena. Creating more robust 3D topography stitching algorithms would allow larger field of measurement to
capture all the richness in morphological signatures of each physical phenomenon separated by scale. An extrem case would be a stitching of the whole surface of an object, but an analysis workflow should be designed, not limited to the statistical analysis of roughness parameters. A segmentation of such large topography, such as with Wolf pruning segmentation, can be used to evaluate heterogeneous surfaces, especially interface;

- Stitching error evaluation: Only a few researchers tried to evaluate stitching errors. It is a difficult matter as most information about stitching errors are lost during blending, which is necessary to achieve a seamless topography but hides errors along the way. Also, pixel-by-pixel comparison is not suitable as the correct pixel alignment cannot be secured. Lastly, stitching aims to minimize errors during registration steps. A statistical approach for the evaluation of in-plane errors is proposed, but stitching out-of-plane registration errors are still to be robustly assessed.

- Algorithmics: Stitching larger surfaces implies a larger run time of the stitching algorithm itself. For very large stitching, algorithmic complexity should be considered. A study of the complexity of the different stitching algorithms would allow to select the most suitable algorithm depending on the nature of the surface morphologies, the expected size of the stitched topography and the allotted run time. Also, it would be interesting to find the maximum number of stitching elementary maps that can be stitched in an allotted time;

- Optimization: Selecting the most suitable stitching algorithm depending on the nature of the morphology is still unexplored. Building a topography database, dedicated to 3D topography stitching algorithm development, would allow the study of the determination of the stitching algorithm. Then, stitching software would be able to select the most adequate algorithm for the whole stitching procedure, or even better, for each local stitching. The latter would improve the ability of a stitching software to stitch heterogeneous topographies, by adjusting the stitching algorithms to morphological heterogeneity;

- Multimap: Using multimap-based 3D topography stitching algorithms can increase the precision of the registration, by using information from topography, image or even reflectance map. This approach should be extended to more types of maps. The design of a multimap measurement apparatus would help both stitching and multimap analysis for a better understanding of physical phenomena.

Also, several challenges will be faced by 3D stitching algorithms. As surfaces are getting more and more

complex with surface functionalization, especially if multiple surface functionalities are piloted on the morphology. Therefore, the number of heterogeneous surfaces or multi-frequency morphological structures can be expected to increase, and 3D topography stitching algorithms would need evolution to correctly handle these cases. On the other hand, stitching algorithm speed would need to be increased, up until real-time stitching in metrological inspection. While acquisition time of elementary maps is not yet real-time, machine learning models inferring topographies from a single image instead of multiple image layers might appear, thus reducing acquisition time drastically. Then, real-time stitching algorithms would be necessary. An increase in algorithm speed would also be profitable for the stitching of large datasets of elementary maps.

Potential avenues for stitching algorithms should focus on stitching accuracy, speed, data ingestion and error correction. Machine learning could then be used in two steps of the stitching procedure. First, machine learning might improve feature detection of topographies, both as a preprocessing or a feature extractor, leading to more accurate and faster registration. Another possibility would be stitching only one image per subregions of the stitching grid, performing the stitching procedure on such images and inferring topography or roughness parameters from the stitched image. This would save acquisition and stitching time drastically, making such stitching strategy closer to real-time stitching in metrology. Also, error correction is a hard challenge, which first requires the design of more relevant stitching error metrics. By introducing better error metrics, determination of the most suitable stitching algorithm depending on surface typology would get easier, as well as designing better error correction and registration algorithms. Finally, most feature descriptors were designed for image registration but are not optimized for 3D topography, as the frequency of high-gradient local regions is far less common in topography than in image. The design of a feature descriptor dedicated to 3D topography would make stitching of large elementary maps more efficient, reducing the need for a large number of elementary maps for the stitching and thus, reducing stitching run time.

## Acknowledgments

This work was supported by the Agence Nationale de la Recherche (ANR) [grant number ANR-22-CE51-0031].

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).
ORCID iDs

Robin Guibert https://orcid.org/0000-0002-2656-2840

Frédéric Robache https://orcid.org/0000-0003-0547-4873

Julie Lemesle https://orcid.org/0000-0001-6119-8725

Christopher A Brown https://orcid.org/0000-0002-9355-4824

Maxence Bigerelle https://orcid.org/0000-0002-4144-245X

References

[1] Yu Y and Peng H 2011 Automated high speed stitching of large 3D microscopic images *2011 IEEE Int. Symp. Biomed. Imaging Nano Macro* 2011 (IEEE) 238–41

[2] Yan L *et al* 2019 Subaperture stitching testing for fine flat mirrors with large apertures using an orthonormal polynomial fitting algorithm *Opt. Lasers Eng.* **120** 49–58

[3] Brown C A *et al* 2018 Multiscale analyses and characterizations of surface topographies *CIRP Ann.* **67** 839–62

[4] ISO 25178-2 2021 Geometrical Product Specification (GPS)-Surface Texture: Areal. Part 2: Terms, Definitions and Surface Texture Parameters. Int Organ Stand Geneva Switz 2021

[5] ASME 2002 B46.1 standard - surface texture, surface roughness, waviness and lay *Am Soc Mech Eng Am Natl Stand* 1–98

[6] Wang R, Li Z-L, Ren M-J and Zhu L-M 2022 A registration-based stitching method for obtaining high-accuracy material removal distribution in the sub-aperture polishing process *Precis. Eng.* **77** 251–62

[7] Laubach S, Ehret G, Riebeling J and Lehmann P 2016 A new form measurement system based on subaperture stitching with a line-scanning interferometer *Adv. Opt. Technol.* **5** 415–22

[8] Pan F, Lu X, Dong B, Ma X and Xiao W 2016 Subaperture stitching interferometry based on digital holography *Opt. Lasers Eng.* **86** 228–35

[9] Chen S, Xue S, Wang G and Tian Y 2017 Subaperture stitching algorithms: a comparison *Opt. Commun.* **390** 61–71

[10] Kong L and Chen S 2020 Subaperture stitching test for large aperture mild acylinders *Opt. Commun.* **455** 124526

[11] Chen S, Xue S, Dai Y and Li S 2015 Subaperture stitching test of large steep convex spheres *Opt. Express* **23** 29047

[12] Chen S, Li S, Dai Y and Zheng Z 2007 Testing of large optical surfaces with subaperture stitching *Appl. Opt.* **46** 3504

[13] Marinello F, Bariani P, Chiffre L D and Hansen H N 2007 Development and analysis of a software tool for stitching three-dimensional surface topography data sets *Meas. Sci. Technol.* **18** 1404–12

[14] Lemesle J, Guibert R and Bigerelle M 2023 A novel 3D topography stitching algorithm based on reflectance and multimap *Appl. Sci.* **13** 857

[15] Huang L *et al* 2019 Two-dimensional stitching interferometry for self-calibration of high-order additive systematic errors *Opt. Express* **27** 26940

[16] O'Donohue S, Murphy P, Fleig J and Devries G 2006 *Stitching Interferometry for Flexible Asphere Metrology Proceedings of the 21st Annual ASPE Meeting (ASPE 2006)* 4

[17] Chen S, Lu W, Guo J, Zhai D and Chen W 2022 Flexible and high-resolution surface metrology based on stitching interference microscopy *Opt. Lasers Eng.* **151** 106915

[18] Seo Y B, Joo K-N, Ghim Y-S and Rhee H-G 2021 Subaperture stitching wavelength scanning interferometry for 3D surface measurement of complex-shaped optics *Meas. Sci. Technol.* **32** 045201

[19] Guibert R, Hanafi S, Deltombe R, Bigerelle M and Brown C A 2020 Comparison of three multiscale methods for topographic analyses *Surf. Topogr.: Metrol. Prop.* **8** 024002

[20] Lu W, Chen S, Zhang K and Zhai D 2022 Characterization of diffractive relief structures over large areas by stitching interference microscopic topography *Measurement* **202** 111850

[21] Deltombe R, Kubiak K J and Bigerelle M 2014 How to select the most relevant 3D roughness parameters of a surface: relevance of 3D roughness parameters *Scanning* **36** 150–60

[22] Deltombe R, Bigerelle M and Jourani A 2015 Analysis of the effects of different machining processes on sealing using multiscale topography *Surf. Topogr.: Metrol. Prop.* **4** 015003

[23] Spichtinger J, Schulz M, Ehret G and Tutsch R 2022 Traceable stitching interferometry for form measurement of moderately curved freeform surfaces *Opt. Eng.* **61** 094102

[24] Leach R (ed) 2011 *Optical measurement of surface topography* (Springer)

[25] Murphy P, Forbes G, Fleig J, Dumas P and Tricard M 2003 Stitching interferometry: a flexible solution for surface metrology *Opt. Photonics News* **14** 38–43

[26] Cimba M *et al* 2004 *Original Optical Metrologies of Large Components* ed R Geyl, D Rimmer and L Wang 322

[27] Mironova A, Robache F, Deltombe R, Guibert R, Nys L and Bigerelle M 2020 Digital cultural heritage preservation in art painting: a surface roughness approach to the brush strokes *Sensors* **20** 6269

[28] Bigerelle M *et al* 2023 Fractal and statistical characterization of brushstroke on paintings *Surf. Topogr.: Metrol. Prop.* **11** 015019

[29] Fu S, Cheng F, Tjahjowidodo T, Zhou Y and Butler D 2018 A non-contact measuring system for *in-situ* surface characterization based on laser confocal microscopy *Sensors* **18** 2657

[30] Montgomery P C, Benatmane A, Fogarassy E and Ponpon J P 2002 Large area, high resolution analysis of surface roughness of semiconductors using interference microscopy *Mater. Sci. Eng.* B **91–2** 79–82

[31] Zhao S and Yang Q 2014 A stitching method for AFM based large scale scanning with high resolution *IFAC Proc Vol 47* 2697–702

[32] Liu M Y, Cheung C F, Cheng C H, Su R and Leach R K 2017 A Gaussian process and image registration based stitching method for high dynamic range measurement of precision surfaces *Precis. Eng.* **50** 99–106

[33] Otsu N 1979 A threshold selection method from gray-level histograms *IEEE Trans. Syst. Man Cybern.* **9** 62

[34] Lei Z *et al* 2016 A novel 3D stitching method for WLI based large range surface topography measurement *Opt. Commun.* **359** 435–47

[35] Kanopoulos N, Vasanthavada N and Baker R L 1988 Design of an image edge detection filter using the Sobel operator *IEEE J. Solid-State Circuits* **23** 358–67

[36] Gao W, Zhang X, Yang L and Liu H 2010 An improved Sobel edge detection 2010 3rd Int. Conf. Comput. Sci. Inf. Technol. (IEEE) 67–71

[37] Zhang Y, Han X, Zhang H and Zhao L 2017 Edge detection algorithm of image fusion based on improved Sobel operator 2017 IEEE 3rd Inf. Technol. Mechatron. Eng. Conf. ITOEC **457–61**

[38] Moretti M, Gambucci G, Leach R K and Senin N 2019 Assessment of surface topography modifications through feature-based registration of areal topography data *Surf. Topogr.: Metrol. Prop.* **7** 025003

[39] Dong B, Xiao W, Pan F and Che L 2018 High-precision spherical subaperture stitching interferometry based on digital holography *Opt. Lasers Eng.* **110** 132–40

[40] King C W and Bibby M 2014 Development of a metrology workstation for full-aperture and sub-aperture stitching measurements *Procedia CIRP* **13** 359–64

[41] Paul P, Knoll A W, Holzner F and Duerig U 2012 Field stitching in thermal probe lithography by means of surface roughness correlation *Nanotechnology* **23** 385307

[42] King C 2008 A new approach to stitching optical metrology data
[43] Kammers A D and Daly S 2011 Small-scale patterning methods for digital image correlation under scanning electron microscopy *Meas. Sci. Technol.* **22** 125501

[44] Rust M A and Todd R I 2007 High resolution surface studies of superplastic deformation *Mater. Sci. Forum* **551–2** 615–20

[45] Xu Z *et al* 2009 Three-dimensional profile stitching based on the fiducial markers for microfluidic devices *Opt. Commun.* **282** 493–9

[46] Pravenaa S 2016 *A Methodical Review on Image Stitching and Video Stitching Techniques* **11** 7

[47] Agier R 2017 Recalage de groupes d'images médicales 3D par extraction de points d'intérêt *Phd Thesis* Université de Lyon

[48] Lowe D G 1999 Object recognition from local scale-invariant features *Proc. Seventh IEEE Int. Conf. Comput. Vis.* 2 (IEEE) 1150–7

[49] Brown M and Lowe D 2002 Invariant features from interest point groups *Procedings Br. Mach. Vis. Conf.* 2002 (British Machine Vision Association) **1**–10

[50] Brown L 2003 Recognising panoramas *Proc. Ninth IEEE Int. Conf. Comput. Vis.* 2 (IEEE) 1218–25

[51] Brown M and Lowe D G 2007 Automatic panoramic image stitching using invariant features *Int J Comput Vis* **74** 59–73

[52] Bay H, Tuytelaars T and Van Gool L 2006 Surf: speeded up robust features *Eur. Conf. Comput. Vis.* (Springer) 404–17

[53] Li L, Liu Q and Huang W 2019 Robust subaperture stitching method based on three-dimensional surface data matching *Opt. Eng.* **58** 1

[54] Valgren C and Lilienthal A J 2007 *SIFT, SURF and seasons long-term outdoor localization using local features* 253–8

[55] Harris C and Stephens M 1988 A combined corner and edge detector *Proc Fourth Alvey Vis. Conf.* 147–52

[56] Leutenegger S, Chli M and Siegwart R Y 2011 BRISK: binary robust invariant scalable keypoints *2011 Int. Conf. Comput. Vis.* (IEEE) 2548–55

[57] Rosten E and Drummond T 2005 Fusing points and lines for high performance tracking *Int. Conf. Comput. Vis.*, Springer 1508–15

[58] Matas J, Chum O, Urban M and Pajdla T 2004 Robust wide-baseline stereo from maximally stable extremal regions *Image Vis Comput* **22** 761–7

[59] Obdržálek D, Basovník S, Mach L and Mikulík A 2010 Detecting scene elements using maximally stable colour regions ed A Gottscheber, D Obdržálek and C Schmidt *Res. Educ. Robot. - EUROBOT 2009* 82 (Springer Berlin Heidelberg) 107–15

[60] Jianbo Shi, T 1994 Good features to track *Proc. IEEE Conf. Comput. Vis. Pattern Recognit. CVPR-94*, Seattle, WA, USA (IEEE Comput. Soc. Press) 593–600

[61] Rublee E, Rabaud V, Konolige K and Bradski G 2011 ORB: An efficient alternative to SIFT or SURF *2011 Int. Conf. Comput. Vis.* 2564–71

[62] Alcantarilla P F, Bartoli A and Davison A J 2012 KAZE features *Comput. Vis.*—ECCV 2012 ed A Fitzgibbon, S Lazebnik, P Perona, Y Sato and C Schmid (Springer) 214–27

[63] Juan L and Oubong G 2010 SURF applied in panorama image stitching *2010 2nd Int. Conf. Image Process. Theory Tools Appl.* 495–9

[64] Zhao C and Burge J H 2008 Stitching of off-axis sub-aperture null measurements of an aspheric surface *Interferom. XIV Tech. Anal.* 7063 (SPIE) 409–15

[65] van Laarhoven P J M and Aarts E H L 1987 Simulated annealing *In: van Laarhoven PJM, Aarts EHL, editors.* *Simulated Annealing Theory Appl.* (Springer) 7–15

[66] Szu H and Hartley R 1987 Fast simulated annealing *Phys. Lett. A* **122** 157–62

[67] Kirkpatrick S, Gelatt C D and Vecchi M P 1983 Optimization by simulated annealing *Science* **220** 671–80

[68] Lensch H P A, Heidrich W and Seidel H-P 2001 A silhouette-based algorithm for texture registration and stitching *Graph Models* **63** 245–62

[69] Holland J H 1992 Genetic algorithms *Sci. Am.* **267** 66–73

[70] Forrest S 1993 Genetic algorithms: principles of natural selection applied to computation *Science* **261** 872–8

[71] Davis L 1991 *Handbook of Genetic Algorithms* (CUMINCAD)

[72] Gendreau M and Potvin J-Y (ed) 2010 *Handbook of Metaheuristics* 2nd ed (Springer)

[73] Mühlenbein H 2003 Genetic algorithms *Local Search Comb. Optim.* ed E Aarts and J K Lenstra (Princeton University Press) 137–72

[74] Xu Y and Sun C 2017 Image stitching method based on genetic algorithm *Adv Eng.* **126** 7

[75] Department of Computer Engineering, Firat University, Elazig, Turkey, Yetis H, Baygin M and Karakose M 2017 A new micro genetic algorithm based image stitching approach for camera arrays at production lines *J Image Graph* **5** 20–4

[76] Wang Y and Xu L 2008 A global optimized registration algorithm for image stitching *2008 Congr. Image Signal Process.* 3, 525–9

[77] Prejmerean V 2006 The stitching of images 2006 IEEE Int. Conf. Autom. Qual. Test. Robot. 2, 386–91

[78] Kennedy J and Eberhart R 1995 Particle swarm optimization *Proc. ICNN95 - Int. Conf. Neural Netw.* 4, 1942–8

[79] Shi Y and Eberhart R C 1999 Empirical study of particle swarm optimization *Proc. 1999 Congr. EComput.-CEC99 Cat No 99TH8406* 3, 1945–50

[80] Bratton D and Kennedy J 2007 Defining a standard for particle swarm optimization *2007 IEEE Swarm Intell. Symp.* 120–7

[81] Poli R, Kennedy J and Blackwell T 2007 Particle swarm optimization: an overview *Swarm Intell* **1** 33–57

[82] Anantathanavit M and Munlin M 2013 Radius particle swarm optimization *2013 Int. Comput. Sci. Eng. Conf. ICSEC* 126–30

[83] Zhu F, Guo P and Chen X 2019 A sub-aperture stitching algorithm for cylindrical surfaces *9th Int. Symp. Adv. Opt. Manuf. Test. Technol. Adv. Opt. Manuf. Technol.* (Int. Society for Optics and Photonics) 1083805

[84] Ma B *et al* 2019 A fast algorithm for material image sequential stitching *Comput. Mater. Sci.* **158** 1–13

[85] Mooney J T and Stahl H P 2005 Sub-pixel spatial resolution micro-roughness measurements with interlaced stitching *In: Duparr½Å* A ed B Singh and Z-H Gu 587801 (https://doi.org/10.1117/12.617651)

[86] Szeliski R, Uyttendaele M and Steedly D 2011 Fast poisson blending using multi-splines *2011 IEEE Int. Conf. Comput. Photogr. ICCP* 1–8

[87] Dessein A, Smith W A P, Wilson R C and Hancock E R 2014 Seamless texture stitching on a 3D mesh by poisson blending in patches *ICIP - 21st IEEE Int. Conf. Image Process.* (https://doi.org/10.1109/ICIP.2014.7025407)

[88] Szeliski R 2006 Locally adapted hierarchical basis preconditioning *ACM Trans. Graph.* **25** 1135–43

[89] Lemesle J, Robache F, Le Goic G, Mansouri A, Brown C A and Bigerelle M 2020 Surface reflectance: an optical method for multiscale curvature characterization of wear on ceramic-metal composites *Materials* **13** 1024

[90] Wyant J C and Schmit J 1998 Large field of view, high spatial resolution, surface measurements *Int. J. Mach. Tools Manuf* **38** 691–8

[91] Conrady A E 1919 Decentred lens-systems *Mon. Not. R. Astron. Soc.* **79** 384–90

[92] Brown D 1966 Decentering distortion of lenses *Photogramm. Eng.* **32** 444–62

[93] Pruss C, Baer G B, Schindler J and Osten W 2017 Measuring aspheres quickly: tilted wave interferometry *Opt. Eng.* **56** 111713

[94] Tian C, Chai X and Shao F 2021 Stitched image quality assessment based on local measurement errors and global statistical properties *J Vis Commun Image Represent* **81** 103324

[95] Wang C, Gómez E and Yu Y 2019 Characterization and correction of the geometric errors using a confocal microscope for extended topography measurement: II. Experimental study and uncertainty evaluation *Electronics* **8** 1217
[96] Pérez P, Syam W, Albajez J, Santolaria J and Leach R 2019 Stitching-free lateral error compensation for focus variation microscopy *Meas. Sci. Technol.* **30** 065002

[97] Gao F, Jiang Z, Zhao Z and Li B 2015 Measurement of aspheric surface combining point diffraction interferometry and annular subaperture stitching *Opt. Eng.* **54** 014102

[98] Rommeveaux A and Barrett R 2010 Micro-stitching interferometry at the ESRF Nucl Instrum Methods Phys. Res. Sect. Accel. Spectrometers Detect Assoc. Equip. **616** 183-7

[99] Bariani P 2005 Dimensional metrology for microtechnology

[100] Chen D, Wang C, Valyukh S, Wu X and Yu Y 2020 Bayesian uncertainty evaluation of stitching interferometry for cylindrical surface *Measurement* **157** 107626

[101] Pawlus P, Reizer R, Wieczorowski M and Krolczyk G M 2023 Study of surface texture measurement errors *Measurement* **210** 112568

[102] Lee S, Choi M, Kim H and Park F C 2007 Geometric direct search algorithms for image registration *IEEE Trans. Image Process.* **16** 2215-24

[103] Rankov V, Locke R J, Edens R J, Barber P R and Vojnovic B 2005 *An Algorithm for image stitching and blending* ed J-A Conchello, C J Cogswell and T Wilson 190

[104] Lemesle J et al 2023 Top-down determination of fluctuations in topographic measurements *Materials* **16** 473

[105] Lemesle J et al 2023 Height fluctuations and surface gradients in topographic measurements *Materials* **16** 5408