SC-AOF: A Sliding Camera and Asymmetric Optical Flow-Based Blending Method for Image Stitching

1. Introduction

Image stitching is a technology that can align and blend multiple images to generate a high-resolution, wide-field-of-view and artifact-free mosaic. It has broad and promising applications in many fields such as virtual reality, remote sensing mapping, and urban modeling. The calculation of the global homography, as an important step in image stitching [1,2], directly determines the image alignment accuracy and the final user experience. However, global homography only works for planar scenes or rotation-only camera motions. For non-planar scenes or when the optical centers of cameras do not coincide, homography tends to cause misalignment, resulting in blurring and ghosting in the mosaic. It can also cause perspective deformation, making the final mosaic blurred and severely stretched at the edges. Many solutions have proposed to solve the problems of parallax and perspective deformation in image stitching, so as to improve the quality of stitched images. But most state-of-art mesh-based [3,4,5] and multi-plane [6,7,8] method are time-consuming and vulnerable to false matches.

In this work, an innovative image stitching method combining sliding camera (SC) and asymmetric optical flow (AOF), referred to as the SC-AOF method, is proposed to reduce both perspective deformation and alignment error. In the non-overlapping area of the mosaic, the SC-AOF method manages to keep the viewpoint of the mosaic same as or one rotation around camera Z axis from those of input images. In the overlapping area of the mosaic, the viewpoint is changed from one input image viewpoint to another, which can effectively solve the perspective deformation at the edge. A global projection plane is estimated to project input images onto the mosaic; After that, an asymmetric optical flow method is employed to further align the images. In the blending, softmax function and alignment error are used to dynamically adjust the width of the blending area to further eliminate ghosting and improve the mosaic quality. This paper makes the following contributions:

• The SC-AOF method innovatively uses an approach based on sliding camera to reduce perspective deformation. Combined with either a global projection model or a local projection model, this method can effectively reduce the perspective deformation.
• An optical flow-based image alignment and blending method is adopted to further mitigate misalignment and improve stitching quality of mosaic generated by a global projection model.
• Each step in the SC-AOF method can be combined with other methods to improve the stitching quality of the other methods. ‘

This article is organized as follows. Section 2 presents the related works. Section 3 firstly introduces the overall method of this article, then an edge stretching reduction method based on sliding camera and a local misalignment reduction method based on asymmetric optical flow are elaborated in detail. Section 4 presents our qualitative and quantitative experimental results compared with other methods. Finally, Section 5 concludes our method.

2. Related Works

For the local alignment, APAP(as-projective-as-possible) [8,9] use weighted DLT (direct linear transform) method to estimate the location-dependent homography and then eliminate misalignment. But if some key-points match incorrectly, the image areas near these key-points may have incorrect homography, resulting in serious alignment errors and distortion. REW (robust elastic warping) [10,11] uses the TPS (thin-plate spline) interpolation method to convert discrete matched feature points into a deformation field, which is used to warp the image and achieve accurate local alignment. TFT (triangular facet approximation) [6] uses the Delaunay triangulation method and the matched feature points to triangulate the mosaic canvas, and the warping inside each triangle is determined by the homography calculated based on the three triangle vertices, so the false matches will lead to serious misalignment. The warping-residual-based image stitching method [7] first estimates multiple homography matrices, and calculates warping residuals of each matched feature point using the multiple homography matrices. The homography of each region is estimated using moving DLT with the difference of the warping residuals as weight, which make the method can handle larger parallax then APAP, but is less robust to the incorrectly estimated homography. NIS (natural image stitching) [12] method estimates a pixel-to-pixel transformation based on feature matches and the depth map to achieve accurate local alignment. In [13], by increasing feature correspondences and optimizing hybrid terms, sufficient correct feature correspondences are obtained in the low-texture areas to eliminate misalignment. The two methods require additional calculations to enhance robustness, but also are susceptible to the uneven distribution and false matches of feature points.

For perspective deformation, SPHP (shape preserving half projective) [14,15] spatially combines perspective transformation and similarity transformation to reduce deformation. Perspective transformation can better align pixels in overlapping areas, and similarity transformation preserves the viewpoint of the original image in non-overlapping areas. AANAP (Adaptive as-natural-as-possible) [16] derives the appropriate similarity transform directly based on matched feature points, and use weights to gradually transit from perspective transform to similarity transform. The transitions from the homography of overlapping area to the similarity matrix of non-overlapping area adopted by SPHP and AANAP are artificial and unnatural, and can generate some "strange" homography matrices, causing significant distortion in the overlapping area. GSP (global similarity prior) [17,18] adds a global similarity priori to constrains the warping of each image so that it resembles a similarity transformation as a whole and avoid large perspective distortion. In SPW (single-projective warp) [19], the quasi-homography warp [20] is adopted to mitigate projective distortion and preserve single perspective. SPSO (Structure Preservation and Seam Optimization) [4] uses a hybrid warping model based on multi-homography and mesh-based warp to obtain precise alignment of areas at different depths while preserving local and global image structures. GES-GSP (geometric structure preserving-global similarity prior) [21] employs deep learning-based edge detection to extract various types of large-scale edges, and further introduces large-scale geometric structure preserving to GSP to preserve the curves in images and protect them from distortion. These methods all require constructing and solving a super large linear equation to acquire the corresponding coordinates after mesh warping.

Based on the above analysis, generating a natural mosaic quickly and robustly remains a challenging task.

3. Methodology

The flow chart of the SC-AOF algorithm is illustrated in Figure 1. The detail about each of its steps are described below.

Step 1: Feature point detection and matching. SIFT (Scale-invariant feature transform) and SURF (speed-up robust feature) methods are generally used to detect and describe key-points from two input images. Using KNN (k-nearest neighbor) method, a group of matched points is extracted from the key-points and used for camera parameter estimation in step 2 and global projection plane calculation in step 3.

Step 2: Camera parameter estimation. The intrinsic and extrinsic camera parameters are the basis of the SC method, and can be obtained in advance or estimated. When camera parameters are known, we can skip step 1 and directly start from step 3. When camera parameters are unknown, they can be estimated by minimizing the epipolar and planar errors, as described in Section 3.3.

Step 3: Sliding camera based image projection. In this step, we estimate the global projection plane first, then adjust the camera projection matrix and generate a virtual camera in the overlapping area by interpolation, and obtain the warped images by global planar projection, as detailed in section 3.1. Misalignment can be found in the two warped images obtained in current step. Therefore, we need to use AOF method in step 4 to further improve the alignment accuracy.

Step 4: Flow-based image blending. In this step, we first calculate the bidirectional asymmetric optical flow between the two warped images; then further align and blend the two warped images to generate mosaic using the optical flow (see Section 3.2 for more details).

3.1. SC: Viewpoint Preservation Based on Sliding Camera

The sliding camera (SC) method is proposed for the first time to solve perspective deformation, and is the first step in the SC-AOF method. For this reason, this section will first introduce the stitching process of this method, and then details how to calculate the global projection plane and the sliding projection matrix required by this method.

3.1.1. SC Stitching Process

In order to ensure that the mosaic can maintain the perspective of the two input images, the SC method is used. That is, in the non-overlapping area, the viewpoints of the two input images are preserved. In the overlapping area, the viewpoint of camera is gradually transformed from $I_1$ to $I_2$.

As shown in Figure 2, the image $I_1$ and $I_2$ are back-projected onto the projection surface $n$, so that the corresponding non-overlapping areas ${\Omega }_1$, ${\Omega }_2$ and overlapping area ${\Omega }_o$ are obtained. Assume that the pixels in the mosaic $I$ are $u^1,u^2,...,u^8$, which correspond to the sampling points $S_1,S_2,...,S_8$ in the projection surface $n$. When the sampling points are within the projection area ${\Omega }_1$ of image$I_1$, the mosaic is generated from the viewpoint of $I_1$. $S_1,S_2,S_3$ are the intersection points of the backprojection lines of $u^1,u^2,u^3$ in $I_1$ and the projection surface $n$. Therefore, $u^i=P_1{S}_i\left(i=\mathrm{1,2},3\right)$, where $P_1$ is the projection matrix of $I_1$. When the sampling points are within the projection area ${\Omega }_2$ of image $I_2$, the mosaic is generated from the camera viewpoint of $I_2$. Similarly, we obtain $S_i$ and $u^i=P_2{S}_i$, where $i$=6,7,8. In the overlapping area ${\Omega }_o$ of $I_1$ and $I_2$, the SC method is used to generate a virtual camera, whose viewpoint gradually transitions from the viewpoint of $I_1$ to that of $I_2$. $S_4$ and $S_5$ are the intersection points of the back-projection lines of $u^4$, $u^5$ in the visual camera and projection plane $n$, respectively. The virtual camera's image is generated from images $I_1$ and $I_2$ using perspective transformation. For examples, pixel $u^4$ of the virtual camera corresponds to pixel $u_1^4$ in $I_1$ and pixel $u_2^4$ in $I_2$, and are generated by blending the latter two pixels.

Global projection surface calculation. In order to match the corresponding pixels $u_1^4$ of $I_1$ and $u_2^4$ of $I_2$, the projection surface $n$ needs to be as close as possible to the real scene point, we can use the moving plane method [7,8,9] or the triangulation method [6] to obtain a more accurate scene surface. Since the SC-AOF method will use the optical flow to further align the images, for the stitching speed and stability, only the global plane is calculated as the projection surface. Section 3.1.2 will calculate the optimal global projection surface using the matched points.

Sliding camera generation. Generally, since the pixel coordinates of $I_1$ and $I_2$ are not uniform, in the mosaic $I$, when $I\left(\tilde{u}\right)=I_1\left(P_{1}S\right)$ in the non-overlapping area of $I_1$, $I\left(\tilde{u}\right)=I_2\left(P_{2}S\right)$ is false in the non-overlapping area of $I_2$, where $S$ is the sampling point on the projection surface. It is necessary to adjust the projection matrix of $I_2$ to $P_2^{\prime }$, so that $I\left(\tilde{u}\right)=I_2\left(P_2^{\prime }S\right)$. As the red camera is shown in Figure 2. Section 3.1.3 will deduce the adjustment method of the camera projection matrix, and interpolate in the overlapping area to generate a sliding camera, and obtain the warped images of $I_1$ and $I_2$.

3.1.2. Global Projection Surface Calculation

The projection matrices of cameras $C_1$ and $C_2$ corresponding to images $I_1$ and $I_2$ are:

$$P_1=K_1\left[I\right|0]P_2=K_{2}R[I|-t]$$

(1)

where $K_1$ and $K_2$ are the intrinsic parameter matrices of $C_1$ and $C_2$ respectively; $R$ is the inter-camera rotation matrix; and $t$ is location of the optical center of $C_2$ in the coordinate system $C_1$.

The relationship between the projection $u_1$ in $I_1$ and the projection $u_2$ in $I_2$ of 3D point $p$ on plane $n$ is:

$${\tilde{u}}_2\sim K_{2}R(I_{3\times 3}+t{n}^T)K_1^{-1}{\tilde{u}}_1=H{\tilde{u}}_1$$

(2)

where ${\tilde{u}}_1$ and ${\tilde{u}}_2$ are the homogeneous coordinates of $u_1$ and $u_2$ respectively; The intersection point $p$ satisfies $n^{T}p+1=0$. $\sim$ means that ${\tilde{u}}_2$ is parallel to $H{\tilde{u}}_1$.

If camera parameters $K_1,K_2$, $R$ and $t$ are known, then we can deduce the following (3) from (2)

$$n^T{\tilde{y}}_1=-\frac{(R^T{\tilde{y}}_2\times t{)}^T(R^T{\tilde{y}}_2\times {\tilde{y}}_1)}{(R^T{\tilde{y}}_2\times t{)}^T(R^T{\tilde{y}}_2\times t)}=b$$

(3)

where ${\tilde{y}}_1=K_1^{-1}{\tilde{u}}_1$ and ${\tilde{y}}_2=K_2^{-1}{\tilde{u}}_2$ are the normalized coordinates of ${\tilde{u}}_1$ and ${\tilde{u}}_2$.

We use (3) of all matched points to construct a overdetermined equation and obtain the fitted global plane $n$ by solving this equation. Since the optical flow-based stitching method will be used to further align the images, the RANSAC method is not used here to calculate the plane with the most inliers. Instead, the global plane that fits all feature points as closely as possible is selected, misalignment caused by global plane projection will be better sloved during optical flow blending.

3.1.3. Projection Matrix Adjustment and Sliding Camera Generation

To preserve the viewpoint in the non-blending area of $I_2$, it is only required to satisfy $I\left(\tilde{u}\right)=I_2\left(N\tilde{u}\right)=I_2\left({\tilde{u}}_2\right)$. $N$ is a similarity transformation between $I_2$ and $I$, and can be obtained by fitting the matched feature points:

$$N_{*}=\underset{S}{min}\sum _{j=1}^n\left|\right|N{\tilde{u}}_1^j-{\tilde{u}}_2^j|{|}^2$$

(4)

where ${\tilde{u}}_1^j$ and ${\tilde{u}}_2^j$ are the homogeneous coordinates of pixels in $I_1$ and $I_2$ respectively.

Therefore, in the non-overlapping area of $I_2$, $\tilde{u}=N_{*}^{-1}{\tilde{u}}_2=N_{*}^{-1}{K}_2{R}_2(S-t)$, where $S$ is the corresponding 3D point of ${\tilde{u}}_2$ on plane $n$. So we get the projection matrix $P_2^{\prime }$=$N_{*}^{-1}{K}_2{R}_2$.

By RQ decomposition, the internal parameter matrix $K_2^{\prime }$ and rotation $R_2^{\prime }$ are extracted from $P_2^{\prime }$:

$$\begin{array}{cc}\hfill N_{*}^{-1}{K}_2{R}_2& =K_2\left(K_2^{-1}{N}_{*}^{-1}{K}_2\right)R_2\hfill \\ & =K_2\left(K_{*}{R}_{*}\right)R_2=K_2^{\prime }{R}_2^{\prime }\hfill \end{array}$$

(5)

where $K_2^{\prime }$ and $R_2^{\prime }$ are upper triangular matrix and rotation matrix respectively; and the third line of both matrices is $\left(\begin{array}{ccc}0& 0& 1\end{array}\right)$.

Compared with $P_2$, $P_2^{\prime }$ has a different intrinsic parameter matrix, and its rotation matrix only differs by one rotation around $Z$ axis, and its optical center $t$ is not changed.

From $I_1$ to $I_2$, the projection matrix transitions from $P_1$ to $P_2^{\prime }$. $K_m$, $R_m$ and $t_m$ in the projection matrix $P_m=K_m\left[R_m\right|-t_m]$ of the virtual camera can be calculated by weighting the corresponding intrinsic parameter matrices, rotation matrices and translation vectors of $I_1$ and $I_2$:

$$t_m=(1-m)0_{3\times 3}+mt$$

(6)

$$K_m=(1-m)K_1+m{K}_2^{\prime }$$

(7)

$$q_m=\frac{\mathit{sin}((1-m)\theta )}{\mathit{sin}(\theta )}{q}_1+\frac{\mathit{sin}(m\theta )}{\mathit{sin}(\theta )}{q}_2$$

(8)

where $q_1,q_2,q_m$ represent the quaternions corresponding to $I_{3\times 3}$, $R_2^{\prime }$ and $R_m$, $\theta$ is the angle between $q_1$ and $q_2$, and $m$ is the weighting coefficient.

The weighting coefficient $m$ can be calculated by the method in AANAP [16]:

$$m=⟨\overrightarrow{k_m{P}_{*}},\overrightarrow{k_m{K}_M}⟩/|\overrightarrow{k_m{K}_M}{|}^2$$

(9)

The meaning of each variable is shown in Figure 3. In the overlapping area, if $u$ corresponds to sliding camera $(K_m,R_m,t_m)$, then the relation between $u$ and $u_i$ in $I_i(\mathrm{i}=\mathrm{1,2})$ can be expressed as:

$$\tilde{u}=K_m{R}_m(I+t_m{n}^T/d)K_1^{-1}{\tilde{u}}_1=H_m^1{\tilde{u}}_1$$

(10)

$$\tilde{u}=H_m^1{H}^{-1}{\tilde{u}}_2=H_m^2{\tilde{u}}_2$$

(11)

(10) and (11)are also applicable to non-overlapping area. Projecting $I_1$ and $I_2$ through $H_m^1$ and $H_m^2$ onto the mosaic respectively to get warped images $I_1^{\prime }$ and $I_2^{\prime }$. Obviously:

$$I_i^{\prime }\left(\tilde{u}\right)=I_i\left(\right(H_m^i{)}^{-1}\tilde{u}\left)\right(i=\mathrm{1,2})$$

(12)

Figure 4 shows the experiment result on two school images used in [10]. Due to the parallax between $I_1$ and $I_2$, blending $I_1^{\prime }$ and $I_2^{\prime }$ will cause ghosting. Therefore, the next section will use an optical flow based blending method (AOF) to further align the images.

3.2. SC: Viewpoint Preservation Based on Sliding Camera

The mosaic generated by the SC method will inevitably has misalignment in most cases. So the optical-flow-based method is further employed to achieve more accurate alignment. This section firstly introduces the image alignment process based on asymmetric optical flow (AOF), and then details the calculation method of AOF.

3.2.1. Image Fusion Process of AOF

$I_1$ and $I_2$ are projected onto the custom projection surface to obtain warped images $I_1^{\prime }$ and $I_2^{\prime }$, which are blended to generate the mosaic $I$ . As the 3D points of the scene are not always on the projection plane, ghosting artifacts can be seen in the mosaic, as shown in Figure 4 in the above section. Direct multi-band image blending will lead to artifacts and blurring. As shown in Figure 5, point $P$ is projected to two points $p_1$ and $p_2$ in the mosaic, resulting in duplication of content. To solve the ghosting problem in the mosaic, the optical flow-based blending method in [22] is adopted.

Suppose $F_{2\to 1}\left(p_2\right)$ represents the optical flow value of $p_2$ in $I_2^{\prime }$ and $F_{1\to 2}\left(p_1\right)$ represents the optical flow value of $p_1$ in $I_1^{\prime }$. If the blending weight of pixel $\tilde{p}$ in the overlapping area is $\lambda$ (from the non-overlapping area of $I_1^{\mathrm{\text{'}}}$ to the non-overlapping area of $I_2^{\mathrm{\text{'}}}$, $\lambda$ gradually transitions from 0 to 1, as shown in Figure 5, then after blending, the pixel value of image $I$ at $\tilde{p}$ is:

$$I\left(\tilde{p}\right)=\left(1-\lambda \right)I_1^{\prime }\left({\tilde{p}}_1\right)+\lambda I_2^{\prime }\left({\tilde{p}}_2\right)$$

(13)

Where ${\tilde{p}}_1=\tilde{p}+\lambda F_{2\to 1}\left(\tilde{p}\right)$ represents the corresponding value of $\tilde{p}$ in $I_1^{\prime }$,${\tilde{p}}_2=\tilde{p}+\left(1-\lambda \right)F_{1\to 2}\left(\tilde{p}\right)$ represents the corresponding value of $\tilde{p}$ in $I_2^{\prime }$. That is, for any pixel $\tilde{p}$ in the overlapping area of the mosaic, its final pixel value can be obtained by a weighted combination of its corresponding values in the two warped images using optical flow.

In order to get a better blending effect, following the method presented by Meng and Liu [23], a softmax function is used to make the mosaic transition quickly from $I_1^{\mathrm{\text{'}}}$ to $I_2^{\mathrm{\text{'}}}$, narrowing the blending area. And if the optical flow value of a warped image is larger, the salience is higher, and the blending weight of the warped image should be increased accordingly. Therefore, the following blending weight $\beta$ can be used as the blending weight:

$$\beta =\frac{\mathit{exp}({\alpha }_s\lambda (1+{\alpha }_m{M}_2))}{\mathit{exp}({\alpha }_s(1-\lambda )(1+{\alpha }_m{M}_1))+\mathit{exp}({\alpha }_s\lambda (1+{\alpha }_m{M}_2\left)\right)}$$

(14)

where $M_1=\left|\right|F_{2\to 1}\left(\tilde{p}\right)\left|\right|$ and $M_2=\left|\right|F_{1\to 2}\left(\tilde{p}\right)\left|\right|$ represents the optical flow magnitude; ${\alpha }_s$ is the shape coefficient of the softmax function; and ${\alpha }_m$ denotes the enhancement coefficient of the optical flow. The larger ${\alpha }_s$ and ${\alpha }_m$ are, the closer $\beta$ is to 0 or 1, and the smaller the image transition area.

Also, similar to multi-band blending, a wider blending area is used in smooth and color-consistent areas, and a narrower blending area is used in color-inconsistent area. And the pixel consistency is measure using $D_c$:

$$D_c=\left|\right|I_1\left({\tilde{p}}_1\right)-I_2\left({\tilde{p}}_2\right)\left|\right|$$

(15)

The final blending parameter $\alpha$ is obtained：

$${\lambda }_d=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(c_d{D}_c\right)$$

(16)

$$\alpha =(1-{\lambda }_d)\lambda +{\lambda }_d\beta$$

(17)

$\beta$ corresponds to a fast transition from $I_1^{\prime }$ to $I_2^{\prime }$, $\lambda$ corresponds to a linear transition from $I_1^{\prime }$ to $I_2^{\prime }$. When the color differs smally, he transition from $I_1^{\prime }$ to $I_2^{\prime }$ is linear, and when the color difference is large, we tend to have a fast transition from $I_1^{\prime }$ to $I_2^{\prime }$.

Then the pixel value of the mosaic is:

$$I\left(\tilde{p}\right)=(1-\alpha )I_1\left({\tilde{p}}_1\right)+\alpha I_2\left({\tilde{p}}_2\right)$$

(18)

The curve in the left panel of Figure 6 shows the curve of $\beta$ changing with $\lambda$ under different optical flow intensities. $\beta$ can be used to achieve quick transition of the mosaic from $I_1^{\mathrm{\text{'}}}$ to $I_2^{\mathrm{\text{'}}}$, narrowing the transition area. In case of large optical flow, the blending weight of the corresponding image can be increased to reduce the transition area. The curve in the right panel of Figure 6 shows the impact of ${\lambda }_d$ on the curve of $\alpha$ as a function of $\lambda$. When ${\lambda }_d$ is small, a wider fusion area tends to be used; otherwise, a narrower fusion area is used, which is similar to the blending of different frequency bands in multi-band blending method.

3.2.2. Calculation of Asymmetric Optical Flow

The general pipeline of the optical flow calculation is to construct an image pyramid, calculate the optical flow of each layer from coarse to fine, and use the estimated current-layer optical flow divided by the scaling factor as the initial optical flow of the finer layer until the optical flow of the finest layer is obtained [23,24,25]. Different methods are proposed to achieve better solutions that satisfy brightness constancy assumptions, solve large displacements and appearance variation [27,28], address edge blur and improve temporal consistency [29,30,31]. Recently, some deep learning methods are proposed. For example, RATF(recurrent all-pairs field transforms for optical flow) [32] extracts per-pixel features, builds multi-scale 4D correlation volumes for all pairs of pixels, and iteratively updates a flow field through a recurrent unit. FlowFormer(optical flow Transformer) [33] based on transformer neural network architecture designs a novel encoder which effectively aggregates cost information of correlation volume into compact latent cost tokens, and a recurrent cost decoder which recurrently decodes cost features to iteratively refine the estimated optical flows.

In order to improve the optical flow calculation speed, we use the method based on optical flow propagation and gradient descent adopted in Facebook surround360 [34] to calculate the optical flow. When calculating the optical flow of each layer, first calculate the optical flow of each pixel from top to bottom and from left to right. From the optical flow values of the current-layer left and top pixels and upper-layer same-position pixel, the value with minimum error represented by (19) is selected as the initial value of the current pixel. Then a gradient descent method is performed to update the optical flow value of the current pixel, and then spread to the right and bottom pixels, as a candidate for the initial optical flow of the right and bottom pixels. After completing the forward optical flow propagation from top to bottom and from left to right, perform a reverse optical flow propagation and gradient descent from bottom to top and from right to left to obtain the final optical flow value.

When calculating the optical flow value $F\left(u\right)$ of pixel $u$, the error function $E\left(F\right(u\left)\right)$ used is:

$$E\left(F\right(u\left)\right)=E_I+{\alpha }_S{E}_S+{\alpha }_T{E}_T$$

(19)

$$E_I\left(F\right(u\left)\right)=\left|\right|\nabla I_1\left(u\right)-\nabla I_2(u+F(u\left)\right)\left|\right|$$

(20)

$$E_S\left(F\right(u\left)\right)=\left|\right|F\left(u\right)-G(u;\sigma )*F\left(u\right)\left|\right|$$

(21)

$$E_T\left(F\right(u\left)\right)=\left|\right|F\left(u\right)*D(1/W,1/H)\left|\right|$$

(22)

where, $E_I$ denotes the optical flow alignment error of the edge image (which is Gaussian filtered to improve the robustness), $E_S$ denotes the consistency error of the optical flow, $G(u;\sigma )*F\left(u\right)$ denotes the Gaussian-filtered optical flow of pixel $u$; $E_T$ denotes the magnitude error after normalization of optical flow, with excessively large optical flow being penalized; $W$ and $H$ are the width and height of current-layer image, respectively; and $D(1/W,1/H)$ denotes the diagonal matrix with diagonal elements $1/W$ and $1/H$.

3.3. Estimation of Image Intrinsic and Extrinsic Parameters

The SC-AOF method requires a known camera parameters of images $I_1$ and $I_2$. When only the intrinsic parameters $K_1$ and $K_2$ of an image are known, the essential matrix ${\left[t\right]}_{\times }R$ between two images can be obtained by feature point matching, and the rotation matrix $R$ and translation vector $t$ between images can be obtained by decomposing the essential matrix. When both intrinsic and extrinsic parameters are unknown, the intrinsic parameters can be estimated by calibration [35,36] firstly, and then the extrinsic parameters of the image can be estimated accordingly. In these cases, both intrinsic and extrinsic parameters of image $I_1$ and $I_2$ can be estimated robustly.

When none of the above methods is feasible, it is necessary to calculate the fundamental matrix from the matched feature points, and restore the camera internal and external parameters.

When camera has zero skew, known principal point and aspect ratio, then each intrinsic parameter matrix has only 1 degree of freedom (focal length of camera). The total degree freedom of camera parameters is 7($t$ has 2 degrees of freedom because the scale cannot be recovered, $R$ has three degrees of freedom, and each camera has 1 degree of freedom), which is equal to the fundamental matrix $F$'s degree of freedom. The internal and external parameters of the image can be recovered using a self-calibration method [37]. But even if these constraints are met, The solved camera parameters by [37] suffers from large errors when the scene is approximately planar and the matching error is large. Therefore, here we use the method of optimizing the objective function in [6] to solve the internal and external parameters of the camera.

In order to get an accurate fundamental matrix, firstly, the feature points need to be distributed more evenly in the image. As shown in Figure 7, a uniform and sparse distribution of feature points can both reduce the computation time and obtain more robust intrinsic and extrinsic camera parameters and global projection plane, which will lead to improved stitching results.

Secondly, it is necessary to filter the matched feature points to exclude the influence of outliers. Use the similarity transformation to normalize the matched feature points. After normalization, the mean value of the feature points is 0, and the average distance to the origin is $\sqrt{2}$.

Thirdly, multiple homographies are estimated to exclude outlier points. Let $F_{cond}$ and $n$ denote all matched feature points and the total number of matched feature points. In $F_{cond}$, The RANSAC method with threshold $\eta =0.01$ is applied to compute homography $H^1$ and its inlier set $F_{inlier}^1$, and the matches of isolated feature points which have no neighboring points within 50 pixel distance is removed from $F_{inlier}^1$. A new candidate set is generated by removing $F_{inlier}^1$ from $F_{cond}$. Repeat the above steps to calculate $m$ homography matrices $H^m$ and corresponding inlier set $F_{inlier}^m$ until $\left|\right|F_{inlier}^m\left|\right|$<20 or $\left|F_{inlier}^m\right||<0.05n$. The final inlier set is $F_{inlier}={\cup }_{i=1}^m{F}_{inlier}^i$. If $m=1$, then there is only one valid plane. In this case, apply the RANSAC method with threshold $\eta =0.1$ to recalculate homography $H^1$ and the corresponding inlier set $F_{inlier}^1$.

After excluding the outliers, for any matched points $\{x_1,x_2\}$ in the inlier set $F_{inlier}$, the cost function is:

$$E(x_1,x_2)=\left(1-\lambda \right)h(r_e,{\sigma }_e)+\lambda h(r_p,{\sigma }_p)$$

(23)

where the $\lambda$ balanced epipolar constraint and the infinite homography constraint, and generally take $\lambda =0.01$. $h$ is a robust kernel function, which can mitigate the effect of mis-matched points on the optimization of camera internal and external parameters. $r_e$ and $r_p$ denote the projection errors of the epipolar constraint and of the infinite homography constraint, respectively:

$$r_e=\frac{1}{\rho }{x}_2^T{K}_2^{-T}[t{]}_{\times }R{K}_1^{-1}{x}_1=\frac{1}{\rho }{x}_2^{T}F{x}_1$$

(24)

$$r_p=x_2-\frac{1}{\omega }{K}_{2}R{K}_1^{-1}{x}_1$$

(25)

where $\rho$ denotes the length of the vector composed of the first two components of $F{x}_1$. That is, assuming $F{x}_1=(a,b,c{)}^T$, then $\rho =\sqrt{a^2+b^2}$. $w$ represents the third component value of the vector $K_{2}R{K}_1^{-1}{x}_1$.

4. Experiment

To verify the effectiveness of the SC-AOF method, The mosaics generated by our method and the existing APAP [4], AANAP [16], SPHP [14], TFT [7], REW [10] and SPW [18] are compared on some typical datasets used by other methods to verify the feasibility and advantages of the SC-AOF method in solving deformation and improving alignment; Next, SC-AOF method is used together with other method to demonstrate its good compatibility.

4.1. Effectiveness Analysis of SC-AOF Method

In this section, various methods of image stitching are compared and analyzed based on three indicators: perspective deformation, local alignment and running speed. The experimental setup is as follows.

• The first two experiments compare typical methods for solving perspective deformation and local alignment respectively, and all the methods in the first two experiments are included in the third experiment to show the superiority of SC-AOF method in all aspects.
• Since the averaging methods generally underperforms linear blending ones, all methods to be compared adopt linear blending to achieve the best performance.
• All methods other than ours use the parameters recommended by their proposers. Our SC-AOF method has the following parameter settings in optical flow-based image blending: ${\alpha }_s=$10, ${\alpha }_m=$100, and $c_d=$10

4.1.1. Perspective Deformation Reduction

Figure 8 shows the results of the SC-AOF method versus SPHP, APAP, AANAP and SPW methods for perspective deformation reduction in image stitching. School, building and park square datasets were used in this experiment. We can see from Figure 8 that, compared with the other methods, our SC-AOF method makes the viewpoint of the stitched image change in a more natural manner and effectively eliminates perspective deformation. As explained below, all other methods underperform our SC-AOF method.

The image stitched using the APAP method has its edges stretched to a large extent. This is because it does not process perspective deformation. This method only serves as a reference to verify the effectiveness of perspective deformation reducing algorithms.

The AANAP algorithm can achieve a smooth transition between the two viewpoints, but results in severely "curved edges". And there is even more severe edge stretching for the park square dataset than that of the APAP method. This is because, when the AANAP method extrapolates from homographies, it linearizes the homography in addition to similarity transformation, causing affine deformation in the final transformation.

Compared with the APAP method, the SPW method makes no significant improvement in perspective deformation, except for the image in the first row. SPW preserves perspective consistency, so multiple-viewpoint method excels in solving perspective deformation compared to single-viewpoint method.

The SPHP algorithm performs well overall. However, it causes severe distortions in some areas (red circled in Figure 8c) due to the rapid change of viewpoints. This is because the SPHP method estimates the similarity transformation and interpolated homographies from global homography. As a result, the similarity transformation can't not reflect the real scene information and the interpolated homographies may deviate from a reasonable image projection.

4.1.2. Local Alignment

Figure 9 shows the results of the SC-AOF method versus APAP, TFT and REW methods for local alignment in image stitching. It can be seen that SC-AOF performs well in all scenes, showing the effectiveness of our method in local alignment.

• The APAP method performs fairly well in most images, though with some alignment errors. This is because the moving DLT method smooths the mosaics to some extent.
• The TFT generated stitched image of excellent quality in planar areas. But when there is a sudden depth change in the scene, there will be serious distortions. This is because large errors appear when calculating plane using three vertices of a triangle in the area with sudden depth changes.
• The REW method has large alignment error in the planar area and aligns the images better than the APAP and TFT method in all other scenes. This is because the fewer feature points in the planar area might be filtered out as mismatched points by the REW method.

The SSIM (structural similarity) is employed to objectively describe the alignment accuracy of different methods. The scores of all methods on the datasets temple, rail tracks, garden, building, school, park-square, wall, and cabinet are listed in Table 1.

• APAP and AANAP have high scores on all image pairs, but the scores are lower than our method and REW, proving that APAP and AANAP blur mosaics to some extent.
• When SPHP is not combined with APAP, only the global homography is used to align the images, resulting in lower scores compared to other methods.
• TFT has higher scores on the datasets except for the building dataset. TFT can improve alignment accuracy but also bring instability.
• SPW combines quasi-homography and content-preserving warping to align images, which adding other constraints while also reducing the accuracy of alignment, resulting in lower scores compared to REW and our method.
• Both REW and our method use global homography matrix to coarsely align the images. Afterwards, in REW and our method, a deformation field and optical flow are applied to further align the images respectively. Therefore, the both methods has higher scores and robustness than other methods.

4.1.3. Stitching Speed Comparison

The running speed is a direct reflection of the efficiency of each stitching method. Table 2 shows the speed of the SC-AOF method versus the APAP, AANAP, SPHP, TFT, REW and SPW methods. The temple, rail tracks, garden, building, school, park-square, wall, and cabinet datasets contributed by [4,10] are used in this experiment. It can be seen that the REW algorithm has the fastest stitching speed. The reason is that it only needs to calculate TPS parameters based on feature point matching and then compute the transformations of grid points quickly. Our SC-AOF method ranks second in terms of stitching speed, and the AANAP algorithm requires the longest running time. Both the APAP and AANAP methods calculate the local homographies based on moving DLT, and the AANAP method needs to calculate the Taylor expansion of anchor points. Combined with the results in Figure 8 and Figure 9, it is verified that the SC-AOF method has the best overall performance, which can maintain a desirable operation efficiency while guaranteeing the final image quality. Our method can have broad applications and promotion significance.

4.2. Compatibility of SC-AOF Method

The SC-AOF method can not only be used independently to generate stitched image with reduced perspective deformation and low alignment error, but also be decomposed (into SC method and image blending method) and combined with other methods to improve the quality of the mosaic.

4.2.1. SC Module Compatibility Analysis

The sliding camera (SC) module in the SC-AOF method can not only be used in the global alignment model, but also be combined with other local alignment models (e.g., APAP and TFT) to solve perspective deformation while maintaining the alignment accuracy. The implementation steps are as follows.

• Use the global similarity transformation to project $I_2$ onto the $I_1$ coordinate system to calculate the size and mesh vertices of the mosaic;
• Use (6)-(9) to calculate the weights of mesh vertices and the projection matrix, replace the homography $H$ in (2) with the homography matrix in local alignment model, and bring them into (12) to compute the warped images and blend them.

Figure 10 presents the stitched images when using the TFT algorithm alone vs. using the TFT algorithm combined with the SC method. It can be found that the combined method is more effective in mitigating the edge stretching, and that it generates more natural images. This shows that the SC method can effectively solve perspective deformation suffered by the local alignment method.

4.2.2. Blending Module Compatibility Analysis

The asymmetric optical flow-based blending in the SC-AOF method can also be used in other methods to enhance the final stitching effect. The implementation steps are as follows.

• Generate two projected images using one of the other algorithms and calculate the blending parameters based on the overlapping areas;
• Set the optical flow value to be 0, replace linear blending parameter $\lambda$ with $\alpha$ in (17) to blend warped images, preserve the blending band width in the low frequency-area and narrow the blending width in the high-frequency area to get a better image stitching effect.
  Figure 11 shows the image stitching effect of the APAP algorithm when using linear blending vs. when using our bending method. It can be seen that the blurring and ghosting in the stitched images are effectively mitigated when using our blending method. This shows that our blending algorithm can blend the aligned images better.

5. Conclusions

In this paper, to solve the perspective deformation and misalignment in image stitching using homographies, a SC-AOF method is proposed. In image warping, a new virtual camera and a projection matrix are generated as the observation perspective in the overlapping area by interpolating between two projection matrices; and the overlapping area transitions gradually from one viewpoint to another to achieve viewpoint preserving and smooth transition of the stitched image and thus solve the perspective deformation problem. In image blending, optical flow-based blending algorithm is proposed to further improve alignment accuracy. The width of blending area is automatically adjusted according to softmax function and alignment accuracy. Finally, extensive comparison experiments are conducted to demonstrate the effectiveness of our algorithm in reducing perspective deformation and improving alignment accuracy. In addition, our algorithm had broad applicability, as its component modules can be used with other algorithms to mitigate edge stretching and improve alignment accuracy.

However, the proposed local alignment method would fail if the input images contain large parallax, which will cause severe occlusion to prevent us from obtaining the correct optical flow. The problem of local alignment failure caused by large parallax also exists in other local alignment methods. Exploring more robust optical flow calculation and occlusion processing method to reduce misalignment in large parallax scene is an interesting research direction of the future work.