A novel algorithm for optimizing bundle adjustment in image sequence alignment

The L-M algorithm is a widely used algorithm for solving non-linear least squares problems, which combines the concepts of gradient descent and the Gauss-Newton algorithm to provide a robust approach for parameter estimation.

Let $\varphi(x)=(\varphi_{1}(x),\varphi_{2}(x),...,\varphi_{m}(x))^{T}$, the objective is to minimize the sum of squared residuals

The update rule for the parameters $x$ is given by

where $J(x_{k})$ is the Jacobian matrix of partial derivatives, $I$ is the identity matrix and $\mu_{k}$ is referred to as damping parameter, which controls the behavior of the algorithm.

A larger $\mu_{k}$ makes the update resemble gradient descent while a smaller $\mu_{k}$ approximates the Gauss-Newton method. This adaptability allows the L-M algorithm to efficiently handle a wide range of problem conditions.

In this approach, we employ the trust-region method to adjust $\mu_{k}$. First, we define a quadratic function at the current iteration point to represent the predicted value

The step $d_{k}$ is then computed based on the current value of $\mu_{k}$, and the ratio of the actual decrease to the predicted decrease is evaluated as

Equation (4) reflects the degree to which the linearized model aligns with the nonlinear objective function.

• •

  When $\xi_{k}$ approaches 1, the quadratic model $P(d)$ provides a good approximation of the objective function at $x_{k}$, suggesting that $\mu_{k}$ should be reduced.

• •

  When $\xi_{k}$ approaches 0, the quadratic model $P(d)$ poorly approximates the objective function at $x_{k}$, indicating that $\mu_{k}$ should be increased.

• •

  When $\xi_{k}$ is neither close to 0 nor 1, $\mu_{k}$ is considered to be appropriately chosen, and no adjustment is necessary.

Typically, the threshold values for $\xi_{k}$ are set to 0.25 and 0.75. Based on this, the update rule for $\mu_{k}$ in the L-M method is as follows

In summary, the complete flow of the L-M algorithm is as follows

An essential component of BA is the camera projection model, which maps the 3D coordinates of objects in the real world onto 2D coordinates in an image. In this paper, we focus on the application of a projection model specifically for cryo-ET image sequence[10]

where $(X,Y,Z)^{T}$ is the 3D coordinates of the fiducial markers, $(u,v)^{T}$ is the 2D projection of the fiducial markers in the image, $s$ is the image scale change, $\textbf{t}=(t_{0},t_{1})^{T}$ is the translation during projection and P represents the orthogonal projection matrix. R denotes the rotation matrix of the projection process, the details of R and P as follows

$\textbf{x}_{c}=(s,\alpha,\beta,\gamma,t_{0},t_{1})^{T}$ is called as camera parameters and $\textbf{x}_{p}=(X,Y,Z)^{T}$ is called as point parameters.

When sufficient projection data for a single 3D point is available, triangulation can be used to approximate the coordinates of that point. However, due to various uncertainties, the initial estimates of the 3D coordinates may not be entirely accurate. When this estimated 3D point is reprojected onto the 2D image, a discrepancy, or error, arises in comparison to the observed point.

Let the observed point be denoted as $\textbf{z}\stackrel{{\scriptstyle\Delta}}{{=}}(u,v)^{T}$ and the point obtained by reprojection can be expressed as follows

where $(\cdot)_{0}$ represents the initial estimation of parameters. Thus the reprojection error can be written as

To account for all observations, we can introduce a subscript to represent the error for each observation. Let $\textbf{z}_{ij}$ denotes the data generated by the $i$-th marker in the $j$-th image. The overall cost function can be expressed as

During the sample collection process, some marker points may become obscured at certain angles and thus will not be visible. To represent the visibility of these points, we introduce a visibility indicator $\delta_{ij}=\{0,1\}$. If the projection of the $j$-th marker is visible in the $i$-th image, $\delta_{ij}=1$; otherwise, $\delta_{ij}=0$.

Solving this least-squares problem involves jointly optimizing both the camera parameters $\textbf{x}_{c}$ and the point parameters $\textbf{x}_{p}$, a process known as Bundle Adjustment.

The effectiveness of the adjustments is evaluated using the following average residual formula

where $||\cdot||_{1}$ denotes the 1-norm, $x_{c}^{*}$ and $x_{p}^{*}$ are the adjusted parameters. A smaller $L_{1}$ value indicates a lower reprojection error.

Suppose $f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a nonlinear function with a continuous second derivative. Consider the following optimization problem

Generally it can be solved by the following iteration

where $d_{k}$ is step size.

If we interpret $d_{k}$ as a control input in a control system, the iteration (12) exhibits a structure analogous to that of the following discrete-time linear time-invariant system

By designing a performance index related to the cost function [4], the optimization problem (11) can be reformulated as the following optimal control problem,

where $f(x_{N+1})$ represents the terminal cost function, $N$ is the terminal time, and $x_{0}$ is given. The positive definite matrix $R$ denotes the control weight.

Based on the Pontryagin Maximum Principle, the optimal trajectory satisfies a set of forward and backward equations. Solving these equations yields the optimal controller for the problem (14)

Thus, we can derive the following iterative algorithm

Although the above iteration is implicit, we can obtain an approximate explicit solution through Taylor expansion. First, take the first-order Taylor expansion of $\nabla f(x_{i})$ at $x_{i-1}$

where $\nabla f_{k}=\nabla f(x_{k}),H_{k}=H(x_{k})$, and $H(x_{k})$ is the Hessian matrix of the objective function $f(x)$. By replacing all $x_{i}$ with $x_{k}$ for $i>k$, we obtain the following optimal control algorithm to address BA

The optimal control algorithm has been demonstrated to exhibit a superlinear convergence property[11], and the algorithm operates through the following stages

The selection of an appropriate weight matrix $R$ is critical for the effective operation of the OCA. For convenience, $R$ is generally defined as $R=\lambda\cdot I$, where $\lambda>0$ is a constant. Under the condition that $(R+H(x_{0}))>0$, smaller values of $\lambda$ lead to faster convergence of the OCA. In the simulation experiments presented in Section 4.2, a set of poorly initialized data was generated, with initial points positioned far from the extremum points. In such cases, a larger $\lambda$ is typically required to ensure algorithmic convergence; however, this also results in a slower convergence rate for the OCA. To resolve this issue, we propose an update procedure based on the bisection method, which allows the weight matrix $R$ to progressively decrease as the number of iterations increases. This approach ensures algorithmic convergence while simultaneously enhancing convergence speed. The detailed update strategy is presented as follows.

It is important to note that this update procedure does not affect the convergence of OCA. Given that $\lambda_{k}$ is non-increasing, it follows that the spectral radius $\rho((R_{k}+H(x_{k}))^{-1}R_{k})$ is also non-increasing. As $\lambda_{k}$ decreases, the convergence rate of OCA improves. For a detailed proof, the reader is referred to [11].

We selected three real-world cryo-electron microscopy datasets, Centriole^††http://bio3d.colorado.edu/imod/files/tutorialData-1K.tar.gz, VEEV^††https://doi.org/10.5281/zenodo.11172321 and Vibrio ^††https://doi.org/10.5281/zenodo.11172858 to evaluate our algorithms. The Centriole dataset consists of a tilted sequence of 64 projections, with the projected images tilted between -61° and +65° at 2° intervals. Each projection is a 1024x1024-pixel image, with each pixel corresponding to 1.01 nm. The VEEV dataset consists of a tilted sequence of 21 projections, with the projected images tilted between -50° and +50° at 5° intervals. Each projection is a 1536x2048-pixel image, with each pixel corresponding to 0.2 nm. The Vibrio dataset[13] consists of a tilted sequence of 41 projections. The tomographic reconstructions of these samples were performed using the incline program from the MarkerAuto software suite[3].

Table 1 presents the experimental results of the real datasets, indicating that the Optimal Control algorithm has a faster convergence speed compared to the L-M algorithm.

Using the VEEV dataset as an illustrative example, Figure 1 presents the error in the parameters to be optimized at each iteration, represented as $x_{k}-x_{k-1}$. Notably, when $k=28$, the OCA satisfied the convergence criterion of $\varepsilon<10^{-6}$. In contrast, the L-M algorithm exhibited oscillatory behavior, which led to an increased number of iterations. Consequently, it only met the convergence criterion at $k=236$.

To obtain a broader range of experimental results, we randomly constructed $n$ 3D points and $m$ projection images at various angles for simulated dataset experiments, where $n=20,40,60$ and $m=21,41,64$. The projection images are sized $1024\times 1024$ pixels, and the 3D points are distributed within an $800\times 800\times 400$ space. To enhance the realism of the simulated dataset, we processed the constructed data as follows

• •

  Step 1 Randomly select 5%–30% of the 2D projection points as invisible.

• •

  Step 2 Add Gaussian noise with a standard deviation of $a\%\times 1024$ to the 2D projection points, where $a=0.2,2,10$.

• •

  Step 3 Add $b\%$ noise to the camera parameters, with noise magnitude based on the average of each parameter (avg), uniformly distributed in the range $[-\text{avg},\text{avg}]$, where $b=5,10$.

• •

  Step 4 Select 5% of the 2D projection points as outliers and add Gaussian noise with a mean of $0.01\times 1024$ and a standard deviation of $0.04\times 1024$.

For convenience, we denote the dataset with parameters $m=21,n=20,a=0.2,b=5$ as “21c_5%_20p_0.2%.” To ensure a fair comparison, the weight matrix $R$ in the OCA and the damping factor $\mu$ in the L-M algorithm were iteratively adjusted to promote faster convergence. The results of this experiment are presented in Table 2.

When large noise is introduced into the simulated dataset, the initial estimates of the camera parameters and the 2D point coordinates exhibit significant deviations. This suggests that the initial parameters are considerably distant from the extremum, necessitating a greater number of iterations for both the OCA and L-M algorithm. It is worth noting that, in the L-M algorithm, the damping factor $\mu$ is continuously updated, whereas the weight matrix $R$ in the OCA remains constant throughout the iterations. If we allow the matrix $R$ in the OCA to be updated gradually similar to $\mu$ in the L-M algorithm, the OCA could potentially achieve faster convergence. The detailed experimental results are displayed in the Table 3.