A New Dynamic Programming Method for UAV 2 Image Stitching Using Optical Flow 3

Dislocation is one of the major challenges in unmanned aerial vehicle (UAV) image 15 stitching. In this paper, we propose a new dynamic programming for seamlessly stitching UAV 16 images using optical flow. Our solution consists of two steps: Firstly, an image-matching algorithm 17 is used to correct the images so that they are in the same coordinate system. Secondly, a new 18 dynamic programming algorithm is develop based on the concept of a stereo dual-channel energy 19 accumulation using optical flow. A new energy aggregation and traversal strategy is adopted in our 20 solution, which can find a more optimal seam line for adjacent image stitching. Our algorithm 21 overcomes the theoretical limitation of the classical Duplaquet algorithm. Experiments show that 22 the algorithm can effectively solve the dislocation problem in UAV image stitching. Beyond that, 23 our solution is also direction-independent, which has more adaptability and robustness for UAV 24 images. 25


Introduction
Image stitching has a long history starting in the early days of computer vision and photogrammetry.With the rise of UAV remote sensing technologies, this research has become paramount to many applications based on UAV survey including 3D reconstruction, ecological farming, disaster emergency management, and photography activity.These are due to UAV remote sensing technologies' strength in low cost, high speed, and easy access [1][2][3][4][5] .However, its low flight altitude and the camera perspective constraints, the coverage area of a single UAV image is small.In many applications mentioned above, in order to expand the image coverage area to capture more information from the target area, multiple UAV images are collected, leading to the need of stitching multiple images to form a mosaic image.What's more, high-altitude wind has a significant impact on the UAV platform due to its light weight [6] , problems such as irregular image overlapping and uneven image exposure are introduced into the adjacent images.Therefore, images captured from an unstable UAV platform will lead to a vulnerable stitched image with ghosting, blur, dislocation, and color inconsistence.Overall, there are many challenges with the problems about the state-of-the-art image stitching methods.In response to these difficulties, this paper proposed a new algorithm to realize the seamless stitching of UAV images through a comprehensive theoretical analysis of the stitching methods based on dynamic programming and optical flow.The algorithm especially solves the dislocation and ghost problem cause by selecting the optimal seam-line in the building-intensive areas.In this new method, we introduce the optical flow to construct the energy function for seam-line searching, it can better considerate the image structure information into the seam line optimization.Its basic idea is to determine the geometric errors introduced by perspective errors, camera distortions, and radiation errors by analyzing the mapping relationships between the left and right images, then using these errors to aid the seam-line search process.

Related Work
There are various methods for seamless stitching of UAV remote sensing images have been investigated [7][8][9][10][11][12][13][14][15][16][17] .Among them, seam line-based methods are intended to reduce grayscale and geometric differences, they look for the least-cost seam-line in the overlapping region of adjacent images by constructing an energy function.In this paper, we will focus on the seam-line search methods based on dynamic programming and optical flow.

A．Methods based on dynamic programming
This is a kind of mainstream image stitching methods.These methods in [11-13] focusing on the energy difference between the images and its effect is superior, but there still exists some problems.
For example, dynamic programming-based methods in [11-13] adopt Dijkstra's shortest path algorithm to search for the optimal seam-line, which addresses the ghosting and dislocation problems because of the movements of the objects and registration errors, but it suffers from low search efficiency and complex weight determination.The ant colony method in [14] also is based on dynamic programming, which can evade the area where the color contrast is larger on the image, while it will easily lead the search processing to the local optimum due to its sensitivity to ants' number.In addition, there are some other methods based on the snake model [15], and some based on morphological model [16, 17].Although these methods can almost guarantee the consistence of the geometric structure and evade the phenomenon of ghosting in the overlapping regions under certain conditions, they are unable to guarantee that ghosting and seams will be overcame at the same timeespecially when there is a significant brightness difference in the adjacent images.Simultaneously, these methods are unable to achieve satisfactory results when there are rich texture structures in image pairs, registration errors, and radiation brightness differences.Furthermore, most of the current seam-line generation methods based on dynamic programming theory rely strongly on image direction that leads to a low robustness with energy functions.

B．Methods based on optical flow
Optical flow is the pattern of apparent motion of objects, surfaces, and edges in a visual scene caused by the relative motion between an observer and a scene [18][19] .The American psychologist Gibson introduced the concept of optical flow in the 1940s [20] .Sequences of ordered images allow the estimation of motion as either instantaneous image velocities or discrete image displacements [21] .
David and Weiss introduce gradient-based optical flow [22] . of measurements [23] .So far, there are many methods to calculate the optical flow, and these methods have great differences.There is still no systematic classification of the existing methods.Here, we divide the optical flow calculation methods into the following four categories: methods based on gradient, methods based on matching, methods based on frequency, and methods based on Bayesian.
Among them, because gradient-based methods are simple computation and effective, they have been widely studied.Lucas-Kanade method and Horn-Schunk method are representative methods, they are used to calculate the motion of a partial pixels of images (called sparse optical flow) and the motion of all pixels of images (called dense optical flow), respectively [24][25] .The energy function constructed in this paper needs the optical flow value of each pixel in the overlapping area of adjacent images, so, we use the dense optical flow method to calculate them.In the image stitching, using the methods based on optical flow to estimate the camera's motion parameters has the following advantages over the methods based on feature matching: do not need to extract image features, are not sensitive to noise, apply to complex scenes, and can do dense optical flow calculation on the entire image without extrapolation of interpolation.
However, the methods based on optical flow has some weaknesses.Specifically, Feature matching-based methods can be applied to the adjacent image with large difference and correct mark the corresponding points of adjacent images.The methods based on optical flow assume that the change between images is continuous and the difference between adjacent images is very small, which greatly limits the application of these methods.For UAV image stitching, the difference between adjacent images may be very large due to the fast flight and illumination changes of UAV.
Therefore, it is difficult to do UAV image stitching only by the methods based on optical flow.
Nonetheless, the optical flow information of pixels in the overlap area of adjacent images can well provide the structural information of the images, which is conducive to search of the optimal seamline [26] .
In this paper, the optical flow information of the pixels in the overlapping area of the adjacent images is used to construct the improved dynamic programming energy function, trying to find the best seam-line between the adjacent images and realizing the seamless stitching of UAV images.The reminder of this paper is organized as follows.In Section 3 we explain the methodology of our proposed new method based on Duplaquet' method in detail.Experiments and results are described and analysis in Section 4. Section 5 discusses and concludes the paper.

Duplaquet' method
In 1958, Bellman proposed an optimization solution that transforms the multi-stage process into a series of single-stage processes, and thus he invented the dynamic programming method [27] .Based on this, Duplaquet proposed a classic method for searching image seam lines.It ensures that the lengths of all alternative seam-lines are equal, and the seam-line with the smallest accumulated energy value is the optimal seam-line.Formula (1) shows the energy criterion defined in the algorithm [28] : where Cdif (x, y) is the mean value of the gray level difference of the pixel in the overlapping areas between adjacent images, Cedge (x, y) is the minimum gradient value of the pixel in the overlapping areas of the image pair, and λ is a weighing factor, which can be used to adjust the proportion of gray difference and structure difference in energy function.

Problems Analysis of Duplaquet's method
The energy criterion in the Duplaquet' method only considers the horizontal and vertical gradients, and compares the pixels in three adjacent directions near the current pixel, as shown in   the Duplaquet' method cannot guaranteed to obtain the best seam-line by using the classical Sobel operator to calculate the approximate gradient of the pixels based on the horizontal and vertical templates (Formula 2) without considering diagonal directions in the calculation process [29] .Figure 3 is a simulation demonstration of the best seam-line search results under different methods.The red seam-line represents the result of Duplaquet's method, it not only crosses the edge of the houses, but also deviates from the ideal seam-line.This will directly lead to the optimal seam-line susceptible to dense buildings.

Our method
We developed a new dynamic programming method for UAV image stitching using optical flow.
The proposed method searches seam-lines by improving gradient guidance direction, energy accumulation directions (energy aggregation direction, energy traversal direction), and the optical flow value of the pixels in the overlapping area are taken into account.

Gradient Calculation
The Duplaquet method only considers the horizontal and vertical gradient in the energy criterion; it often fails to obtain the optimal seam-line.In order to solve this problem, we develop a new gradient operator based on the classical Sobel operator, by considering eight directional neighborhood information of current pixel and the similarity of its surrounding structure [30] .The new approach of gradient calculation is as follows: (3)   line consists of neighborhood pixels with smallest energy value.In our method, the longest seam line is the optimal seam line.

Calculation of Optical Flow Value of Pixels in the overlapping Area
In this paper, in order to better take into account the image structure information in the overlapping area of the adjacent images, we use the optical flow value of pixels in the overlapping area as a constraint condition for the construction of seam-line energy function.According to Section 2 mentioned, this paper will use a dense optical flow method for optical flow calculation.The H-S method proposed by Horn and Schunck, it is a very popular dense optical flow method, and easy to calculate [25] .The pixel displacement in the overlapping area can be calculated by H-S method.
Formula (4) is its objective function. (4) In formula (4), Eflow is the value of optical flow, u and v is the displacement in the x and y axis directions.T is the reference image, I is the current image.λ is a weight factor.

The Energy Function
Based on the analysis of the theoretical model, we constructed a mathematical abstract expression of the theoretical model.Assuming that image f1 and image f2 are an original image pair to be stitched, the energy function is defined as follows: (5) In formula (5), B(x, y) determines whether the current pixel P (x, y) is in the boundary of the overlapping area, when B(x, y) =1, it means that it is not in the boundary region, and when B(x, y) = 10, it is in the boundary region.δ(*) is the Gaussian smoothing term, which uses the information in the local window to enhance the local influence of the current pixel.f1(•), f2(•) are pending images to be stitched.O is the overlapping area.d(*) represents the gradient function of one of the eight directions.N(x, y) is the energy value of the invalid area, which is the constant term, and the value is 100 times larger than the maximum value of O. Eflow(x, y) is the optical flow constraint item.

Computation Procedure
Because the overlapping area of adjacent UAV images is often irregular, in order to facilitate the calculation, that need to be handled properly.As Figure 6 shows, the irregular overlapping area is in  The image energy matrix is E, it can be calculated by formula (5).The algorithm procedure is as follows: ○ 1 Extract the corresponding points from the adjacent images in order to correct the left and right pending matching images to the virtual unified reference image, so that images are in the same coordinate system.○ 4 Fill the energy matrix E according to the results of ○ 3 .
○ 5 Reestablish two energy aggregation channels: the C1 and C2 matrices, which have the same size as E; each pair of corresponding elements in these two matrices hold two scalar numbers representing the current aggregation value and the current path direction of the seam-line.
○ 6 For the first row of the matrices C1 and C2 assigned with the first row of E as the initial value, and set them corresponding direction as zero.
○ 7 The energy aggregation channel matrices start to make a difference from the second row, which are divided into two aggregation processes from left to right and from right to left.For the energy aggregation channel C1, its aggregation process is from the left to the right; the current pixel only considers the directions of 1, 2, 3, 5, 6, 7and 4. For the energy aggregation channel C2, its aggregation process is from the right to the left, and the current pixel only considers the directions of 1, 2, 3, 5, 6, 7, and 8.
○ 8 When the aggregation is finished, the minimum energy value are found from the last row in C1 and C2 respectively, and then an optimal stitching path is found based on the direction information stored in the matrixes.
In order to ensure that the seam-lines start and end at the intersection points, we select two special intersection points that have the smallest energy value (see that in ○ 3 ), so that the seam-lines can be guided and adsorbed.

Experimental Environment and Data
To verify the effectiveness of our method (Our-flow-DP), we not only utilized the UAV images from different regions with different flight altitudes and cameras, but also compared the experimental results with the classic Duplaquet's dynamic programming method using three search directions (Duplaquet3-DP), dynamic programming method based on Duplaquet using four search directions (Duplaquet3-DP), and the dynamic programming methods from OpenCV (OpenCV-DP). In

Comparison of Stitching Results with Four Groups of Image Pairs
In order to verify the effectiveness of our proposed method, the methods of Duplaquet3-DP, Duplaquet4-DP, OpenCV-DP, and Our-Flow-DP were used to search the optimal seam-lines of image pairs in Figure 7 with irregular overlapping areas.

Comparison of Four Methods under the Condition of Image Rotation
Firstly, the images in Figure 7 were rotated from the horizontal to the vertical.Then, we used the four methods to search the optimal seam-lines for vertical and horizontal images respectively.Figure 13-16 show the results of them.In Figure 13, 15, 16, the partially enlarged pictures illustrated that the optimal seam-lines searched by Our-Flow-DP had basically no change before and after rotation, they always were good at avoiding the ground building and tall trees in the overlapping areas of adjacent images.In Figure 14, our seam-lines changed a little, but they were less affected by the cars and tall trees in the overlapping area, and the directions and movements of the seam-lines basically stayed the same.In contrast, the seam-lines of the other three methods all crossed the edges of the buildings in different places before and after rotation, and the directions and movements of the seam-lines have an obvious changed in Figure 13, 15, 16.In Figure 14, the seam-lines of Duplaquet4-DP is more susceptible to tall trees than Duplaquet3-DP and OpenCV-DP.From the above results analysis, Our-Flow-DP is more direction independent than the other three methods, and it can best avoid houses and tall trees for the best seam-lines searching when there are a large number of buildings and tall tress distribution in the images, this is crucial for finding the most suitable seam-lines for adjacent images.Therefore, due to the specific improvements to the above issues of dynamic programming mentioned in Section 2.3, our method has advantages in adaptability and robustness for different UAV images.The minimum value of our energy function is almost no relationship to the direction of energy aggregation and traversal, and it can better take into account the structural information of the adjacent images.

Efficiency Comparison of Our Energy Accumulation Processing
In the previous experimental results, our method found the most satisfying seam-lines.
Since Our-Flow-DP is based on the classical Duplaquet method, this section will compare the is O(z 8 ), where z is the total number of pixels within the minimum exterior rectangle of the overlapping area, z can be expressed as the product of m and n, m is the width of the minimum exterior rectangle and n is the length of the minimum exterior rectangle.Both m and n are measured by unit pixel.However, because the local energy minima exists in the energy function of the Duplaquet3-DP and Duplaquet4-DP, their results in a lot of time consumption.Therefore, the above assumption is invalid, that is to say: they can not get the same optimal seam-lines.
We selected four experimental image pairs in Figure 7 to verify the above conclusions.In order to speed up the calculation, it is generally necessary to zoom the image at a certain scale.
So, the size of the overlapping area is not same to the size of the original image overlapping area.
The experimental results can be seen in Table 1.The efficiency of our method is more than 37-148 times that of the other three methods.It proved the convergence speed of our energy function is faster than others.In addition, it is further point out that the theory and the results of the proposed method are obviously different with the classic Duplaquet method.The theoretical improvement and experimental comparisons have proven that this paper proposed a global and non-direction optimization method, which not only has the best seam-lines, but also has better time efficiency.
Table 1.Time efficiency comparison of energy accumulation processing

Conclusions
This paper analyzed and evaluated the current mainstream stitching algorithms for the problems of ghosting, dislocation, and seams in UAV image stitching firstly.Then, it selected the essential problems of dynamic programming algorithms for seam line searching, and carried out a detailed theoretical study and a large number of UAV image stitching experimental verifications.At last, this paper developed a new dynamic programming method using optical flow to search the optimal seam lines through several improved key problems of the classical Duplaquet method.
Meanwhile, the superiority and efficiency of the method proposed in the paper are verified by the credible experiments of four image pairs with irregular overlapping area.The stitched results are better than the contrast methods.The proposed method is proved to be direction-independent, more efficient, and more robust.However, there are still some problem with our proposed method, we can see in Figure 14.In the future, we will continue to improve our existing deficiencies to achieve a more perfect and robust method of rapid UAV image stitching.

Figure 1 .
Figure 1.P is current pixel, m, n respectively present the pixel width and height of the overlapping area.When the overlapping area has dense high-rise buildings or tall trees, the seam-lines output from the Duplaquet algorithm are likely across the edge of the buildings due to the inconsistent deformation from image point to the roof point.Therefore, it is easy to find the dislocations in the stitched image leading to a large matching error (as in Figure 2).

Figure 2 .
Figure 2. The stitching results using the Duplaquet method for two datasets.(a) The experimental result with UAV images shot in Wuhan; (b) The experimental result with UAV images shot in Paris.As shown in Figure2.There are two experimental results based on the Duplaquet' method.The seam lines across the houses are prone to ghosting and dislocation in stitched images.In addition, some researchers believe that the energy function is poorly fitted, making it difficult to find the optimal seam-line.For this reason, these researchers attempt to modify the energy function in dynamic programming; however, they overlook the optimality of the corresponding model.This also happens in the algorithms included in the OpenCV library.One of reasons for these problem is that

( 2 )
Specifically, the method of Duplaquet has the following main problems: (1) The gradient guidance direction of the energy function does not support omnidirectional searching.(2) The energy function is direction-dependent; the energy aggregation considers only three directions, and the direction of energy traversal is limited from left to right, as well as from top to bottom.(3) The energy function easy to get a local optimal solution due to the impact of the two factors mentioned above.

Figure 3 .
Figure 3.The sketch map of search results for best seam-lines.

PreprintsFigure 3 ,
Figure 3, the improved purple seam-line 2 is a simulation demonstration result by introduce the fourth direction in energy accumulation.It closer to the ideal seam-line than the red seam-line 1, despite its obvious insufficiency.

Figure 4 .
Figure 4.A Schematic diagram of energy accumulation by improving energy guidelines.Not only the directions of energy aggregation, but also the directions of energy traversal affect the optimal solution of the energy function.Therefore, we redefined the energy criterion and added the new aggregate directions to our dynamic programming algorithm with a stereo dual-channel energy accumulation strategy.It improves the searching scheme of optimal seam line based on energy accumulation.As shown in Figure5, there is a schematic diagram of our optimal seam-lines search strategy, which optimizes the seam-line search criteria by detecting the eight pixels (contain the horizontal direction) surrounding the current pixel.

Figure 5 .
Figure 5. Schematic diagram of our search strategy.In Figure5, P is the current pixel.We redefine the nine related directions surrounding P as follows: 0 (initial invalid direction), 1 (top-left of P for energy channel 1), 2 (top of P for energy channel 1), 3 (top-right of P for energy channel 1), 4 (left of P for energy channel 1), 5 (top-left of P for energy channel 2), 6 (top of P for energy channel 2), 7 (top-right of P for energy channel 2), 8 (right of P for energy channel 2).Seam-line searching is an aggregation process of minimum energy.Each seam-

Figure 6 (
Figure 6(b), it can be extended to a regular area by using the minimum exterior rectangle of the overlapping area.let's say that the overlapping area is m * n.

○ 2
(www.preprints.org)| NOT PEER-REVIEWED | Posted: 15 January 2018 doi:10.20944/preprints201801.0129.v1Define the overlapping area of the adjacent images to be O, the buffer area of O boundary is W (set its width is 20 pixels, and W is an empirical value, the invalid area is N (extend area), and the boundary intersection J. ○ 3 calculation and filled O, W, N according to formula (3), where W Є [1, 10], the closer to the boundary, the larger the value is, and set N = 100 * max (O), J = −1000 * max (O).At the same time, the optical flow item can be calculated by formula (4).

Figure 7 .
Figure 7. Four groups of experimental UAV images: (a) The first image pair; (b) The second image pair; (c) The third image pair; (d) The fourth image pair.

Figure 8 .
Figure 8.A flow chart of our method.
energy accumulation time efficiency of Duplaquet3-DP, Duplaquet4-DP and Our-Flow-DP.Firstly, we assumed that Our-Flow-DP, Duplaquet3-DP and Duplaquet4-DP could find the same optimal seam-lines.Their time efficiency difference can be quantitatively analyzed from the method' complexity.In this paper, the direction of energy aggregation form three aggregation directions increased to eight is mainly an improvement.setting that the time complexity of the Dulapquet3-DP is O(z 3 ), the time complexity of the Dulapquet4-DP is O(z 4 ) and Our-Flow-DP Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 15 January 2018 doi:10.20944/preprints201801.0129.v1

density Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 15 January 2018 doi:10.20944/preprints201801.0129.v1
John, David and Steven provide a performance analysis of a number of optical flow techniques.It emphasizes the accuracy and