Assessment of surface topography modifications through feature-based registration of areal topography data

Surface topography modifications due to wear or other factors are usually investigated by visual and microscopic inspection, and—when quantitative assessment is required—through the computation of surface texture parameters. However, the current state-of-the-art areal topography measuring instruments produce detailed, areal reconstructions of surface topography which, in principle, may allow accurate comparison of the individual topographic formations before and after the modification event. The main obstacle to such an approach is registration, i.e. being able to accurately relocate the two topography datasets (measured before and after modification) in the same coordinate system. The challenge is related to the measurements being performed in independent coordinate systems, and on a surface which, having undergone modifications, may not feature easily-identifiable landmarks suitable for alignment. In this work, an algorithmic registration solution is proposed, based on the automated identification and alignment of matching topographic features. A shape descriptor (adapted from the scale invariant feature transform) is used to identify landmarks. Pairs of matching landmarks are identified by similarity of shape descriptor values. Registration is implemented by resolving the absolute orientation problem to align matched landmarks. The registration method is validated and discussed through application to simulated and real topographies selected as test cases.


Introduction
The quantitative assessment of surface topography modifications as a consequence of wear and/or other damage events is of paramount importance in many tribological studies related to product, process and material characterisation. Currently, the most common way to quantitatively assess topographic modifications is using profile measurement and consists of computing surface texture parameters on profiles measured before and after the modification event. The most frequently used parameter is the mean arithmetic roughness Ra as defined in the specification standard ISO 4287 [1]. Several other ISO 4287 profile parameters, such as Rz, Rq, and the parameters related to the material ratio curve (also described in ISO 4287), are also often considered when studying surface modifications [2].
The use of state-of-the-art areal surface topography measuring instruments to investigate surface modifications generated as a consequence of wear or other damage, has been explored in recent literature [3][4][5][6][7][8]. Areal topography measurement allows the reconstruction of entire portions of the measured surface, and introduces a wide range of opportunities for quantification through the set of areal field parameters, as defined in ISO 25178-2 [9][10][11][12]. Some areal field parameters are a direct evolution of their profile counterparts (e.g. Sa is the areal equivalent of the Ra profile parameter). Other field parameters capture properties which are visible only on areal datasets (e.g. Std-surface texture direction). In addition to providing multiple viewpoints and useful information on how a topography has changed, areal field parameters about a topography, and are typically unable to inform about topographic modifications pertaining to individual features, such as a specific peak being eroded, or a cavity being filled. Such an assessment of localised, quantitative topographic modification could shed further light on wear and other surface modification phenomena. The quantification of shape changes pertaining to individual surface formations is in principle now possible, due to the increasingly accurate metrological rendition of topography achievable by state-ofthe-art areal topography measuring instruments. However, one of the biggest obstacles to such a comparison is the need to correctly relocate the two topographies (i.e. before and after the modification event) in the same coordinate system, which is a precondition for the quantitative computation of several local topography differences, e.g. the computation of volumetric modifications.
Accurate repositioning/relocation of the topography datasets in the same coordinate system, often referred to as registration, is required because areal topography measuring instruments return topography information (height maps, point clouds or triangulated meshes) in their own coordinate system. Repeatable sample fixturing can offer limited help, as the sample datum surfaces may change if the sample is physically modified. Fixture-provided positioning accuracy may also be inadequate, considering the small scales of the topography modifications which often need to be assessed.
In the absence of a reliable, external relocation reference, information directly extracted from the measured datasets may be useful for relocation. The general idea is that the registration problem may be solved by aligning topographic regions that have remained invariant during the surface modification process.
In marker-based alignment, the operator is asked to manually place marker points to establish correspondences. Implicitly, while searching for correspondences, the operator tends to select regions that at least look visually similar, i.e. somewhat invariant across measurements. Once multiple correspondences are set (at least three, non-collinear, to solve registration with six degrees of freedom (6 dof)), the alignment transformation (translation and rotation parameters) can be found by solving the absolute orientation problem [13] (also referred to as the Procrustes superimposition problem in statistical literature [14]). The absolute orientation problem has closed-form solution [13], thus it can typically be solved in a fraction of a second on current hardware. The disadvantage is that markers are placed by an operator, relying on subjective visual assessment, resulting in a process which is often not accurate, traceable or reproducible [5]. Marker placement can be facilitated by particularly visible topographic landmarks, which may have been created on purpose (e.g. micro indentations), or may exist naturally (e.g. pores, scratches or any other visible singularities). For a topographic landmark to be suitable for registration, it should be easily recognisable, and conserved across measurements. However, if the surface modification being studied is subtle and spread across the area under observation, suitable landmarks may not be easily found by simple observation.
An approach to address some of the problems associated with manual marker placement has been recently presented elsewhere [5]: the position of manually placed markers is refined by studying the cross-correlation function in the surroundings of the marker position, the position of the maximum crosscorrelation value corresponding to the optimal shift to maximise local similarity. Limitations of such an approach are: (a) cross-correlation can only operate on the x, y translation plane, thus any rotational misalignment will possibly degrade the performance in the identification of matching regions; (b) because of the inability to handle rotation, local cross-correlation around a marker position is only suitable for small displacements, as larger displacement may generate multiple maxima in the cross-correlation function. Finally, the approach does not remove the issue that an initial marker placement is still required, which again implies that some landmarks should still be visible.
When topographic landmarks useful for registration are not easily identifiable by visual inspection, global alignment can be attempted instead. In global alignment, an alignment transform is sought, which attempts to maximise the degree of overlapping of both the entire topography datasets.
The two most common algorithmic approaches to global alignment are based either on computing the cross-correlation function over the entire datasets, and finding its maximum [8], or on applying the iterative closest point (ICP) algorithm [15,16]. As mentioned earlier, the disadvantage of cross-correlation is that it reduces the search space to translations only. The ICP method on the contrary operates in 6 dof [17]. This provides greater freedom while investigating the solution space. However, the biggest limitation of ICP is that it is a fine alignment solution, i.e. it requires a previous, coarse alignment to operate effectively. This is because ICP operates as a minimisation process, with scarce capabilities to verify premature convergence to local minima [15,17].
The fundamental problem with global alignment is that it tries to maximise overlapping everywhere, including the regions that have undergone modification and should not be aligned. Global alignment thus works primarily in those scenarios where invariant regions are prevalent over modified ones, thus modified regions weight little in the alignment process.
Clearly, a selective alignment solution, based on attempting to maximise the overlapping of invariant regions only, would be preferable. To cope with visual identification challenges and to ensure better repeatability, the identification of suitable matching regions across datasets should be algorithmically driven. However, algorithmic approaches to implement the automated identification and alignment of landmark regions have been scarcely explored in the literature on surface metrology. Some work can be found on computing local topography descriptors within moving windows [18], and on using local topography descriptors to obtain coarse alignment results [16,19]. In this work, a selective alignment approach based on the automated identification of reference landmarks and their alignment is presented, where a transform originally devised for image processing, the scale invariant feature transform (SIFT) [20], is adopted for the identification of landmarks.

Test cases
Two simulated test cases and one real test case were selected. The simulated test cases featured known misalignments, thus could be used to quantitatively assess the accuracy of the registration. The real test case was used to investigate the behaviour of the method when confronted with real-life measured datasets.

Simulated test cases
The simulated test cases consisted of surface topographies algorithmically generated in their original state, and then subjected to simulated modification processes. The two resulting topographies (before and after the modification process) were then subjected to simulated measurement ten times each, with simulated measurement error on local height determination and sample localisation error (position of the field of view on the sample surface). Relocation error between each pair of 'before and after' measurements was obtained by combining the sample localisation errors known from simulation.
The first simulated test case (referred to as test case n. 1) consisted of a step-like feature subjected to visible, localised erosion damage (see figure 1). The second simulated case (referred to as test case n. 2) featured a base topography comprised of a mixture of deterministic features (parallel high-spatial frequency pattern) and random features (hills, dales and scratches) subjected to crater-like, material removal phenomena (see figure 2).
To generate the simulated datasets: -Deterministic topography elements (e.g. the steplike feature of test case 1, pattern in test case 2) were procedurally generated by analytical equations. Random components were added as multi-scale overlays of spatially correlated random noise at different wavelengths and amplitudes.
-The modified topographies (i.e. after the modification events) were obtained by simulated material subtraction of further sets of deterministic and multi-scale, spatially interpolated random components.
-Final results (pre and post the modification event) were saved as triangle meshes.
-Simulated measurements (either on original or modified topography) were implemented by raster scanning, i.e. by intersecting the triangulated mesh with z rays located along the rows and columns of a x, y grid of 600×600 points and (1.5×1.5) μm point spacing, consistent with typical lateral resolutions and fields of view achievable with current areal topography measuring instruments (e.g. focus-variation, confocal, and coherence-scanning interferometry instruments [21]). To simulate error in local height measurement, each intersection point (between the z ray and the triangulated mesh-deterministic result) was disturbed by a random component along the z axis, sampled from a normal, random function N(0, σ z ) with σ z =0.2 μm, which represents a worst-case scenario in terms of the vertical precision of current instruments. To introduce localisation error, the position and orientation of the raster scanning grid was randomly displaced on the x, y plane, and rotated about the z axis. The maximum allowed displacement was 10% of the diagonal of the raster scanning grid, maximum rotation was 5°. Maximum lateral and angular displacements were found consistent with manual placement of a sample under the measuring instrument.
-Ten measurements were simulated on each test case (pre and post surfaces), each time with new simulated measurement error and new lateral and angular displacement. With two test cases, this led to a total of twenty registration problems of known solution.
-The results of measurement simulation were recorded as height maps, i.e. gridded sets of height values), the typical format of raw measurement data obtained by using areal topography measuring instruments. Measurements from topographies pre and post modification were referred to as 'pre' and 'post' datasets respectively.

Real test case
The real test case consisted of a block of AISI 316 L steel, with top surface obtained by face milling ('pre' status). The surface modification ('post' status) was imparted by mechanical interaction with a manual tool, in the form of a series of localised compressive deformations. The surface was measured in its 'pre' and 'post' status by means of an Alicona Infinite Focus G5, an areal topography measurement instrument based on focus-variation technology [21], using a 20×objective at 1×Zoom, NA 0.40, FOV (0.42×0.42) mm, lateral resolution 0.43 μm, acquiring regions of variable size, in stitching mode (i.e. by collating individual fields of view). The sample was physically moved away from the focus-variation microscope to be mechanically modified, and then placed back again under the microscope for the second measurement, without particular care for relocation, by only visually estimating consistency with the previous placement. In figure 3, height maps obtained by measurement are shown for regions cropped to an equal area of (2.54×2.80) mm. The real test case is particularly interesting because many portions of the machined surface are self-similar, making it difficult to visually identify corresponding landmarks in the pre and post topographies.

Proposed registration method
The proposed registration method consists of a series of steps, illustrated in the following paragraphs with the help of the simulated test case n. 1.

Topography data pre-processing
Firstly, the height maps are levelled by subtraction of the least-squares mean plane (the most common F-operator, as prescribed by ISO 25178-2 [9]).

Computation of local topography descriptors
A modified version of the SIFT algorithm [20], adapted to operate on height maps, is used to identify correspondences (referred to as keypoints) in the pre and post datasets. This is accomplished through the following steps. Height maps are converted into pseudo-intensity images by encoding height values as grayscale information. Then, a scale-space decomposition [20] is obtained as follows (figure 4): the intensity image is convolved with a Gaussian kernel multiple times, each time doubling the size of the kernel. This is roughly equivalent to halving the resolution of the image at each step (i.e. following a dyadic series), i.e. increasing the scale by one octave. Each level of the decomposition L x y , , s ( ) is thus defined as:  and L x y , , 2 s ( ) is used to obtain the levels of an approximated Laplacian pyramid as follows (figure 4, right): Consistent with the original SIFT method [20], initial candidates for the role of keypoint are identified as those local maxima in the Laplacian pyramid levels (i.e. high-curvature points) which would also be local maxima if considering the corresponding points in the above and below levels (figure 5).
The first set of tentative keypoints is then subjected to a series of post-processing operations, again following the SIFT method [20]: -The keypoint positions are refined to sub-pixel accuracy by locally fitting the region surrounding each keypoint by Taylor series expansion, and then by finding the new local maximum. The Taylor series expansion is obtained by considering the Laplacian pyramid as a function of three variables ) so that the expansion is: and the refined positions x corresponding to local maxima (i.e. x, y, σ triplets) can be found by setting the first-order derivatives to zero, i.e.: -Low contrast keypoints are eliminated by thresholding the derivatives of the curvature computed across scales.
-Keypoints corresponding to edges are eliminated by checking curvature across the two principal directions: an edge point will have larger curvature across the edge, and almost zero along the edge.
A topography descriptor is computed for each remaining keypoint, as illustrated in figure 6, again consistent with the SIFT method [20]. Firstly, the local gradient vector at the octave, and the position and level where the keypoint is located, is computed on the Gaussian scale-space decomposition. An example gradient vector is shown by the black arrow in figure 6(a). Magnitude and orientation of the gradient vector (m x y , ( ) and x y , q ( )) are obtained by means of finite differences as follows: The surroundings of the keypoint are then resampled over a grid of specified width and resolution, oriented as the gradient vector ( figure 6(a)). Within each cell of the resampling grid, local gradient vectors are computed (magnitude and angle), as shown in figure 6(b). The magnitudes are then weighed by a Gaussian function covering the entire resampling grid, so that gradient vectors located far from the keypoint are attenuated (see again 6.b). Finally, a three-dimensional binning of the gradient vectors is performed as shown in figure 6(c). Binning is on grid cells (a 4×4 configuration as shown in figure 6(c)) and on orientation values (eight bins is the default, as shown in figure 6(c)). A three-dimensional histogram is finally built with the average intensities (magnitudes) of the binned gradient vectors (in figure 6(c) this leads to a total of 4×4×8 histogram values). The histogram values are stored as the final SIFT descriptor value, together with the local orientation of the resampling grid computed at the keypoint.

Identification of initial correspondences
Keypoints computed from a pair of pre and post datasets, located at the same octave and pyramid level, and with similar SIFT descriptors, are selected as candidates for forming matched pairs. Their similarity is an indication of them being suitable to act as correspondences, as they likely belong to invariant portions of topography, except for measurement error. Therefore, in ideal conditions, the identification of a geometric transformation that maximises the degree of overlapping of the keypoints of each matched pair should correspond to the ideal registration solution. In reality, however, not all correspondences may be equally valid, and simple ranking by similarity of descriptors may not suffice. An example of such a condition is shown in figure 7, where an initial set of correspondences has been identified on a pair of pre and post datasets belonging to one of the datasets belonging to the simulated test case n. 1. Note that the keypoints are shown overlaid to the original datasets, whilst in reality they may belong to different octaves and levels of the associated pyramids.

Pruning of the correspondences and identification of the registration transform
In the conventional image processing scenarios where the SIFT is typically applied, i.e. the creation of panorama images via stitching of individual photographs [20], iterative methods based on random sample consensus (RANSAC) [22] are used to identify a homographic transformation which is consistent with the highest number of correspondences. The homographic transformation (non-rigid) is adopted because in stitching one wants to compensate for slight differences between images due to a change of scale and perspective, due to the different points of view from which the images were taken. Differently from stitching, in this work, the alignment of correspondences is used to achieve co-localisation for metrological purposes, thus only rigid transformations are allowed. For this purpose, a new variant to the original SIFT method is proposed, consisting of using RAN-SAC in combination with a rigid transformation, as  2. The 'model' is a rigid transformation leading to the maximum alignment of the keypoints belonging to the selected correspondences. The transformation is comprised of a rotation about the z axis, and a translation in the x, y plane, because the keypoints are defined in the x, y plane. The model is estimated starting from the selected correspondences, by solving a least-squares minimisation problem through the Procrustes superimposition method [14].
3. Once the model (rigid transformation) has been identified, all the available correspondences are tested against it. A correspondence fits well with the model if its two keypoints are aligned after the application of the transform (in which case the correspondence is said to agree with the model, and it is classified as 'inlier'). If the two keypoints are not aligned, the correspondence is considered as 'outlier', i.e. not in agreement with the model. To incorporate a level of acceptable error in the alignment, a threshold on Euclidean distance is set to determine whether a correspondence is an inlier or an outlier. Controlling parameter for this step: OTV (outlier threshold value).
4. The quality of the model can be assessed by computing the total percentage of inliers; the higher the percentage (larger consensus), the better the model. In figure 8, the results of pruning by RANSAC are shown on the same dataset that had originated 127 initial correspondences, as shown in figure 7. The final result amounts to 65 remaining correspondences (51.18% inliers). RANSAC parameters were: NCM: 3; OTV: 6 (pixels); MNI: 100.

Generation of registered pre and post datasets
The final registration transform identified by RAN-SAC pruning is comprised of a 2×2 rotation matrix R (or equivalently, a rotation angle about the z axis) and a 2×1 translation vector T (translation in the x, y plane). During pruning, the transformation is used to rotate and translate the keypoints of one of the datasets to align them to the keypoints of the other. This is a transformation in 3 dof, on the x, y plane, applied to two-dimensional point clouds. However, after pruning is complete, the final transformation must be applied to the entirety of a dataset in order to align it with the other. This implies that the 'moving' dataset (a height map) must be subjected to resampling by interpolation (if translation is not exactly a multiple of pixel width, and if there is any rotation). Thus, the moving dataset is resampled by linear interpolation at the same locations of the associated 'fixed' dataset, so that one-to-one comparison of height points between aligned height maps are possible. An example of registered pairs of height maps is shown in figure 9. The height maps are shown in false colours to better highlight the quality of alignment of the invariant

Results
In this section, the results of the application of the proposed registration method are illustrated for the simulated test cases n. 1 and n. 2, and for the real test case. A comprehensive, quantitative assessment of performance is only possible for the simulated test cases, as the extents of the surface modifications and the imposed initial misalignment are defined by design. However, all the test cases, including the real one, can still be visually inspected to appreciate many aspects of the registration result.
All the results have been obtained through an original Matlab [23] implementation of the registration method. All the code is native except for the computation of the SIFT descriptors, which relies on the VLFeat open source library [24].
The registration parameters for the entire procedure were the same for the simulated and real test cases, except for RANSAC OTV (see below). The full set of parameters is the following: The difference in OTV is likely imputable to a larger amount of self-similar regions in the simulated surfaces, better treated with a less strict pruning of correspondences.

Analysis of the localisation of correspondences
The first interesting element of investigation is to determine whether the correspondences identified by SIFT, and surviving the RANSAC pruning process, are located in regions of the topographies that have not undergone significant modification in the transition from the pre to the post condition (i.e. invariant regions, except for measurement error). This can be qualitatively assessed on any test case where modifications are clearly visible, but can be reliably assessed in a quantitative way only for simulated datasets where the regions affected by modification are clearly marked. In figure 10, a binary map indicating where the modifications are located for one of the datasets of the simulated test case n. 1 is shown along with the positions of the SIFT keypoints, which survived pruning. The dataset in figure 10 is the same as that previously illustrated in figures 7 and 8.
The analysis of the entire set of ten registration problems for the test case n. 1 and ten problems for the test case n. 2 indicated that, on average, 112 correspondences (pairs of keypoints) survived RANSAC pruning for case n. 1 and 685 for test case n. 2. Of the surviving correspondences, 99.92% (mean for test case 1) and 99.97% (mean for test case 2) fell within 'invariant' regions (i.e. both keypoints of the correspondence must lay on invariant regions for the correspondence to be counted as valid). Thus, the likelihood of encountering misplaced correspondences between modified and unmodified regions is expected to be low.

Registration performance and dependency on amount of initial misalignment
The results of the twenty simulated registration problems for the simulated test cases n. 1 and n. 2 were plotted against the simulated displacements to assess the overall performance in terms of residual alignment error, and to investigate whether registration performance would get worse or improve with the amount of initial misalignment. Separate plots were generated for translation (x and y displacement combined into a vector, vector length considered) and for rotation (angular displacement about the z axis). Lateral displacement values were expressed in pixel units (because all the simulated datasets have the same pixel width/spacing), whilst angular displacement values were expressed in degrees. Each plot was fitted to a   straight line by regression, to identify possible trends. In figure 11, 'initial' displacement values are those set by the simulation, whilst 'final' values are the residual errors remaining after application of the registration transform.
For the simulated test case n. 1, the analysis of figure 11(a) shows that the residual, translation displacement error after registration is on average one or two orders of magnitude smaller than the applied initial displacement. For a maximum simulated displacement of 77.2 pixels, the maximum residual lateral displacement after registration was 0.63 pixels (i.e. sub-pixel), with an average residual displacement of 0.43 pixels. For a maximum angular error generated by simulation (4.87°) ( figure 11(b)) the maximum final angular displacement error after registration was 0.07°, with an average residual displacement of 0.03°, i.e. two orders of magnitude smaller.
The results for the simulated test case n. 2 are shown in figure 12. The analysis of figure 12(a) shows that the residual, translation displacement error after registration was on average two or three orders of magnitude smaller than the applied initial displacement. For a maximum simulated displacement of 61.7 pixels, the maximum residual lateral displacement after registration was 0.03 pixels, with an average residual displacement of 0.02 pixels. For a maximum angular error generated by simulation (4.49°) ( figure 12(b)) the maximum final angular displacement error after registration was 0.005°, with an average residual angular displacement of 0.003°i .e. three orders of magnitude smaller.
The regression lines shown in figures 11 and 12, despite the less than ideal fitting results (as indicated by the R 2 coefficients reported in the figure captions), indicate little dependence of registration performance on the amount of imposed, initial displacement (almost horizontal regression lines and large residuals with no apparent trend). Interestingly, in all cases, registration performance improved with larger amounts of initial displacement.

Registration performance and dependency on the test case
As already hinted at in figures 11 and 12, the observed registration performance for the simulated test cases n. 1 and n. 2 was different, with registration for the test case n. 2 performing consistently better. This observation was confirmed by the box plots shown in figure 13, illustrating the distributions of the final lateral and angular errors (i.e. residual displacement after alignment). Not only the mean error was significantly smaller for the test case n. 2, but also the scatter of the performance values was smaller.

Registration performance for the real test case
In figure 14, a cross-section of the aligned topography datasets is shown. The surfaces are rendered in different colour to highlight the regions where one surface is above or below the other. Despite it being impossible to obtain an accurate, quantitative assessment of the quality of the registration result (as the real initial misalignment is unknown), a simple visual inspection of the overlapping between the machining marks of the two datasets (see figure 14(b)) reveals a very good result. Moreover, the inspection of the local surface heights not only reveals the depressed regions corresponding to the mechanical modification processes, which were immediately visible even before the   pre-post comparison, but also highlights the presence of raised regions surrounding each depression, consistent with the the surface modification process.
A successful relocation allows for individual surface topography modifications to be quantitatively assessed. For example, in figure 15 it is shown how the void volume of an individual depression can be quantified.

Discussion
The results presented in this work show that the proposed registration method based on automated identification of correspondences is successful at relocating topography data measured pre and post a modification event. The analysis of the results on the simulated datasets show that, in most cases, the residual displacement errors (i.e. the residual displacement after registration) are at the sub-pixel level. Assuming measurements by optical instrumentation, in most cases the observed errors would be below the optical resolution limits [21]. The results also appear to be qualitatively confirmed for the real test case, despite the lack of a proper reference for a comprehensive, quantitative assessment. However, a number of questions remain to be addressed.
The biggest, looming issue to address, pertains to whether the topographies generated by simulation and the real test case investigated, are sufficiently representative for the wide array of diverse scenarios which may be observed in real industrial practice. How would the performance change with different topographies? What about smaller or larger modifications? The simulated test cases n. 1 and n. 2 feature topographies with different amounts of determinism and randomness, modifications that have different spatial distribution and volumetric relevance, presence of additional topography features which may weigh-in more or less negatively in the determination of reliable registration landmarks. The results for the two simulated test cases show important differences, in particular concerning the performance related to compensation of linear displacement. Despite the fact that one real test case was presented, showing that the method can be successfully applied (at least, up to what is visually appreciable), it would be prudent to test the method on a larger number of measured surfaces, and this work will form the basis of future publications.
Testing on real surfaces is far from being straightforward though. Major challenges remain unsolved on how to actually evaluate the 'goodness' of a registration result. Real measurements would lack associated knowledge of the ideal registration transform, which on the contrary was fundamental in this work to assess performance on the simulated test cases. If one can only rely on the information contained within the measured datasets, and no external common reference for co-localisation can be adopted, then the same metrics used to solve the minimisation problems underlying registration could be adopted as approximations. In other words, one may operate under the assumption that the real problem of minimising the distance between a dataset and its ideal position (relocation problem) can be approximated by either minimising the distance between visible landmarks that have been recognised as occupying the same position, or by minimising the global distance between the two datasets. The latter approach had been discussed in section 1 under the name of 'global registration', and is demonstrated to fail in the presence of significant modifications between the pre and post datasets. Regardless, even operating by aligning landmarks is marred by approximation, thus further work is needed to set up experimental conditions where an actual, external co-localisation reference can be set up to validate the proposed registration methods on real measurements.
Another issue is whether the current method could be further improved, either by operating on its many control parameters (e.g. the SIFT parameters, or the RANSAC pruning parameters), or by pairing it with a further fine registration method, such as one derived by ICP, as already investigated on surface topography [15,16]). Experimental investigations are in progress, but again the problem is that the ICP method is intrinsically global, and must be turned into a selective technique if one wants to avoid overfitting of significantly modified regions. As the current results already lead to sub-pixel accurate registration results, it may very well be that further refinement is not needed, as it would be difficult to evaluate actual improvement when operating below the measurement resolution limits.
Finally, the problem of how to define and compute uncertainty associated to relocation is of fundamental importance, and has not been addressed yet. Both because it is still unclear how uncertainty should be associated to areal topography datasets such as those serving as inputs to the method described here (see [24][25][26] for a discussion of the challenges of building an uncertainty budget for a measured areal topography dataset), and because it is still unclear how such information should propagate through the relocation method, finally affecting the registration result. Clearly, knowledge of relocation-associated uncertainty is fundamental to later estimate the uncertainty associated with any topographic comparison, and consequently, the traceability of the whole data processing and analysis chain.

Conclusions
In this work, a method was proposed for the geometric registration (i.e. co-localisation within the same coordinate system) of surface datasets obtained by areal topography measuring instruments. The registration method is useful in many surface measurement scenarios where the interest is the quantitative assessment of surface topography evolution as a consequence of modification events, e.g. due to functional life or mechanical testing.
As opposed to global alignment solutions, the method proposed in this work is based on the algorithmic identification of invariant regions which can act as fiducial landmarks to be superimposed in a selective alignment process. Correspondences between landmarks are algorithmically recognised by computing the similarity of candidate regions encoded through a dedicated shape descriptor adapted from the scale invariant feature transform (SIFT); a popular descriptor originally devised for image processing. The rigid transformation (rotation and translation) leading to alignment is then computed by solving an absolute orientation/Procrustes superimposition problem.
The performance of the proposed procedure was quantitatively assessed through application to simulated test cases of known solution. The performance of the method was also qualitatively assessed on a real test case involving a machined surface subjected to mechanical modification and measured by focusvariation microscopy. For the simulated test cases, the registration results indicate residual errors at the sub-pixel scale, usually below the resolution limits of optical areal topography measuring systems. Visual inspection of the real test case also indicates good relocation results, with a visually accurate alignment of visible surface features (machining marks).
Future work will involve additional refinements to the registration method based on understanding performance when applied to a wider array of surface types, the development of experimental solutions to validate the method on real measurement results and the development of methods to estimate the associated measurement uncertainty.