Fast image registration algorithm for automated inspection of laser micromachining

We investigated the possibility of fast quality control of laser micromachining of surface based on optical measurements. The main objective is to match CAD-model and 3D geometry of the machined surface. We found that standard matching algorithms have either low performance or are ineffective at high noise level or at presence of different types of distortion. The algorithm based on Ciratefi algorithm, which previously wasn’t used for image registration, was developed. The performance of algorithm was increased by using an iterative search of optimum based on image pyramid. The testing of algorithm on height maps of objects formed by laser methods showed its high accuracy and performance.


Introduction
Laser micromilling is one of the advanced techniques to obtain D micron-resolution products. Material machining is performed in layers using focused beam laser based on CAD-CAM methods. Typical structures are formed in material volume of ~ mm 3 with a resolution of 1-10 microns; however, the machined surface may have a high wall angle to Z-axis and a high degree of roughness (> micron). In order to obtain structures with high quality characteristics the inspection of fabricated objects is required [1][2][3][4]. The process of automated inspection consists of three steps: measurement and reestablishment of 3D geometry, regeneration (matching) of CAD-model and reconstructed 3D geometry, retrieval of quality characteristics and geometrics benchmarking [5]. Required accuracy and performance of object measurements in the area of 11 mm 2 is 1 micron and 1 minute, respectively. Such characteristics are now provided only by optical methods. The purpose of this paper is to create an algorithm and software tools for quality control of laser micromachining. To do this, we had to develop algorithm to match height maps of CAD-model and the machined surface. As height maps we call 16-bit greyscale images, in which intensity of pixel determines a depth of laser removal of material. Height maps contain several million pixels and have submicron resolution. It is necessary to ensure registration of such height maps with an accuracy of 1 pixel during 1 min that will enable to control quality of lots of test objects during the period compared with the time of their formation (less than 1 hour.)

Problem description
Let us consider characteristics of structures formed in surface laser micromachining and particular qualities of their measurement using optical methods. In contrast to machining (milling), the depth of a laser removed layer depends not only on axial movement of machining tools with regard to rough pieces, but also on a number of process parameters (pulse duration, wavelength, pulse energy, etc.) Therefore at not optimal micromachining modes, the thickness of the removed layer may vary nonlinearly from layer to layer, and thus, the depth of structural elements may differ from that one set in CAD-model. Besides, the machined surface may contain defects (recast layer, debris), i.e. interconnected spaces with the area of more than  μm 2 , the average depth of which is different from the depth of the respective layer by more than 10 microns. When measuring structures with high wall angles to Z-axis using optical methods, the returned and scattering radiation is not completely gathered on a photosensor due to input aperture limits that results to significant deterioration of signal-to-noise ratio. Moreover, the registered signal may exceed a dynamic range of the optical sensor due to the fact that material reflection coefficient can be changed by several orders, when subjected to laser emission (for example, in machining, on clear glass or smooth metal surfaces.) We may use image registration techniques for registration of CAD-model and reconstructed D geometry [6]. In image registration we conceive that scale factors are known (from hardware settings), however displacements in three axes and rotation in XY-plane are unknown since it's not always possible to place a measuring sample in the same way as it was placed in laser machining. In certain lateral displacements and rotations in XY-plane the axial displacement is performed in such a way as to match corresponding points of the un-machined or the least-machined surface that enables to determine depth mismatch. Thus, the problem of image registration has three independent variables: two displacement parameters and the angle of rotation. Figure 1 gives height maps of the test CAD-model and the structure formed in the process of laser milling and measured by a confocal laser scanning microscope. For height maps H 1 (x,y) and H 2 (x,y) the problem of image registration is formalized as follows: it requires to discover spatial transformation g and light intensity transformation f of points (x,y) so as to have the following: In view of aforementioned properties of laser micromachining and the measurement method, the problem of image registration has some peculiar features. First, images may differ by their light intensity and contrast. Secondly, there are distortions at the height map of the machined surface, i.e. interconnected areas and contours, which are not present at height map of CAD-model. Thirdly, height map of the machined surface is noised and contains statistical outliers, the occurrence of which is associated with measurement errors. Therefore, for height maps of CAD-model and reconstructed D geometry the equation (1) is approximately fulfilled.
For accurate image registration it is required to develop algorithm which is stable to the above mentioned image properties. In this case, it's not always possible to add special marks around the model to simplify evaluation of mismatch parameter; therefore image registration is to be performed through search of matching either between fragments or between image characteristic points.

Image registration techniques
For image registration of CAD-model and reconstructed D geometry in tasks of automated inspection of machining operation, the standard [7,8] is supposed to be the Iterative Closest Point algorithm (ICP) [9] and its optimizations [10]. Two point clouds, i.e. of CAD-model and reconstructed D geometry, are supplied to algorithm inputs. Matching between points of CAD-model and reconstructed D geometry is set up at each iteration by the nearest-neighbor criterion; a summary distance between two point clouds and transformation parameters, which allow to minimize mismatching, are computed using singular value decomposition of cross covariance matrix; one of the point clouds is transferred and then the next iteration is performed, and so forth, until changes between iterations become minimal. The algorithm converges to global extremum if point clouds are initially roughly converged. The simplest approach for determination of mismatching parameters consists in correlation algorithms. To evaluate registration accuracy of images I 1 and I 2 we use various measures such as the cross-correlation function: , .
x y R x y I x y I x x y y In all possible image displacements (x,y) the summary of light intensity of respective points (x,y) is computed. The maximum correlation value corresponds to the best matching.
To evaluate mismatching parameters we also use Fourier methods. The phase-correlation method allows us to determine lateral displacement. Image rotation in XY-plane and its scaling duplication may be determined by the method [11] and Fourier-Mellin transformation. There are a lot of different methods, in which characteristic image points (angles, borders, lines, etc.) and image descriptors, e.g. image moments, are used in matching. Lateral displacement may be determined using zero-and first-order moments.
The image X-rotation angle () can be calculated by formula: where x y , are mean values of light intensity in x and y directions. In [12] a research of registration accuracy of portrait images has been conducted using Fourier-Mellin algorithm and the moment algorithm. Mismatching parameters are as follows: X and Y displacement, rotation, scaling duplication. It has been established that Fourier-Mellin transformation enables to precisely compute mismatching parameters and provides better matching than the moment method in small image mismatching (change of rotation angle is within , rescaling -within %, displacement -within 10% from the original image size). In [13] Fourier-Mellin transformation and the moment method have been used to determine image mismatching parameters. It is shown that the moment method has high performance but is sensitive to noises. In [14] the fast correlation matching algorithm has been developed to match space-displaced images with a quasiregular structure. In [15] image registration has been carried out in lines which were preliminarily detected on images using Hough transformation. We have tested the majority of the above algorithms on test objects (including Fig. 1 and Fig. 4) and have analyzed registration accuracy. We found that Fourier methods based on fast computation of the correlation function and ICP (Iterative Closest Point) were ineffective due to image differences in light intensity and contrast. Moment methods and methods based on preliminary selection of characteristic points do not also provide accurate registration since they are sensitive to noises, statistical outliers and different distortions. For example, the key-point detection method SIFT [16] identifies peculiar properties in images of reconstructed D geometry around areas corresponding to defects and establishes false registrations, respectively. To evaluate registration accuracy it was decided to use the normalized cross-correlation function which is presented for two vectors v and w in the following way: where v v v , w w w ; v and w -are mean values v and w, respectively; . -L 2 normа. In our case applying the normalized cross-correlation function makes the algorithm sustainable to mismatching in light intensity and contrast, and less sensitive to noises and distortions due to usage of a large number of statistics. To our opinion, when using the normalized cross-correlation function, the most suitable matching algorithms are Brute-force algorithm and Ciratefi algorithm.

Brute-force algorithm
In Brute-force algorithm we test all possible transformation parameter sets. The image of height map of CAD-model (H 1 ) is moved by a sliding window with regard to the image of height map of D geometry (H 2 ), and at points (x, y) we calculate normalized cross-correlation (R (x, y)) by formula: where H 1 and H 2 are, respectively, mean values of overlap point intensities of CAD-model and 3D geometry (x,y) in matching. Function (6) is calculated for height maps turned towards each other by different angles with predefined increment angle. After that the maximum value R max is determined which corresponds to the best matching of height maps.
The correctness of algorithm was tested on the test object (height map resolution was over  pixels) shown in Fig. 1 when one-pixel registration accuracy of height maps was required. When pixels displacements is possible, the registration is correct, however, the computational speed is low (> 1 hour on PC Intel Core i5 CPU, 3,2 GHz, 4GB RAM). Computation of the normalized cross-correlation function may be optimized using Fourier transformation and function expansion (6) in the amount of rectangular-shaped functions [17] or due to holding in tables a walking amount and a sum of squares of point intensity, for which two images are registered [18]. However, since height maps of reconstructed D geometry contain statistical outliers, and given that zero-intensity completing points are added when rotated, the aforementioned solutions are considered to be inefficient.

Ciratefi Algorithm
Ciratefi algorithm (Circular, Radial and Template Matching Filter) [19] previously wasn't used for registration of two high-resolution images but was used in the original publication for matching of high-resolution images -scenes with low-resolution imagestemplates. The algorithm does not require allocation of characteristic points, is invariable to displacement, rotation, scaling, light intensity and contrast of images. The algorithm consists of three stages: Cifi, Rafi and Tefi. Let us denote scene points, which can be com- The normalized cross-correlation is computed by formula (5) in different relative cycle shifts of the arrays v and w(x,y). The relative rotation angle (x,y), at which correlation has got its maximum value, is the expected rotation angle of the template at point (x,y). After that the Tefi stage is performed: the normalized cross-correlation is computed for all points {М} by formula (6); in this case computations are performed for each point (x,y) only once for the expected angle (x,y). After that we determine for {М} the correlation maximum, the maximum point position (x m ,y m ) and the template shift, respectively, with regard to the scene. The rotation angle (x m ,y m ) is the image mismatch angle. Specificity of algorithm is the possibility to considerable speedup computations using point filtration after the Rafi stage (removal from the set {М}). To improve registration accuracy at the Tefi stage we calculated correlation not only for the rotation angle x,y), but also for other two neighbor angles: (x,y) - and x,y) + . The testing of algorithm for different height maps has shown that registration accuracy determined by Ciratefi is the same as that one in the Brute-force algorithm, however, the computation time is much less. Nevertheless, registration is made within several dozens of minutes that is unacceptable, therefore it was decided to optimize Ciratefi algorithm.

Optimization of Ciratefi
We have suggested two different possible approaches to improve the performance of Ciratefi algorithm. First, the speedup can be first achieved through rough searching by angles (for example, with the increment  = ), and then through searching with high angle discretization ( = ) for points with the maximum value of the normalized cross-correlation. Secondly, we may search transformation parameters first for low-resolution images and then for high-resolution images using found parameters. We have implemented both approaches and identified that the first approach does not give an advantage since the Tefi stage requires the longest time in Ciratefi algorithm; however, its performance is essentially independent of angle discretization since the rotation angle is determined for each point at the Rafi stage. In this case the speedup at the bined with a template center, by a set {M}. At the first stage -Cifi -the template center and all points of the set {M} are announced in turn as centers of concentric rings. Through simple mathematical computations the expected scaling factor is determined for each point {M}. We will skip this stage since in our task scaling factors are known from hardware settings. We will show the operation principle of the second and third stages for height maps given in Fig. 2 and Fig. 3. The height map of CAD-model is the template, and the height map of reconstructed 3D geometry is the scene. For de-scriptive and simplicity reasons we'll assume that the required registration accuracy is 1 pixel and 15 degrees. At Rafi stage, radials with the angle increment  = 15 degrees, which corresponds to the required registration accuracy, are drawn for the template center. The same principle is used to draw a set of radials for each point of the set {M} highlighted in Fig. 3 with a black square. Summary light intensity is computed for each radial line. Thus, we obtain the set of 24 values of light intensity v for the template and an array w(x,y) for each scene point.
Rafi stage is compensated by the time required for optimization stage. Experimental test has shown that performance of registration is improved when using image pyramid. Image pyramid is a sequential order of images, in which each subsequent image is obtained from the preceding one by its filtration that enables to suppress high-frequency noises and also by decimation of each second report. Therefore, image pyramid L(x,y,i) may be presented in the following way: S( , , ) ( , , ) ( , ), ( , , ) ( , , ), where * -is a convolution operator, i -is an image index in the sequential order, L(x,y,0) -is the original image; S(x,y,i) -is the result of image convolution with a kernel h(x,y). Gaussian function has been used as a convolution kernel; the minimum image size in the sequential order is  pixels. Image pyramid is constructed for height maps of CAD-model and reconstructed D geometry. And then, starting with the lowest resolution image, optimal mismatch parameters are searched by Ciratefi algorithm. A range of interest (ROI) is determined all round optimum (e.g.  pixels) for which we search parameters on higher-resolution images, etc. Mismatch parameters found for initial resolution images are considered as required.
Note that the use of image pyramid to speed up computations in image processing and registration tasks is a well-known approach [20; 21]. However, it previously wasn't used in combination with Ciratefi algorithm for fast registration of noisy images at presence of distortions.

Experimental information
To form test objects by CAD-model we used a system of laser micromachining based on the complementary scanning principle [22]. The system contains Multiwave impulse fiber laser, a scanning head with scanners CT  and telecentric lens, a portal-frame mechanism for template positioning and MarkKey software. Test CAD-model consists of 25 steps the depth of which is gradually increased in spirals. The size of e It is shown that standard algorithms based on image moments ach step is 1001006 μm 3 . Material is removed in layers so that the central step is not stripped off by laser, the step to the right of the center is removed on the first layer, and the lower right step is removed on the 24 th layer. The test object has been formed on a magnesium-based layer and measured by the confocal laser scanning microscope Carl Zeiss LSM  with axial increment between images of optical sections of  microns; D geometry of the structure has been reconstructed using the standard mass center method [23]. Fig.  4 shows CAD-model and reconstructedD geometry; resolution is 16181618 pixels and  pixels, respectively. Registration was performed using the developed algorithm with the following parameters: filtration keeps 50% of points after the Rafi stage, the size of search area was reduced to the area of  pixels round the maxi-mum, registration accuracy  = 1 degree and 1 pixel. As a result of the performance of algorithm we found that the following conversion parameters are the most suitable: x = , y = , =  degrees. Fig. 5 shows registration results of height maps of CAD-model and reconstructed D geometry; a depth deviation scale is given in the top right angle. Positive and negative values indicate, respectively, how much less and more material than required has been removed.

Fig. 5. Registration of height map of CAD-model and reconstructed D geometry.
Analysis of registration results allows a system operator to determine the presence of depth deviations and sizes of areas with defects in laser micro-machining. Fig. 5 shows that for first steps the layer depth corresponds to the depth defined in CAD-model, and starting with the seventh layer the depth is less than specified and grows nonlinearly. Besides, it is clear that a boundary between layers needs to be improved. Based on registration results qualitative characteristics of the structure can be identified (e.g., Fig. 6, 7) which shows a high practical value of the developed solution.

Benchmarking of registration algorithms
We have compared the developed algorithm, Ciratefi algorithm and the Brute-force algorithm in registration of height maps given in Fig. 4;-is the angle discretization: (1) The Brute-force algorithm ( =  degrees); (2) Ciratefi algorithm (=  and  degree); (3) The optimized Ciratefi algorithm using image pyramid (=  degree); (4) The Brute-force algorithm using image pyramid ( =  degrees). The size of possible registration area is pixels; in (3)(4) the size of search area is reduced to pixels all round the maximum found at the previous stage. In algorithms (2) and (3)  The testing has shown that when  =  degrees Ciratefi algorithm provides the same registration accuracy as the Brute-force algorithm. The performance is herewith 30 Image processing, pattern recognition 349 times higher. When increasing the angle discretization (= degree), the computation time of Ciratefi is insignificantly increased, whereas it would be increased fifteen-fold for Brute-force algorithms. The developed algorithm provides the same registration accuracy as Ciratefi algorithm, but the performance is higher (less than one minute for registration). It is shown by comparison that pyramid search also accelerates exhaustive search, however the computation time is more than in our algorithm even with low angle accuracy of registration (= 15 degrees).

Conclusion
It this paper we have solved the problem of registration of height maps for CAD-model and the laser-machined surface at high noise level and distortions. It is shown that standard algorithms based on image moments, Fourier transformations and search of key points are inefficient either due to the presence of machining defects and measurement errors. The Brute-force algorithm has low performance. We have developed a new algorithm based on Ciratefi algorithm which previously wasn't used for image registration. Its performance has been improved using image pyramid, i.e. iterative search of optimum -from rough to particular. The algorithm may also be suitable for automated quality control in other solutions of surface micromanufacturing and structuring in which D geometry is reconstructed by optical measurement data.