Automatic drift elimination in probe microscope images based on techniques of counter-scanning and topography feature recognition

An experimentally proved method for the automatic correction of drift-distorted surface topography obtained with a scanning probe microscope (SPM) is suggested. Drift-produced distortions are described by linear transformations valid for the case of rather slow changing of the microscope drift velocity. One or two pairs of counter-scanned images (CSIs) of surface topography are used as initial data. To correct distortions, it is required to recognize the same surface feature within each CSI and to determine the feature lateral coordinates. Solving a system of linear equations, the linear transformation coefficients suitable for CSI correction in the lateral and the vertical planes are found. After matching the corrected CSIs, topography averaging is carried out in the overlap area. Recommendations are given that help both estimate the drift correction error and obtain the corrected images where the error does not exceed some preliminarily specified value. Two nonlinear correction approaches based on the linear one are suggested that provide a greater precision of drift elimination. Depending on the scale and the measurement conditions as well as the correction approach applied, the maximal error may be decreased from 8–25% to 0.6–3%, typical mean error within the area of corrected image is 0.07–1.5%. The method developed permits us to recover drift-distorted topography segments/apertures obtained by using feature-oriented scanning. The suggested method may be applied to any instrument of the SPM family.


Introduction
The precision of surface topography measurement, physical dimensions of the elements varying in the range from several angstroms to several tens of nanometres, is mostly defined by the value of the drift of the scanning probe microscope (SPM). As a rule, the instrument drift includes two main components: one is caused by the thermal deformation of mechanical units of the device, and the other results from the creep of piezomanipulators applied [1]. Elimination of the drift can either be done by means of active compensation [2] during the measurement or image correction [3,4] after the measurement or a combination of both methods.
Image-correcting methods, in particular the method [5] proposed in the given paper, as compared to compensating methods, have the advantage that no modifications are required to be done to the microscope in order to eliminate distortions.
Moreover, unlike the passive approaches, active drift compensation introduces additional noises, which prevents it from being used for measurements at the highest microscope resolution.
The drift correction method developed is based on a simple linear system of equations composed for topography features found on each of the counter-scanned images of a pair. A more complex version of the method, the linear drift correction being carried out by two simultaneously obtained pairs of counter-scanned images, permits us to reach notably better results. Finally, transition to nonlinear methods built on the base of the found linear solutions ensures the greatest accuracy of the drift correction. Application of the nonlinear methods is especially effective in the case of a strong nonlinear component of distortion produced by creep of microscope piezomanipulators.
The principal concern in the work is devoted to using the proposed methods in feature-oriented scanning (FOS) [6], namely for drift correction in segments and in apertures. FOS is a method intended for a high-accuracy topography measurement using surface features as reference points of the microscope probe. With this method, during successive passings from one surface feature to another one located nearby, what are being measured are the relative distance between these features and the topography of their neighbourhoods-segments (apertures). That permits us to scan the required area on the surface by parts and then reconstruct the whole image by the obtained fragments.
The suggested drift correction methods have been tested so as to check the operation on different types of probe instruments, different types of surfaces, in various scales, for several drift velocities and with different proportions of thermal and creep drift components. Comparing the obtained results with the previously achieved ones [3] shows that the distance between neighbouring carbon atoms on a graphite surface may now be measured with the error lying in the interval 0.01-0.33% against 5%. Increase in drift correction accuracy is achieved due to the use of analytical expressions allowing for all components of linear raster distortion, transition to two pairs of counter-scanned images, application of nonlinear correction methods as well as methods of automatic topography feature recognition.

Linear drift correction by one pair of counter-scanned images
2.1.1. Building system of equations. Analysis of distortions caused by drift of a microscope probe relative to a sample surface proves that lateral-plane drift would lead to image stretching/contraction along x and y raster axes as well as to picture skewness due to a shift of the image lines/columns relative to each other. The same is observed in the vertical plane with respect to the topography height. Here, the height differences are imaged incorrectly and an additional spurious surface tilt appears. Assuming the drift velocity to change slowly enough while scanning a small-sized image [5,6], the described distortions may be represented in the form of linear transformations as follows: x(x, y) = x + (K x − 1){x + [(k + 1)m x + 1]y}, y(x, y) = y + (K y − 1){x + [(k + 1)m x + 1]y}, wherex,ȳ,z are coordinates of points in the corrected image; x, y, z are coordinates of points in the drift-distorted image; K x , K y , K z are linear transformation coefficients (LTCs); k is the ratio of probe velocity v x in the forward scan line versus probe velocity in the backward scan line; m x is the number (but one) of points in a line of the distorted image, which defines the range of the variable x = 0, . . . , m x . Besides the condition of an invariable drift velocity, equations (1) are written with the assumption that single steps of the microscope in the lateral plane are equal x = y and so are the movement velocities v x = v y . The K x − 1, K y − 1, K z − 1 factors of transformations (1) describe a displacement caused by the drift along x, y, z, respectively, while moving the probe by one step along x or y. The parameter k is to consider the probe displacement along x, y, or z accumulated during the retrace sweep.
Thus, the (K x − 1)x, (K y − 1)x, (K z − 1)x terms allow for drift-induced probe displacement along x, y, z, respectively, that is occurring while moving the probe along the current line. The (K x − 1)(k + 1)m x y, (K y − 1)(k + 1)m x y, (K z − 1) × (k + 1)m x y terms represent probe displacement along x, y, z, respectively, that took place when moving the probe by the previous lines. The (K x − 1)y, (K y − 1)y, (K z − 1)y terms are responsible for probe displacement along x, y, z, respectively, that occurred when moving the probe between the raster lines.
To find the unknown LTCs K x , K y , and K z , one may scan a surface by the trajectory schematically shown in figure 1(a). As a result, a pair of images is obtained (see figure 2(a), position 1) where the lines of the images are scanned in opposite directions and line-to-line movements occur in opposite directions too. Scanning velocities are set equal for these images. Such images will be referred to as counter-scanned images (CSIs). Typical of CSIs is the existence of a point common for both images (see figure 1). That point will be referred to as a coincidence point (CP). The CP is the end point of the raster trajectory of the first direct image and it is the start point of the raster trajectory of the second image counter to the first one. Moving away from the CP, because of drift the images become more different from each other, namely, the distinctions in positions of the same features become more noticeable and the features themselves undergo mutually-opposite transformations (stretchings/contractions and skewnesses).
Provided the same surface feature is present in each image of the obtained pair (see figure 2(a), position 2), then the following system of equations may be composed by using the lateral coordinates (x 1 , y 1 ), (x 2 , y 2 ) of the featurē where digits 1 and 2 in the designations show that the quantity is attributed to the first (direct) or to the second (counter) image, respectively; m y is the number (but one) of points in a column of the distorted image, which defines the range of the variable y = 0, . . . , m y . The items m x and m y in equations (2) provide a transformation of the first image coordinate system to the second image coordinate system (the origin of the second image coordinate system is CP).
Generally, to correct the drift-induced distortions, it is sufficient to only reveal one feature on the CSIs and to determine its lateral coordinates. Since real SPM images have a finite resolution, are noisy, and contain corrupted regions, for the correction parameters to be determined more precisely, it is desirable to use all the features available on the surface excluding, may be, those located along edges of the images.
The fact is that the edges of SPM images are usually distorted nonlinearly by creep [7], by hysteresis [8], as well as by motion dynamics of the piezomanipulators. As a rule, provided the edge distortions are quite distinct, the margins of the scanned image are just discarded after correction.
Those features more distant from CP provide a better precision while determining LTCs. On the other hand, the probability of change in drift velocity is also increasing while moving away from CP. Therefore, selecting a particular feature is each time a subject of compromise: on small-sized scans, all the features are usable, while on large-sized scans, some of the features should be sometimes 'sacrificed' (see section 3.1.1).
Thus, representing the totality of features by its 'gravity centre' with coordinates (x 1 , y 1 ) and (x 2 , y 2 ) in the respective CSI and taking into consideration transformations (1), equations (2) may be rewritten as In the CSIs, the relation between coefficients K1 and K2 is very simple: If one of the coefficients in a pair were to stretch the image, then the other one would contract it and vice versa. If either coefficient translates the image with no distortion, i.e., is equal to 1 (drift is absent), then the other one will also translate the image with no distortion, i.e., will also be equal to 1.
By substituting the coefficients K2 x and K2 y from (4) in equation (3), the sought LTCs K1 x , K1 y for the first image are found (position 3) Then, the LTCs K2 x and K2 y for the second image may be determined through the relationships (4). After that, using the obtained coefficients, the CSIs 1 and 2 are corrected in the lateral plane by means of transformations (1) (position 4). Strictly speaking, one should distinguish the coefficients K z for positive K + z and for negative (K − z ) topography differences z. Therefore, the last expression from (4), which shows the general connection between CSI coefficients K z , should be represented as Whence the following relations may also be derived: In accordance with the definition of the K z coefficient for height differences at point (x,ȳ) located inside the overlap area of the laterally-corrected images 1 and 2 (position 5), the following system of equations may be written: where z 1 , z 2 are distorted height differences in the laterally-corrected images 1 and 2, respectively; z is the true height difference of the topography. By substituting coefficient K2 − z from (6) in the obtained linear system, the solutions could be found as follows: The obtained solution for z shows that the true topography height difference is defined by a half-sum of the distorted height differences. It should be noted that for various height differences z various coefficients K z will be obtained, since the vertical drift induces the same displacement during time t of lateral step execution (coefficients K z used in formula (1) correspond to difference z = 1).
It is convenient therefore to manipulate with the value of vertical displacement (K z − 1) z, which is the same in all points of the overlap area, rather than with the height difference z and corresponding coefficient K z . Using formulae (8), it is easy to determine that the vertical displacement is equal to a half-difference of the height differences z 1 and z 2 . The vertical displacement is calculated at each point of the overlap area and is then averaged.
Thus, the distortions in the vertical plane may be corrected applying formula (1), where the obtained mean vertical displacement should be used instead of K z − 1. However, this correction method is not accurate since the drift is unable to noticeably distort the height difference z during time t of lateral step execution.
For points belonging to the overlap area of laterallycorrected CSIs (position 5), the following equation may be composed:z In the lateral plane, transformations inverse to (1) are the following: Using identityz(x, y) ≡z(x,ȳ), where coordinatesx,ȳ of the overlap area relate to coordinates x, y of the rectangular raster of the corresponding CSI in accordance with transformations (10), and expression forz from system (1), equation (9) may be represented in the expanded notation as Shifts m x , m y of coordinates x 1 , y 1 , respectively, provide a return from the coordinate system of the counter image to the coordinate system of the direct image. The plus sign at factor (K2 − z − 1) is used because the sign of difference z in the counter image is opposite to the sign of difference z in the direct image (here, | z| = 1). The item K1 + z − 1 {m x + [(k + 1)m x + 1]m y } in the right hand side of equation (11) allows for the vertical displacement of the origin of coordinates of the counter image that occurred during direct image scanning.
By using relation (6) between the coefficients K1 z and K2 z , it is easy to express the coefficient K1 z from equation (11) as The coefficient K1 z is calculated at all (x,ȳ) points of the overlap area and then is averaged out (position 6). As the coefficient K z is known, each image may be corrected in the vertical plane according to the expression forz of transformations (1) (position 7). Finally, the obtained images are matched in CP and then the topography is averaged within the overlap area (position 8). Thus, at the output we have got a corrected surface image free of distortions caused by x, y and z drifts, noise level reduced.

Feature recognition.
A search for the features in CSIs and determination of their lateral coordinates may be carried out manually. However, the use of a recognition procedure [6,9] would enable completely automatic correction of driftdistorted SPM images. With that procedure, hill-or pit-like topography elements are taken as the features. Since the surface features are defined in the most general form, in most cases in practice it turns out quite feasible to find a suitable feature in the image.
The recognition procedure supports operating with hills only, with pits only, or with both feature types at once. In respect of increasing precision of LTC determination, advantage of the latter option should be taken as it provides the maximal number of features to be engaged. Coordinates of the 'gravity centre' of a feature serve as coordinates of the feature position. This is quite acceptable owing to the linear character of the distortions. Before recognition, it is recommended that the mean surface tilt be removed and the picture be smoothed. Note that the above manipulations are only executed upon the image duplicates, the originals being subjected just to drift correction.
Since the method developed implies recognition of a scanned image, the topography features should be understood in the broad sense. Physically, they can refer not only to a topography but also to magnetization domains, to places of localized electric charge, and so forth depending on the type of probe microscope used.

Iterative search for feature pairs.
Because of driftinduced distortions, the positions of the same features in CSIs do not match, the degree of mismatch would increase moving away from CP. Neither is the number of features in CSIs equal (see section 3). Moreover, some of the features recognized in one of the images may be unrecognized in the other image because of scanning faults. Thus, after feature recognition in CSIs, an iterative process should be carried out in order to establish the fact that a feature with coordinates (x 1 , y 1 ) in image 1 and a feature with coordinates (x 2 , y 2 ) in image 2 are both the same feature with coordinates (x,ȳ) in the corrected image.
In order to do so, from the image 1 feature list, a feature is chosen such as to be first met while moving from CP along the counter scan trajectory. Then, in the image 2 feature list, a feature is searched for such as to be the nearest to the feature chosen in image 1 and to lie within its certain circular neighbourhood. While searching, only those features are considered that have the same type (hill/pit) as the feature chosen in image 1. If no suitable feature is found in the image 2 feature list then the feature chosen from the image 1 feature list is removed from that list and a new iteration is started.
The coordinates of the features chosen from the lists and intended to detect the next pair are corrected preliminarily by means of transformations (1)  Once a suitable feature has been detected in image 2, the coordinates may be determined of the gravity centres of the feature sets selected by this moment (if a substantial change in the drift velocity has occurred then not all selected features should be used but only part of them obtained for the last several iterative cycles). Then, by formulae (5) (4). The iterative process described is repeated until all features from the feature list of image 1 are analysed.
Since the method suggested implies matching of the gravity centres of the feature collections detected in the images 1 and 2 during LTC determination, the sum of square deviations of the corrected feature positions will be the least possible. The advantage of the proposed algorithm becomes obvious when applied to highly ordered surfaces consisting of identical elements. Provided the number of defects is small, there is no other way of distinguishing one surface feature from another except the method suggested.
To reveal the feature pairs successfully, the scanning area should be located so that several features lie in the vicinity of a CP. When all the features of a measured surface area are concentrated at the bottom part of the CSI and, in addition, are very little distinguished from one another by size, by form and by mutual position, in order to detect the feature pairs correctly a top-down → bottom-up counter-scanning should be carried out rather than the bottom-up → top-down scanning depicted in figure 1.

Linear drift correction by two pairs of counter-scanned images
The accuracy of correction of the surface topography may be improved by applying a drift elimination mode based on two CSI pairs. To implement that mode, the movement velocity at the retrace is set equal to the movement velocity at the forward trace (k = 1) and scanning is carried out by the trajectory shown in figure 1(b). After the LTCs have been found in the way described above, the images of each pair are corrected (see figure 2(b), positions 2-8).
If the matching error for the first image pair is small, i.e., the microscope drifts with practically constant velocity, then determining the second pair coefficients may be simplified since the following equalities are valid: The inverse statement is also correct. These equalities may serve to additionally check the drift correction scheme suggested. Thus, if coefficient ). If K1 y < 1 then K3 y < 1 as well and if K1 y > 1 then K3 y > 1 as well (similar relationships occur for K z ).
The obtained pair of corrected images is matched again and the topography is averaged out within the overlap area of these images (position 9), which results in an additional damping of the noise level. As each CSI pair has a CP of its own (see figure 1(b)), the matching of the corrected images is carried out by matching the gravity centres of the features (here, for noisy images the least mean squares criterion would work again). Before calculating positions of the gravity centres, the features should be excluded from the feature lists that are contained in one CSI pair and are absent in the other one.

Determination of x, y, z components of drift velocity
If LTCs K x , K y , K z are known, the mean velocityv for corresponding drift components is determined by formulaē where z is the minimal length change of the Z manipulator. Having the drift componentsv x ,v y , it is possible to find the modulus of the drift vector in the lateral plane |v xy | = √v 2 x +v 2 y . The vertical drift velocityv z may also be determined as follows. Let us write down the last term of the last equation of transformations (1) for the point with coordinates (m x , m y ). The obtained value will correspond to the topography ascent/descent caused by the vertical drift component that occurred during the scanning time t of the entire image (direct or counter). By calculating the obtained value versus scanning time ratio, it will give the mean vertical drift velocity as Comparing formulae (13) and (14), the following identity may be written: where the values on the left and on the right of the identity sign are total numbers of steps in a raster. In general terms, the third method of determination of the vertical drift velocity consists in measuring the height difference between the last and the first points of the direct scan, measuring the height difference between the last and the first points of the counter scan; then calculating the half-sum of the obtained differences followed by dividing the found value by the value of scanning time t for one image (halfsum calculation permits us to exclude the own mean tilt of the surface).
Since with CSIs the last point of the counter scan does not coincide with the first point of the direct scan because of influence of the lateral drift components, to exclude the effect of the topography upon the final result, the topography in each CSI should be replaced with the mean plane using the least mean squares method. Then, the pointed height differences may be determined by those planes.
When drawing the mean planes, only the points from the overlap area should be used since the topography is represented by the same features there. It is erroneous to calculate the height difference immediately by the topography of the overlap area since the drift-induced distortions within this area have being acting for different times at the direct and the counter regions.
The advantage of the described method is that in order to estimate the vertical drift velocityv z the coefficient K z is not required. Moreover, that coefficient itself may be found by the value of this velocity according to formula (13) (or formula (14)).

Nonlinear drift correction
For surface scans which are characterized by a large number of features distributed quite regularly over the image area, the following nonlinear correction method could be suggested. First, applying the linear approach described above, the pairs of CSI features are determined; then the local LTCs K x , K y are calculated by each feature pair. As a result, a distribution is obtained such that each feature with coordinates x, y has corresponding lateral LTCs K x (x, y), K y (x, y).
Using transformations (1), local displacements D x (x, y) =x(x, y) − x, D y (x, y) =ȳ(x, y) − y in the lateral plane are found for every image feature with the real coordinates x, y. The displacements corresponding to the integervalued coordinates of the points of the distorted image are determined by regression surfaces drawn through the obtained displacements. Finally, the image corrected in the lateral plane is acquired by applying the appropriate local displacements to the points of the distorted SPM image.
To correct drift in the vertical plane, the CSIs corrected in the lateral plane should be superimposed, one on another, by matching the gravity centres of the feature sets being used. By calculating the local coefficient K z and then the difference between the corrected and the distorted topography heights at each point (x,ȳ) of the CSI overlap area, the distribution of local displacements across the image field in the vertical plane D z (x,ȳ) =z(x,ȳ) − z(x,ȳ) may be found. Building a regression surface through the found local displacements, the working distribution is determined, which will help perform nonlinear image correction in the vertical plane.
Constructing regression surfaces also permits us to reduce the influence of the error of determining feature lateral position and topography height on the results of nonlinear correction. The regression surface order is chosen by the residual mismatch of the feature positions (topography heights), so that the mismatch is minimal. In the case of large nonlinear distortions, the regression surfaces may also be used during iterative search for feature pairs.
Another nonlinear correction scheme may be suggested. First, some square neighbourhood (a segment) is cut out around a feature in every corrected image. Then the cut-out topography fragments are put in position, which is the mean of the corrected positions of this feature in the corresponding CSI. Correction of the images and the feature positions may be carried out by either linear or nonlinear methods as described above. Finally, within the segment overlap areas, the topography is averaged. Thus, the corrected topography is getting nearer to the actual one not only because of averaging in the vertical plane but also because of averaging of the feature positions.
The basis of the described nonlinear correction methods is that the true feature position lies somewhere in the segment between the corrected feature positions, most probably gravitating to the middle.

Correction of drift-distorted topography segments in a feature-oriented scanning method
The described drift correction method yields the maximal effect when used with the FOS approach [6], since in that case it would enable drift correction in images of an arbitrary large size. The point is that, beginning with a certain scan size, the main assumption of invariability of the drift velocity during scanning time will necessarily cease to work (see section 3). Although the same assumption must be satisfied for scans obtained by the FOS method, the contradiction in that case is eliminated because a large area is scanned by parts, i.e., by small segments (square neighbourhoods of the surface features) and all movements occur within short distances from one feature to another located nearby.
To make sure of the above, one should compare 1.5 min approximately scanning time of (35 × 35)Å 2 atomic graphite surface (see section 3.1) with 300 ms approximately scanningrecognition time of (4.5 × 4.5)Å 2 segments of the 'Next' and the 'Current' carbon atoms in one skipping cycle [6]. Skipping is a basic measurement operation in FOS intended for accurate determination of relative coordinates of neighbouring features and acquisition of topography segments.
Passing from the atomic scale to surfaces with typical feature dimensions and distances between the features of tens and hundreds of nanometres, the visible evidence of the drift thermocomponent in the image would subside, though the nonlinear creep component of the drift, in contrast, would become more apparent (see section 4). Nevertheless, the negative creep effect may be substantially reduced provided that the topography is measured by parts using small segments and mutually opposite probe movements forming the entire hierarchy of counter movements are applied throughout the apertures (an aperture is an auxiliary scan of the current feature surroundings containing several neighbouring features), in the segments, between the neighbour features (the skipping) as well as in feature lines [6].
Moreover, in the case that a noticeable change in drift velocity occurred (the drift is being monitored continuously during FOS) the measurement process is automatically suspended, the corrupted local data are discarded, and the microscope waits for the drift velocity to become constant, at that executing periodical probe attachments to the current surface feature or inserting idle skipping cycles. Once the drift velocity has become stable, the work is resumed and the interrupted local measurement is executed over again. Thus, because of applying the pointed set of the methods, the total drift turns out to be a slowly changing process again and, therefore, it can be linearized as well.
It should be noted how simple it is to detect the same feature in CSIs with the FOS approach. The fact is that a segment, as a rule, contains one feature only. When a segment includes several features (usually two or three), the features located closer than the others to the CSI centres will correspond to the same current feature in the corrected segment since the main sign of the current feature in the segment is its proximity to the centre of the square raster [6].
Recognition of two/four CSIs may be carried out simultaneously. Recognition of hill-like features and pitlike features in each image may also be processed in parallel [6].
Moreover, the process of direct image recognition and the process of counter image scanning may be implemented concurrently. Thus, the calculation throughput may be substantially increased by applying a two/four processor computer. Maximal advantage can be taken of a multiprocessor computer when the counter-scanning method is used within the FOS approach, where the recognition should be carried out in real-time.
It should be noted in conclusion that, applying the described drift correction method directly, the number of image averagings is restricted to 4 at most, whereas there are no restrictions at all on the number of image (segment) averagings if the proposed method is used within the FOS approach.

Experimental results
In order to verify the operation of the suggested drift correction method, a counter-scanning of an atomic surface of highly oriented pyrolytic graphite (HOPG) and a porous alumina surface was carried out. The graphite surface was measured with a scanning tunnelling microscope (STM) and the porous alumina surface was measured with an atomic-force microscope (AFM). In both cases an SPM Solver TM P4 (NT-MDT Co.) was employed, the measurements were performed in ambient conditions, the sample was moved relative to a fixed probe.
In order to provide a smoother transition of the piezomanipulators from the quiescent state to the scanning state, several tens of 'training' probe passages along the first line of the raster were carried out just before the scannings started. In this way, it allows us to substantially decrease the creep-caused distortions at the beginning of the scan. As well, during the training, the actual scanning speed v x was determined. Feature recognition in CSIs was executed in the course of virtual FOS [6].

Correction of drift-distorted scan of the atomic surface of pyrolytic graphite
In figure  The results of recognition of the obtained graphite CSIs are shown in figure 4. The atoms of carbon (hills) and the interstices (pits) are used as features of the current topography. In table 1, for each CSI presented are the number of features found, the value of mean lattice constant a, and its relative measurement error determined during the recognition (the standard value for the HOPG lattice constant is equal to 2.464Å). figure 5, marked with a '+' sign are those features detected in both CSIs during the iterative search for feature pairs. There were 424 common features detected in the first CSI pair and 430 in the second   one. Searching among the detected pairs for such a pair whose features are located the widest apart, for example along the x coordinate, it is possible to determine the absolute drift-induced measurement error in the lateral plane x max . Maximal relative error δ max for one image is calculated by the formula

Linear drift correction. In
The relative error along the y coordinate is calculated likewise. The largest of the two found errors is then selected. The error for the first CSI pair is 7.8%, for the second 8.0%. The lateral LTCs obtained after applying the suggested method are given in table 2. The greater the distance from a feature to CP, the less sensitive are the found LTCs to the error of determining this feature position. Therefore, it is possible to find the number of significant digits by  The results of drift correction in the first and in the second CSI pairs are presented in figures 6(a) and (c). Some decrease in the peak-to-peak height difference (see the vertical scale) at the corrected images as compared to the original ones points out that the topography has become smoother after correction (due to residual mismatch and noise level reduction and because the sizes of the corrected surface area are less than the original ones). The topography image presented in figure 7(a) is a result of matching the corrected CSIs shown in figures 6(a) and (c) followed by topography averaging within the overlap area. Here, one should take note of an even lower noise level despite the fact that there was no additional 'artificial' topography smoothing. Matching precision of the CSI pairs may be estimated by sight at figure 7(b). In digits, the mean residual mismatch of the corrected feature positions is 0.05Å. Figure 8 shows the dependence of the average of mutual mismatches of the corrected positions of the feature pairs upon pair serial number n, which was obtained during the iterative search for feature pairs. Remember that every feature receives its number according to its distance away from the   If the drift velocity stayed unchanged while scanning, the error would approach zero as new features were involved. However, it is seen in the plot that after a minimum is reached (C point) the error starts to grow (section C-D). Thus, the observed increase in error points out the fact that the drift velocity has changed during the counter-scanning.
As the error threshold is specified within the CD section, those dependences will help to determine the features whose residual mismatch does not exceed the error. Then, by using  the features found, it becomes possible to calculate the LTCs that provide drift correction with the specified error. The resulting image is obtained after discarding the part of the corrected image located beneath the position of the last of these features.

Nonlinear drift correction.
In figure 9 are shown the third order regression surfaces drawn through the local D1 x (x, y), D1 y (x, y) displacements of the features of the first CSI. The results of the nonlinear correction in the lateral plane for each CSI pair are given in figure 10. In figure 11 is shown the image of figure 10 having undergone matching followed by averaging in the vertical plane. The residual matching errors typical of the nonlinear correction method suggested are given in table 3. The obtained results would point out the fact that the nonlinear approach, as compared to the linear one, provides higher correction precision of topography due to a more accurate matching of the features being used (compare figure 10 with figure 6). Pictures corrected with the linear method are typically blurred after matching, mostly in the upper and in the lower parts. In contrast, the image reconstructed from separate  fragments as shown in figure 12 has no blur at all (compare with figure 7(a)). However, it is only the FOS approach to realize the idea of element-wise correction in full measure. The segmented structure of the image in figure 12 is only noticeable by the jagged edges; no other segment-specific artefacts could be found in that picture. It should be noted that the peak-to-peak height difference in figure 12 (see vertical scale) is greater than the peak-to-peak height difference in figure 7(a), which points out a more adequate topography restoration in the vertical plane.

Lateral calibration of scanner. Mean lattice constant
a calculated by data presented in figure 12 is equal to 2.113Å. The obtained value corresponds to a relative systematic measurement error of 14%. The above error results from the last incorrect calibration that has been done to the microscope XY scanner disregarding the drift distortions. Applying the automatic calibration method described in [9], the correctional coefficientsK x = 1.1896,K y = 1.1429 and the obliquity angle α = 0.4 • were determined (the latter is to allow for a certain nonorthogonality between X and Y piezomanipulators). This scanner calibration was carried out by the equilateral triangles found in figure 12. The only triangles having sides five times the lattice constant of graphite were chosen from the totality of triangles. 184 suitable triangles were revealed among the carbon atoms and 175 among the interstices.
Numerically, the lateral calibration coefficients found turned out to be very close to each other:K x x = 0.306Å, K y y = 0.307Å. This fact indicates a practical identity of the scanning channels of the microscope being used (a thorough identity of the channels is provided by the manufacturer during a special adjustment process). It should be noted that the identity of the lateral calibration coefficients points out the absence of disproportion in the obtained images; the identity may also be regarded as a confirmation of the validity of the performed drift correction.
After correcting scale and obliquity in the image shown in figure 12, the graphite atomic surface will look as presented in figure 13(a). The mean lattice constant a in the obtained image is equal to 2.4638Å, which corresponds to the relative measurement error of 0.01%, approx. A tiny difference of 0.001Å between calibration coefficients is caused by both errors of the scanning channel adjustment and errors of the drift correction method. Thus, another estimate of measurement error for a may be found, about 0.33%. The obtained error values are in accordance with the errors presented in table 3.
In figure 13(b), a stylized image of the graphite surface is demonstrated where carbon atoms are shown symbolically as spheres. Using stylized images allows us to represent schematically, as a generalized view, the structure of a surface consisting of the same elements, as well as to visually reveal surface defects that are indistinguishable in a real image. For instance, near the upper edge of figure 13(b), a descent of the atomic relief is well noticeable but it is absolutely unnoticeable in figure 13(a). In this case, an edge artefact occurs, which is typical of the Fourier smoothing method applied (see figure 4). To get rid of the flaw, the edges of the image should be cut off. The stylized images also enable easy observation of an offset of the interstitial positions (not shown in figure 13(b)) relative to the atom positions, which may point to a crystal structure peculiarity and/or a defect but most often is a consequence of probe tip asymmetry.
It should be noted that determination of the calibration coefficients and the obliquity angle also allows confirmation of the validity of the proposed drift correction approach. The above parameters should vary insignificantly near the scanner's origin of coordinates with various drift values and directions, different image sizes and scanning velocities. The pointed out criterion has been practically confirmed by executing a series of test measurements.

Calculation of drift velocity.
Drift velocities in the lateral plane are numerically determined by formulae (13). In order to take account of the calibration error discovered, the calculated velocitiesv x andv y were multiplied by the coefficientsK x andK y , respectively. The obtained values are presented in table 4.  In the case that the atomic graphite surface is being scanned in the constant Z mode, the drift velocity in the vertical plane cannot be determined. The reason is that the tunnel current slow changes induced by the vertical drift component are compensated by the servosystem in this mode. As a result, the current image is practically undistorted by the vertical drift.

Correction of drift-distorted scan of porous alumina surface
In figure 14, two CSI pairs of a quasiordered porous surface of alumina are shown. The pores in alumina are made by anodizing an aluminium foil in an aqueous solution of oxalic acid [10]. To better discern the morphology details, mean surface tilts were preliminarily subtracted from the images presented hereafter. On a defect-free surface, each pore is surrounded by six neighbouring pores, which form a regular hexagon. Around the orifice of every pore, six little hills are located forming a regular hexagon as well. The features of porous alumina surface detected during recognition are shown in figure 15. The primary statistics of the found features are presented in table 1.
It was observed during the experiments that the porous alumina surface adsorbs air moisture very well, i.e., it is hydrophilic. The presence of a water layer on the surface makes resolution of the pores practically impossible. An attempt was made to prick this layer by increasing pressure force on the surface with a cantilever probe (the cantilever from NT-MDT Co. was used) and then to measure the solid topography, which did not however lead to good results. Neither was it successful obtaining a high contrast by increasing amplitude of free cantilever oscillations. Observation of the pores along with a finer structure-a rim of six hills-only became possible after the adsorbed water layer has been removed by heating the sample to temperature 70-80 • C for a few minutes. At room temperature and relative humidity of 50-70%, the 'lifetime' of the surface once the heater is switched off is as little as 10-15 min. After that time the surface again gets covered with a water film making it difficult to perform a high resolution scan. It was also discovered that the scanning process itself would stimulate moisture redistribution/condensation on a porous surface since repeated scanning of areas having once been scanned yielded a worse contrast in comparison with the neighbouring areas where the scanning was carried out later for the first time. A similar effect on a hydrophobic HOPG surface was observed in [11].
As the employed microscope had no heater imbedded into the sample holder, the heating was done outside the microscope. Then, the hot sample was mounted in the microscope and the scanning was carried out during the next 5-10 min. Because of the strict time restrictions, the counter-scanning was implemented immediately after probe engagement and finding a suitable area of the porous surface.
Since a long-time microscope relaxation is impossible, besides the creep appearing when scanning and the thermodrift typical of the microscope being used, the following additional drifts took place during the alumina topography measurement: approach-caused creep of the Z manipulator; creeps of the X, Y manipulators as a result of probe offset towards the scan start position; as well as the thermal drift arising from nonuniform cooling of the holder and the sample after their heating to the mentioned temperature. Hardening of the double-sided adhesive used for sample mounting purposes probably also led to an additional motion.
Acting together, the mentioned negative factors distort the image so strongly that the counter scan contains only half (!) of the surface area obtained in the direct scan. Despite such significant distortions, the given images could be corrected with a quite acceptable error at least by applying the linear model, as will be shown later on.

Linear drift correction.
In figure 16, marked with a '+' sign are those alumina surface features that were detected in both CSIs during the iterative search for feature pairs. Driftinduced maximal relative measurement error in the lateral plane δ max is 24.6% for the first CSI pair and 24.4% for the second CSI pair. In spite of the significant distortions, the feature pairs were revealed without errors. By inserting the found coordinates of the features in expressions (5), LTCs K x , K y were calculated (see table 2). It should be noted that the relationships between the inequalities for K x coefficients of CSI pairs given in section 2.2 got violated apparently because of a considerable change in x component of drift velocity (take note of the strong curvature of the regression surface given below). LTCs K z found by formula (12) during vertical correction are given in table 2.
The corrected topography and the residual feature mismatching for the first and the second CSI pairs are shown in figure 17. The lateral and the vertical mismatch errors as well as the relative correction errors are presented in table 3. In figure 18(a), the result of linear drift correction by two CSI pairs is presented; corresponding residual mismatch is shown in figure 18(b).
To avoid image blur causing a loss of useful topographical data, the matching of the corrected CSIs (figure 2, positions 8, 9) may be skipped. In that case, the microscopist should decide which of the four images to select that best characterizes the surface under investigation.
In figure 19 are presented the dependences of the average of mutual mismatches of the corrected positions of the feature pairs on the pair's serial number, which were obtained during the iterative search for feature pairs. Drift velocity components, which were active while counter-scanning the alumina, are given in table 4. Figure 20 shows some example regression surfaces drawn through local lateral feature displacements and local topography height displacements of the first CSI. Nonlinear correction of the first and the second CSI pairs as well as corresponding residual mismatchings of the features are given in figure 21. Nonlinear drift correction carried out by two CSI pairs is presented in figure 22. As compared to the linear correction, the nonlinear correction provides a much sharper image of the porous alumina surface (compare figure 21 with figure 17, also compare matching errors indicated in table 3). In figure 23(a) is shown free of drift distortions topography of porous alumina reconstructed during nonlinear correction of feature segments. Here, it also can be easily seen that the image built is sharp not only in the central part, as was the image obtained above in figure 18(a), but all over the scan area. By using the data given in figure 23(a), the mean distances between hills (44.9 nm) and pores (69.9 nm) were determined. A model of the porous alumina surface is shown in figure 23(b).

Discussion
Despite a substantial thermal drift component produced by sample heating and creep caused by the probe moving to the starting point of the scan, the main contributor to the image distortion of the alumina surface is the piezomanipulators' creep generated during line scanning. That is proved by the distortions along the 'slow' scan direction typical of that kind of creep: the direct image is contracted and therefore contains more features than the counter image, which is in contrast stretched out (see figures 14-16 as well as tables 1, 2 and 4).
The observed picture is a complete inverse of the one obtained on graphite (see figures 3-5 as well as tables 1, 2 and 4). Therefore, the thermodrift should be considered as the main distortion factor at the atomic scale, as has been initially supposed. The described distortions of graphite and alumina were qualitatively observed repeatedly for different samples, scan sizes and scan velocities (sizes and velocities were set to half or twice as much as the values mentioned above).
The fact that all four images of the porous alumina surface turned out to be matched precisely enough in a single image (see figure 18(b) and table 3), despite significant residual mismatching in the pairs (see figures 17(b), (d)), points out that the distortion is developed the same way in the first and second CSI pairs (see tables [1][2][3][4]. The obtained result proves experimentally the validity of the correction scheme applying two CSI pairs.
Because of built-in redundancy, FOS productivity is considerably less than the productivity of conventional  scanning whereas the precision of the measurements in contrast is much higher [6]. For a number of practical tasks very high measurement precision is not required; instead what is important is the functionality provided by the FOS approach retaining a satisfactory scanning productivity. In that case, the skipping operation should be declined, the surface being collected aperture-wise instead of segment-wise. Implementing the aperture counter-scanning followed by drift correction by the suggested methods, it is possible to obtain the true topography in the aperture and so determine the true distance between the current and the next features of the chain. Possessing the corrected apertures and the relative distances between them, it is easy to reconstruct the whole surface.

Conclusions
The suggested method ensures drift correction based exclusively on the information contained in CSIs; no information about drift velocity is required.
Application of the counter-scanning allows for (1) easy detection of feature pairs due to the existence of a CP whose neighbourhood actually contains the same topography; (2) increasing drift correction precision due to a greater difference in position coordinates of the features that compose a pair; (3) ensuring a preset measurement error within the bounds of a certain image area in accordance with the actual change in drift velocity that occurred during the counter-scanning; (4) decreasing the creep produced by the piezoscanner during FOS by means of inducing the counter creep while scanning topography segments/apertures.