3D image combination in aviation computer vision systems

Abstract

The relevance of 3D image formation in computer vision systems being applied in aviation is caused by the need to solve many problems using helicopter and unmanned aviation. These problems include mounting of various technical installations with the help of helicopter, rescue mission in complex terrain as well as landing helicopter on unknown sites. The work considers the possibility to use two alternative methods to form 3D images of the underlying surface from the sequence of 2D images obtained with the help of stereo pair placed on aircraft. The first method is based on applying the algorithm of flat image combination to three-dimensional case with the help of homography matrix. The mathematical apparatus developed for the case of flat homography with arbitrary number of key points is generalized to a three-dimensional case. Homography matrix takes into account all kinds of point cloud deformations (shift, rotation, scale change) and allows to obtain the final solution in one step. The second method of point cloud combination is based on applying an iterative procedure to sequentially refine the transformation of one point cloud into the other. In every iteration the key points are found using FAST algorithm. Correspondence between them is established using non-dense optical flow algorithm, optimum estimate of rotation matrix is found in accordance with Kabsch algorithm.

1. Introduction

The article presents the research results of two methods to reconstruct 3D images of underlying surface in the Earth's plane based on the sequence of flat video images of this surface obtained from aircraft (AC) with the help of stereo pair. Image processing is supposed to be performed in an onboard computer in real-time mode. The reconstruction process of underlying surface 3D images is associated with the solution of several interrelated tasks: point clouds sequence formation based on depth maps and matrix with camera parameters; determination of key points on a pair of neighbouring images and establishing the correspondence between them; combination of current and previous point clouds. Conversion of one point cloud to the other involves two transformations: performing point cloud rotation in space and shift along the vector. The transformations required can be performed either simultaneously during the process of transformation with the help of homography matrix, or sequentially, usually after several iterations. and transformation F be determined: relative to the first cloud. In the first and second point clouds there is a certain set of k k  minn, m points that are images of the same points of a real scene. These points are called appropriate. The G  M j , M j kj 1 set of appropriate points is the basis to determine transformation F parameters, which allows two surfaces (two point clouds) to be "sewn together". The required X  FX  includes rotation matrix R  rij i3, j 1 and shift vector t  t x , t y , t z T to ( 1 ) ( 2 )  r11 X  R  X   t   r21   r31 of appropriate points from point clouds being sewn.

Solving the problem of combining two point clouds requires solving two particular problems: 1) forming the set G  M j , M j kj 1 of appropriate point pairs in coincident point clouds; 2) determining the transformation parameters F (rotation matrix R and shift vector t ).

The most difficult one is the first task. The known methods to solve it can be divided into three groups: 1) surface feature identification of compared point sets [1]; 2) usage of spin image signatures [2, 3]; 3) complete point permutation in clouds being compared [4].

Correct determination of the relationships between appropriate parts of point cloud pair allows them to be combined in one step. Otherwise, it is necessary to use multi-step registration procedures, which sequentially reduce error, but require correspondingly large cost of computational time for their implementation.

3. Methods of the research

The authors have developed an algorithm based on the usage of spin image signatures [2, 3] to form the set G  M j , M j kj 1 of appropriate points. A general scheme of the method proposed for a regular grid is as follows: at every point M i0 of the first and second point clouds, the least squares method is used to find estimates aˆ, bˆ of coefficients a, b of equation of the plane z  ax  by passing through the given point. Optimum estimates aˆ, bˆ are found as a result of minimizing the deviation of the required t plane from some set of points M i i 1 , t  4ss  1 taken from the square neighbourhood of point M i0 sized 2s  1 2s  1 . These estimates aˆ, bˆ are necessary for normal vector n  aˆ, bˆ, 1 formation at every point. All normal vectors are normalized n0  n  .

 n  elements of which are triples of numbers  j , j , j  - point descriptors M i0 . Where: Each point M i0 from the first and second point clouds is associated with vector  j , j , j tj 1 ,  j - the corner between normal vector ni0 of the current point M i0 and vector n0j j 1, 2,...,t of 0 point M j from point neighborhood M i0 ;  j algebraic projection of vector

M i0 M j on vector point M i01 from the first cloud, and  j2, j2,  j2 tj1 - the set of descriptors of some arbitrary point M 2 from the second cloud, checked for correspondence to point M 1 . For each point M 2 i i0 i from the second cloud, the proximity condition for each parameter is checked first, namely, in the cycle of j from 1 to r, the following conditions are checked:  j1  j2  1 ,  j1  j2   2 ,  j1   j2   3 where 1, 2 , 3 are the given numbers. If at least one of three conditions is not met for at least one value of j index, then transition to a new point of the second cloud occurs. Otherwise, the sums are found   M i2  jt1 i01j  ij2 ,   M i2  jt1 i01j  ij2 ,   M i2  jt1 i01j   ij2 , i  I.

For each sum in the composition ( 3 ) there are points minimum of corresponding sum from ( 3 ) is reached.

M i12  argMmi2in   M i2 , M i22  argMmi2in   M i2 ,

M s2, s  1, 2, 3 , in each of which the

M i32  argMmi2in   M i2  ( 3 ) ( 4 ) If i1  i2  i3  i , then point M 2 from the second cloud is identified as appropriate to point i M i01 from the first cloud. Otherwise, point M i2 is recognized as inappropriate to point M i01 .

The described algorithm to form pairs of corresponding points was tested on a model point cloud pair. The first point cloud V is the surface of a hyperbolic paraboloid

V  x, y, z z  x 2  y 2 , x  10,20 ; y  8,8  .

 50 75 

The second point cloud V ' was obtained by a similar formula and with the same coordinates along Oy axis, and along Ox axis in the interval from 10 to 30.

V   x, y, z z  x 2  y 2 , x  10,30 ; y  8,8  .

 50 75 

The first and second point clouds have a common part (overlapping area) in amounts of 187 points: on coordinate x from 10 to 20, and on coordinate y - from -8 to 8. Both point clouds are shown in figure 1.

The points of the first cloud are marked with circles of blue, and the second cloud - with red squares. The intersection of point clouds is clearly visible in figure 1 as a horizontal strip from combined points of the first and the second clouds in the interval from 10 to 20 along Ox axis.

After point cloud formation, the second cloud underwent shear and rotation transformations. It was shifted by 25 units on coordinate z , so the shift vector t has t  0, 0, 25T form. Then it was rotated relative to each of three axes Oz, Oy и Ox, by angles in accordance with 3, 3, 6 . The rotation was carried out in the following sequence: firstly- in the plane of roll by    30 angle, then in the plane of pitch by   

60 matrix ( 5 ) was obtained.

angle and in the end - in the plane of course by    angle. As a result, a rotation 60

The result of the second cloud points conversion by formula ( 1 ) is shown in figure 2. Now it's necessary to perform appropriate point search in a cloud pair using algorithm ( 3 ) - ( 4 ). Algorithm results are shown in figure 3. points. In this case, as it can be seen below, this number of points is enough to perform proper surface combination.

By analogy with the problem of combining two-dimensional images using homography matrix, the problem of combining (sewing) two point clouds V  xi , yi , zi i1 n  and V  xi , yi , zi in1 can be solved by homography in R 3 space. We will search a matrix of the following type:  h11 h12 h13 h14  H   h21 h22 h23 h24  ,  h31 h32 h33 h314   h41 h41 h41 that would lead to fulfill the following matrix equalities in homogeneous coordinates xi   h11 yi    h21 zi   h31     h41 ( 6 ) ( 7 ) ( 8 )  h11 In this matrix H 44 , submatrix R   h21  h31 h12 h22 h32 hh1233  is responsible for point cloud rotation, and submatrix T  h14 h24 h34T - for point cloud shift.

Having eliminated parameter  from the system of equations ( 6 ), we obtain a system of linear algebraic equations (SLAE) for each pair of appropriate points Mi xi , yi , zi  and Mixi , yi , zi , i  1, k from the point clouds combined  xih11  yih12  zih13  h14  xi xih41  yi xih42  zi xih43  xi ,  xih21  yih22  zih23  h24  xi yih41  yi yih42  zi yih43  yi ,  xih31  yih32  zih33  h34  xi zih41  yi zih42  zi zih43  zi . or in matrix form A GH  B . Where: A - block matrix of 3k 15 order:  Z  A3k15  O  O

O O Z

O O Z

W1 

 W2  W3  with a basic matrix

With

blocks in matrix form

Z k4  xi , yi , zi ,1ik1 ,

W1   xi xi ,  yi xi ,  zi xi ik1 , W2   xi yi ,  yi y,  zi yik1 , W3   xi zi ,  yi zi ,  zi zi ik1; , O - zero matrix of order k  4 , GH -column vector of unknowns of size 15 1: GH  (h11 h12 h13 h14 h21 h22 h23 h24 h31 h32 h33 h34 h41 h42 h43)T ; B

B  x1 x2 ....xk y1 y2 ...yk z1 z2 ... zk T - column vector of free terms of size 3k 1 . Matrices W1, W2, W3 are of k 3 dimWi  k 3, i  1, 2,3 order.

Due to measurement errors, matrix equality ( 6 ) may not be fulfilled for some points from point clouds. We will look for matrix H , which minimizes the distance between points of one set and preimages of the other one. Problems to find homography matrices in transformation of one set of points to another set of points with an arbitrary number of points in them, both in two-dimensional case and in three-dimensional case, are formulated identically Where  - Euclidean norm.

The solution of task ( 7 ) with three-dimensional homography by least square method leads to normal SLAE [5]

AGH  B 2  min .

GH AT AGH  AT B

The right side

of AT B   Z T B1 

Z T B2

Z T B3  Z T Z   O AT A   O  W1T Z  normal 3 T WiT Bi  . i1 

O44 Z T Z

O W2T Z

O

O Z T Z

W3T Z equation

Z T W1  Z T W2  Z TW3  . 3 WiT Wi  i 1  ( 8 ) has estimate of rotation matrix, coincides with a given matrix chosen   105 . Submatrix Ht  2.581010 7.16 1012 25.0 also gives almost accurate estimate of shift vector t. The result of point cloud sewing through the use of obtained homography matrix is almost ideal.

Let's now explore the influence of errors while determining point coordinates in point clouds on the accuracy of surface combination with homography matrix. We shall add normally distributed random component with zero mathematical expectation and given noise level to coordinate z in each point cloud. Let's perform two experiments. In the first   0.01 , and in the second   0.1 , Homography matrix

R , ,  ( 5 ) within the accuracy 0.104 0.993 24.985  2.70 105 0.0001 1  found from noisy data at noise level  0.01 , insignificantly varies in comparison with homography matrix ( 9 ) obtained from undistorted data. In rotation submatrix H R , the values of elements have changed by no more than 0.001. Elements of the fourth row of the matrix have changed to a greater degree, but not so much to significantly affect the results of point cloud transformation.

Errors of surface combination substantially increase with increasing noise level. Homography matrix found at   0,1 error level in the original data has the following form

Substantial distortions appeared in submatrix H R that is responsible for the rotation. Elements h13, h31 values increased fourfold in comparison with ( 9 ), diagonal elements of this submatrix changed absolutely in the range from 0.062 to 0.262 (by 26%). Some elements of the fourth row of the matrix also changed. They also affect the result of point clouds transformation. Transformed point cloud distortions can be observed in figure 4 in point cloud projection on Oxy plane.

The research allows to make the following conclusion: 3D Surfaces combination based on one point cloud transformation to the other using homography matrix gives good results at precisely known point coordinates in point clouds and requires low computational costs. However, even low level of errors in point coordinates in clouds being combined can lead to significant distortions of combining results.

Unstable results of 3D image combination with the help of homography matrix force to look for other combination methods. A method intended to form 3D images of underlying surface in an onboard computer of aircraft must meet serious requirements for both speed and precision of results. These requirements are generally satisfied by sequential refinement method to transform one point cloud to the other. In the algorithm proposed, point cloud formation is made from a sequence of 2D image pairs of the underlying surface obtained from stereo pair placed on aircraft. Point cloud sewing is performed through the use of Kabsch algorithm [8]. The algorithm consists of several interrelated blocks, each of which solves its own independent task, namely:

- Block for depth maps formation from a sequence of 2D images pairs of the underlying surface from the stereo pair;

- Block for point cloud sequence formation based on depth maps and a matrix containing camera parameters; - Block for key points determination on i-th and (i +1)-th images using FAST algorithm [9]; - Block to establish the correspondence between key points in a point cloud pair using LucasKanade non-dense optical flow algorithm [10];

- Block to combine i-th and (i +1)-th point clouds with the search of rotation matrix R using Kabsch algorithm.

In the process of combining the images containing uniform surfaces (table, floor, flat areas of the earth's surface without pronounced features, etc.), the standard version of ICP algorithm gives poor results. On such surfaces, it is difficult to find key points and establish correspondences between them. Using FAST key point detector and optical flow algorithm in a stereo pair scheme allows finding key points on adjacent frames of flat images and establishing the correspondence between them in twodimensional space. The correspondences established between points on flat images allow finding correspondences between key points of points in three-dimensional space.

The choice of FAST algorithm as a detector of key points is caused mostly by its high performance. Unlike other detectors using partial derivatives calculation, this algorithm is based on direct comparison of investigated pixel brightness and therefore has a minimum number of computational operations. FAST algorithm constructs a circle of 16 pixels inscribed in a square with the side of 7 pixels. The brightness of each of these pixels is compared with the brightness of a central pixel and is referred to one of three classes (darker, comparable, brighter ones). As a result, set P of all checked pixels is divided into three subsets, the decision tree is constructed in accordance with algorithm [9]. The decision tree built is used to identify singular points.

Depth map is constructed on the basis of determining displacement value between corresponding pixels of left and right images. The search of a corresponding pixel in a row is made by a window of specified size passing through image and determining objective function maximum. To improve the quality of a depth map received, cameras are calibrated, image rectification is performed [13], and image filtration and improvement for example using the method described in [14] can also be made. Figure 5a shows depth map in pseudo coloring, where values of disparities are respectively large (red) and small (blue) from red to blue.

On the second step, point clouds Pi and Qi 1 are formed for i-th and (i+1)-th frames by displaying corresponding depth maps to three-dimensional space using matrix H containing camera parameters. Figure 5b shows a three-dimensional point cloud formed at this stage.

The third step of the algorithm is designed to search multiple key point pairs. Key points on i-th image are found by FAST detector. For key points found, using optical flow, appropriate points on (i +1)-th image (Figure 5c) are found. In the method considered, optical flow is realized using LucasKanade algorithm [10] and its enhancement in [15].

For neighboring frames of time arrival from stereo pair, ordered arrays Pi0  pi01, pi02,..., pi0k and key points Qi01  qi01,1, qi01,2 ,...,qi01,k 

are formed. In these arrays, point pi0j  x j , y j , z j  with number j from i-th point cloud corresponds to point qi01, j  xj , yj , zj  with the same number from cloud i+1.

The next 8 steps of the algorithm implement an iterative procedure of two point clouds combination in accordance with the logic of ICP algorithm. In this case, rotation matrix R is found using Kabsch algorithm [8, 11]. The following computational procedures are sequentially realized. Pi0 arrays mass centers are calculated Qi01 pi0  1 k 0

 pij , k j 1

0 qi 1  1 k 0

 qi 1, j k j 1 ( 10 )

The coordinates of mass centers found are subtracted from the coordinates of each point of a corresponding array p1ij  pi0j  pi0 1 0 0 qi1, j  qi1, j  qi1 j  1, k , or, in a matrix form

Then, a covariance matrix is calculated

On the last 11th step of the algorithm, approximation error of converted i  1 point cloud Qˆi 1 to i-th cloud Pi0 at r  1 iteration is calculated  ir11  Qˆi1  Pi0  k xˆ j  x j 2  yˆ j  y j 2  zˆ j  z j 2 .

j 1

If in the equation  ir11  ir1   is fulfilled, then the decision to complete an iterative process of successive point cloud approximation i  1 to i-th cloud is made. Otherwise, calculation cycle is repeated.

Let us consider the example of applying the described algorithm for the formation and combination of point cloud sequence based on video images of the underlying surface obtained with video cameras pair placed on hexacopter suspension. Figure 6a shows two point clouds (arrays Pi0 and Qi01 ) obtained at the output of the algorithm's second block. The result of their combination is shown in figure 6b - the result of their combination. Figure 7 shows point cloud “sewing” results obtained from two stereo sequences for two variants of shooting the same section of the underlying surface. (11) (12) (13) (14) (15) (16)

4. Conclusion

The computational formulas are obtained and the study on model images of 3D image combination algorithm based on homography matrix usage is carried out. On the basis of this approach, the algorithm to find appropriate points in combined clouds is proposed. The algorithm is based on spin image signatures application. Homography matrix allows taking into account all kinds of possible point cloud transformations. The algorithm can operate in real time and it provides good quality of surface combination in the absence of errors in coordinate assignment of points in clouds being combined. At the same time, in the process of research this method was shown to be very sensitive to perturbations of point coordinates in point clouds being combined.

A modified version of ICP iterative algorithm, which provides the solution to the task of 3D images of underlying surface reconstruction in real time with good quality, is proposed. Main differences between the proposed variant of an iterative algorithm and its classical variant are the application of Lucas-Kanade algorithm of non-dense optical flow together with FAST algorithm to form a set of key point pairs as well as Kabsch algorithm to find rotation matrix.

Stein F and Medioni G 1992 Structural Indexing: Efficient 3-D Object Recognition IEEE Transactions of Pattern Analysis and Machine Intelligence 14( 2 ) 125-145 Bergström P and Edlund O 2014 Robust registration of point sets using iteratively reweighted least squares Computational Optimization and Applications 58( 3 ) 543-561 DOI: 10.1007/s10589-014-9643-2 Bouguet J Y 2017 Camera Calibration Toolbox for Matlab Computational Vision at the California Institute of Technology (Access mode: http://www.vision.caltech.edu /bouguetj/calib_doc/) (16.10.2017) Nikonorov A V, Petrov M V, Bibikov S A, Kutikova V V, Morozov A A and Kazanskiy N L 2017 Image restoration in diffractive optical systems using deep learning and deconvolution Computer Optics 41( 6 ) 875-887 DOI: 10.18287/2412-6179-2017-41-6-875-887 Lyubutin P S 2015 Development of optical flow computation algorithms for strain measurement of solids Computer Optics 39( 1 ) 94-100 DOI: https://doi.org/10.18287/0134-2452-2015-39-194-100

References

[1] Stein F and Medioni G 1992 Structural Indexing : Efficient 3- D Object Recognition IEEE Transactions on Pattern Analysis and Machine Intelligence 14 ( 2 ) 125 - 145

[2] Johnson A and Hebert M 1997 Surface Registration by Matching Oriented Points Proc. 3 DIM 145-152

[3] Chernikov I S and Krylovetsky A A 2009 Three-dimensional reconstruction and recognition by sets of range-finding data Computer modeling 2009: works of international scientific and technical conference (St . Petersburg) 284 - 287

[4] Chen C , Hung Y and Cheng J 1999 Ransac-Based Darces : A New Approach to Fast Automatic Registration of Partially Overlapping Range Images IEEE Transactions on Pattern Analysis and Machine Intelligence 21 ( 11 ) 1229 - 1234

[5] Demidovich B P and Maron I A 1966 Fundamentals of computational mathematics (Moscow: Science) 664 p

[6] Voevodin V V and Kuznetsov Yu A 1984 Matrices and Calculations (Moscow: Science) 320 p

[7] Andrushevsky N M 2008 Analysis of stability of solutions of systems of linear algebraic equations (Moscow: MAKS Press) 71 p

[8] Kabsch W 1976 A solution for the best rotation to relate two sets of vectors Acta Crystallographic 32 ( 5 ) 922 - 923

[9] Nikolajczyk K and Schmid C 2004 Scale and affine invariant interest point detectors International Journal of Computer Vision 60 ( 1 ) 63 - 86

[10] Bouguet J-Y 2017 Pyramidal Implementation of Lucas Kanade Feature Tracker Description of algorithm Intel Corporation Microprocessor Research Labs (Access mode : http://robots. stanford.edu>cs223b04/algo_tracking) ( 11 . 10 . 2017 ).