Least squares optimization is a standard optimization tool across many branches of science, due to the fact that Gauss’ approximation to the Hessian makes second-order Newtonian optimization computationally efficient. The use of least squares optimization in photometric costs dates back to at least 1981, when Lucas & Kanade applied it to image registration . Since then, inverse approaches have further improved the efficiency [2, 3] and robustness  of least squares methods on certain registration problems. Applications of photometric least squares optimization have also increased, covering such diverse tasks as visual odometry  and active appearance models .
Similarly, cross correlation, a standard statistical tool across many branches of science, is used to measure the similarity between two signals. The normalized version (NCC), sometimes called zero-mean normalized cross correlation, has been used in image registration as far back as 1972 
. This version normalizes the means and variances of the data before applying cross correlation, making the measure robust to changes in gain and bias. Common applications of NCC today include multi-view stereo[8, 9] and 2-d keypoint tracking , both of which are instances of local image registration.
These two methods are complementary; the first offers efficiency and the second provides robustness to changes in intensity. However, they have not yet been successfully combined. Indeed, of the methods that use NCC, current tracking systems tend to employ an exhaustive search , which does not scale well to higher dimensional search spaces, while popular multi-view stereo pipelines either optimize using BFGS [11, §6.1] with numerical gradients  (computed using finite differences), or a sampling approach , both of which are inefficient.
Gradient-based methods, which include least squares optimization, scale well with dimensionality by traversing the state space one direction (of high gradient) at a time. While previous works have optimized NCC costs using such methods, this work is the first to incorporate NCC into a least squares optimization, demonstrating that it is simple, efficient, and converges well. In addition, we introduce a sparse formulation that both handles local variations in intensity and improves robustness to occlusions.
The following subsections review existing gradient-based NCC optimizers, and methods of registration with local intensity variations, particulary in least squares frameworks. Section 2 then presents the NCC least squares framework, section 3 presents our results, and we conclude in section 4.
1.1 Gradient-based optimization of NCC
NCC has been optimized using gradient-based methods with analytic gradients, particularly for image registration: Irani & Anandan  use a second-order Newtonian framework, requiring the computation of second derivatives, which increases computational cost and amplifies image noise. Brooks & Arbel  extend inverse compositional optimization to arbitrary cost functions (including NCC), using the BFGS optimizer: a first order optimizer which builds an approximation to the Hessian over time. Evangelidis & Psarikis  define a least squares formulation of NCC (eqn. (13), employed here), but then maximize the standard formulation (eqn. (12)) as follows: they compute the Jacobian for the zero-mean data, then at each iteration compute a scalar “perturbation” that is applied to the reference data to account for normalization. The resulting step approximates the least squares one. However, we note that none of these approaches has been widely adopted.
1.2 Registration with local intensity variations
Registration methods can be broadly split into direct and feature-based approaches. Feature-based approaches, such as that presented in Brown & Lowe’s Autostitch , pre-compute features in each image, match them across images, then minimize a geometric error on the correspondences. This approach has the benefit of broad convergence, but is computationally expensive, and does not minimize the true measurement errors, which occur in the domain of pixel intensities, not feature locations. Rather than requiring correspondences, direct methods minimize an error on pixel intensities or some transformed pixel space (see below). Therefore, these methods exploit more image data, including edges, and have been shown to be more accurate than feature-based methods for direct odometry 16] have revolutionized 2-d tracking , but do not extend to higher-dimensional registration problems such as active appearance models  and direct odometry . These latter problems use optimization-based approaches, which are efficient but require a good initialization. As NCC is a direct error metric, and non-linear least squares is an optimization approach, we limit further discussion to direct, optimization-based methods.
Several categories of direct registration methods already exist to handle data with local intensity variations. We enumerate these with regard to image-based alignments, where scene effects such as non-Lambertian reflectance and a changing viewpoint, lighting effects (e.g. shadows moving with the time of day), and camera effects (e.g. lens vignetting) all create local changes in intensity.
The first category is image transformation. Prior to least squares optimization, images are converted pixelwise into a lighting invariant space. Some methods compute distance transforms over edge images [18, 19], converting intensity into a form of geometric, rather than photometric, error. Other methods compute a multi-channel descriptor per pixel, based on the local image texture, some of which are hand designed [20, 21, 22], others of which are learned . It should be noted that these descriptors are often not invariant to in-plane image rotations,111Invariance of a pixelwise transformation to an in-plane rotation requires , where is an image. so registration in such cases will fail.
The second category models intensity changes generatively. Lucas & Kanade  originally proposed a two-parameter gain and bias model, which can only model a global change in intensity. Silveira & Malis  extend the gain model to independent regions, but retain the global bias. This allows for more local intensity variations, at the expense of a higher-dimensional search space.
The final category consists of lighting invariant scores between registered patches, such as NCC and mutual information [25, 26]. Mutual information supports a wider range of lighting changes than the linear-invariance of NCC, but nevertheless requires a globally consistent transformation of intensities. In addition, a least squares formulation does not exist.
2 NCC least squares
Our task is formulated as that of finding a warp matrix that defines the registration between two data samples, and , by minimizing an NCC least squares cost between the registered samples (symbols defined below):
We assume single channel data for brevity (least squares trivially generalizes to multiple channels ). Therefore, the data represents a scalar field, of dimensions. The data could be -d (audio), -d (image), -d (volumetric), or higher. In this work, we only apply the formulation to image data, such that .
represent coordinates within the scalar field of the source data. These coordinates are the sample points whose error is to be minimized during registration. The function , transforms the coordinates from the source to target frame, as follows:
where the warp matrix encodes the registration,222 is the identity matrix/transformation.
is the identity matrix/transformation.so is kept on the manifold of allowed registrations (see sec. 2.3.1), and applies any non-linearities present in the measurement process to each column of ; in the case of images, this is a projection onto the image plane, and correction for camera calibration and lens distortion. Our evaluation tasks require only a simple projection: . Finally, the warp function samples the target data at the coordinates
The function transforms a sampled vector to one of zero mean and unit length.333 is a vector of ones; similarly, is a vector of zeros. Though different from the standard one, this NCC formulation, previously given by Evangelidis & Psarikis , has an identical cost function shape (see ncc_formulation). Our first contribution is demonstrating that this change in formulation allows NCC to be straightforwardly and efficiently (see ncc_one_pass) incorporated into a least squares framework.
2.2 Sparse formulation
Our second contribution is to introduce the following sparse, robust, locally normalized formulation, inspired by the success of sparse features for direct odometry :
The cost function consists of a sum of local NCC costs, making this cost invariant to local (not just global) variations in intensity. Furthermore, each cost is robustified by a function444Suitable robustification functions must satisfy and , and increase monotonically for positive values . Note that robustifying a single cost term, such as equation (1), is redundant. Since increases monotonically, the minima occur at identical warps.
, that downweights large errors, such that costs close to converging have more influence than costs far from converged. The robustification of costs to outlier data, such as occluded regions, is a standard tool in vision literature[28, 27].
2.3.1 Warp update parameterization
To support inverse  and ESM  approaches (sec. 2.3.3), we update the warp via a compositional approach . Given a vector of the change in variables , which is computed at each iteration of the optimization (see below), the warp is updated as follows:
where converts an update vector into a warp matrix, such that the set of warps is a group . In particular, must be the identity matrix, and, due to the projection from 3-d to 2-d in the case of images, the generators of each dimension must have zero trace . We use the warp update parameterizations
in our experiments, where encodes 2-d translations, and encodes homographies [29, eqn. 85–87].
2.3.2 Iterative update
Any standard non-linear least squares optimizer (Gauss-Newton, Levenberg-Marquardt, ) can be used to optimize equations (1) & (5); here we use Gauss-Newton. The per iteration update for equation (5) (of which equation (1) is a special case) is computed thus:
where is an update vector of zeros. The scalar value and matrix are Triggs’ correction factors [27, eqn. 11], functions of which account for robustification. In our implementation, this update is repeated until , or the least squares cost fails to go below the minimum found for three consecutive iterations.
2.3.3 Jacobian computation
The following three approaches to Jacobian computation for compositional warp updates have been proposed in the literature.
Jacobians can also be computed in the source image, at the identity warp [3, 2], , such that they are constant. When Gauss-Newton is used, the pseudo inverse can also be precomputed, resulting in a much faster update. However, Triggs’ second order correction for robust kernels cannot be used with a precomputed pseudo inverse, so the optimization must use iteratively reweighted least squares  instead: which results in slower convergence . Experiments here use the second order correction and therefore recompute the pseudo inverse each iteration.
Efficient Second-order Minimization (ESM)
Taking the average of the above two Jacobians,
, provides a more accurate estimate of the Hessian, improving both the rate and speed (number of iterations) of convergence.
Only the formula for the Jacobian of is given herein (see ncc_derivative). Formulae for other Jacobians are available in prior work . We also note that modern auto-differentiation tools (employed in our implementation) can automatically compute analytic derivatives at runtime.
Our experiments validate the value of the two contributions of this work. The first is a new method to optimize the well-known NCC cost, so we evaluate its performance relative to other such optimizers. The second is a sparse, robust formulation which improves performance in situations with local intensity variations and partial occlusion, so we compare this approach against other methods invariant to local intensity changes, as well as on partly occluded image regions.
We run quantitative experiments on homography-based alignment. The reasonably high dimensionality of this problem’s state space (8-d) indicates performance on the kind of problems for which photometric least squares methods remain popular, such as active appearance models  and direct odometry .
3.1 Costs evaluated
We evaluate three costs in our experiments:
Dense, globally normalized cost
This cost is defined by equation (1), where is a dense grid of pixel coordinates sampled at the pixel corners555Sampling (under bilinear interpolation) at pixel corners means that a sample must shift at least half a pixel in any direction before the image intensity changes non-linearly. In contrast, the intensity at a pixel center is non-differentiable. within a region of interest.
Sparse, locally normalized cost
This cost is defined by equation (5), without robustification (, , ). Each set of coordinates has 8 sample points on a 24 grid: spaced 2 pixels apart, centered on and aligned perpendicular to features called edgelets. Regions of constant image gradient produce the same NCC cost. Therefore, a change in gradient, an image edge, is required for registration. Edgelets are extracted at sub-pixel local maxima of gradient magnitude, in the direction of maximum gradient, as shown in Figure 2. Blocks have 8 samples, as vector operations on modern CPUs work well with this number. The height of 2 gives them some sensitivity to the orientation of the edge, while the width of 4 allows them to converge from a reasonable distance. The set consists of the blocks on all edgelets extracted within a region of interest.
Sparse, robustified, locally normalized cost
To ensure that results are broadly applicable, experiments should be run on a range of image textures, with real (not synthetic) data, locally varying intensities across images of the same scene. We therefore created a new dataset of 7 sets of 4 images, shown in Figure 1, with real-life changes in lighting due to capture at different times of day. We have also computed ground truth homographies between each pair of images within each set, using a standard feature matching followed by RANSAC approach .
3.2.1 Test case generation
We randomly select 100 50
50 pixel regions from each image. All regions contain at least 50 edgelets, and are fully visible in the 3 other images of the same scene. We generate a perturbation for each region, by applying a random shift, drawn from a normal distribution, to each corner separately, then scaling the perturbation so that the mean shift per corner is 0,..,10 pixels, creating starting distances of increasing magnitude for each region. Each perturbation of each region is warped into each image in the set (including the source image), using the ground truth homography, from which the initial warp from source to target is then computed[32, alg. 4.2]. This creates a total of 10011447 = 123,200 test cases, 92,400 of them between different images and 30,800 between identical images.
We have made the images, ground truth homographies and test cases publicly available .
3.3 NCC optimization comparison
Of the three pre-existing gradient-based NCC optimization methods [12, 13, 14], we compare our non-linear least squares (NLS) method against the two most recent of these: the first-order BFGS method , and the second-order approximation to NCC optimization , which the authors call enhanced correlation coefficient maximization (ECCM). The Newton-based method  employs a full second-order differentiation of the NCC cost (not given in the paper), requiring a significantly more complex implementation, which is not publicly available. For this reason, we omit comparison to that method here. However, we note that our least squares update is the Gauss-Newton approximation to the Newton update of the standard NCC cost. Therefore, we would expect our method to provide similar accuracy with a simpler, faster implementation.
We run all three optimizers on each of the three costs described above, using the warp update of eq. (8) (homography fitting) with each of the three Jacobian values described in section 2.3.3, (FWD), (INV), and (ESM), over every test case. We extend the baseline methods to support ESM, as well as our two sparse costs.
3.3.1 Rate of convergence
The primary goal of registration optimization is to reach an alignment that matches the ground truth; we call this convergence. Since matching ground truth exactly is unlikely, we allow an error of up to 1 pixel in the position of each corner in the target image for an optimization to be considered converged. The proportion of tests which converged is plotted against the mean corner starting distance, as shown in Figure 3(top row) for the case that source and target images differ. It should first be noted that the rate of convergence is not 100% for a starting distance of 0, and is also different between methods. The ground truth homography for some test cases does not coincide with the minimum of the NCC cost, and not all optimizers reach the optimum equally well. As a result, when the starting position is inside the converged region, worse optimizers actually get a higher rate of convergence. In the case that source and target images are identical, the ground truth and NCC optimum do coincide, giving all optimizers 100% convergence for a starting distance of 0, as shown Figure 4.
The following trends are visible in both Figures 3(top row) and 4. In the dense case (a), NLS performs marginally better than ECCM and significantly better than BFGS; in the sparse cases (b) & (c), NLS performs significantly better than both ECCM and BFGS; ESM provides better convergence overall for both NLS and ECCM, though marginally in the dense case; for BFGS, the FWD scheme performs best; when the sparse cost is robustified (c), the rate of convergence improves at all starting distances.
3.3.2 Optimization time
All tested methods are implemented in the same langauge (). As a result, optimization time is a strong overall indicator of the relative performance of the different methods, if not absolute achievable performance. Figure 3(bottom row) shows the convergence time of each method, across different starting distances and cost functions. It is important to note that NLS methods are faster, often significantly, than their BFGS or ECCM counterparts. In this high dimensional space, BFGS requires many more function evaluations. Meanwhile, the pertubation factor computation of ECCM is not inexpensive, slowing, in particular, the INV variant significantly (beyond the ESM variant, which uses fewer iterations, for sparse costs). It should also be acknowledged that for NLS, ESM is faster than FWD and converges more often, but INV is considerably faster still. Consequently, there is a trade-off between speed and efficacy when choosing between INV and ESM.
3.3.3 Final cost
The goal of image alignment is to converge to ground truth. However, when comparing optimizers, one should also consider the lowest cost achieved. As we have seen, the two don’t always coincide. Therefore, we also show the relative cost achieved by each method, compared to the lowest cost found across all methods, in Figure 5. The results are similar to the relative convergence rates: for the dense case, ECCM performs slightly worse than NLS, and for the sparse cases, NLS methods perform significantly better than all methods except BFGS FWD, which is still worse (and significantly slower).
3.3.4 Qualitative analysis of costs & optimizers
The causes of the relative performance of the costs and optimizers tested are illuminated by the qualitative plots shown in Figure 3, which visualize the cost surfaces generated by each of the above costs, for a 5050 pixel region centered on the image shown in Fig. 2(a). Trajectories computed by the ESM variant of each optimizer tested are also shown, starting from 2-d shifted positions (black dots) and using the warp update of equation (7).
Fig. 3(a) indicates that the dense cost has a relatively wide basin of convergence. The optimization trajectories of BFGS show it takes large steps in directions of maximum gradient, then makes abrupt changes in direction. ECCM and NLS both take smooth paths to the minimum.
The sparse, locally normalized cost surface (Fig. 3(b)) has a narrower basin of convergence, and more local minima. NLS and BFGS behave similarly to with the dense cost, though NLS takes shorter steps. ECCM behaves erratically, but in this 2-d space nevertheless often converges in the end.
The robustified sparse cost surface (Fig. 3(c)) has an even narrower trough around the minimum, but flatter plateaus and fewer local minima either side, compared to the unrobustified sparse cost. The optimizer trajectories are similar in nature to those of the unrobustified sparse cost.
The smoother behaviour of the NLS optimizer leads to an improved convergence rate in the higher dimensional space of homography fitting, where a greater number of local minima trap erratic optimizers.
3.4 Robustness to occlusion
Figure 7 shows the rate of convergence for each of our 3 NCC costs, using the ESM NLS optimizer, for the case that one quadrant of the source patch has been replaced with random texture, emulating occlusion of the scene by another object. The graph shows that the robustified sparse cost performs significantly better than the other costs, except at starting distances greater than 7 pixels, where the dense cost, with it’s broader convergence basin, still performs better.
Rate of convergence (%)
3.5 Lighting invariance compared to other costs
We compare our new sparse, robust cost formulation against two other locally lighting invariant photometric costs: descriptor fields  (first order version) and the census bit planes transform . Both are examples of image transformation approaches, and showed state of the art results when published.666 A deep learned transformation
A deep learned transformation has appeared since, but the proprietary model has not been released publicly. We implemented and evaluated both within our framework, using Gauss-Newton ESM to optimize both, and compare to our NLS ESM method in Figure 8. The results show (Fig. 8(a)) that convergence rates for our method are best overall. However, these convergence rates are slightly worse than census bit planes at small starting distances, which then drop off rapidly, and slightly worse than descriptor fields at larger starting distances, possibly due to the latter’s image blurring. In terms of optimization time (Fig. 8(b)), census bit planes is significantly slower, using 8 channels compared to NCC’s 1. Descriptor fields uses 4 channels, but is comparable in speed to our method. This implies convergence in fewer iterations, which could again be due to its image blurring.
This work shows that NCC optimization in a non-linear least squares framework is not only possible, but straightforward to implement, efficient, and effective at optiming NCC. It also shows that the sparse, robust formulation introduced here better registers images with local intensity variations than both the standard NCC formulation, and other locally lighting invariant methods tested, as well as being more robust to occlusions.
Our experiments focused on optimizations at a single resolution, whereas in practice such methods are generally used in a multi-resolution framework in order to broaden convergence. However, since each level of a multi-resolution framework is an instance of a single resolution optimization, our approach and results are directly applicable to these frameworks also.
We note that our NCC least squares framework can be extended to the following scenarios: data fields of arbitrary dimension and/or multiple channels, and methods which don’t just solve for registration parameters, but also solve for other variables, such as the data coordinates , as in direct odometry .
Appendix A NCC least squares optimization
a.1 Formulation equivalence
The standard formulation of NCC [33, eqn. 8.11] is the following score between two patches, which should be maximized:
We now show that it has the same shape (negated, up to scale and ignoring a scalar offset) as the least squares NCC cost between two patches [14, eqn. 4], which is minimized:
Therefore not only is the optimum of each score achieved at the same inputs, but the gradients of the score functions are the same (up to a scale factor of -2) for all inputs.
a.2 NCC Jacobian
The analytic Jacobian can be computed first by subtracting the mean, then by applying the length normalization using the quotient rule, as follows:
. Applying the chain rule and some rearrangement (using), this gives
a.3 One-pass inverse compositional optimization
In an inverse compositional framework, the Jacobian is pre-computed. Therefore, only the patch (not derivatives) require NCC normalization at each iteration. In the case where , it is more efficient to apply this normalization to than to . This section shows how this can be achieved for the formulation of equation (1).
First note that a matrix whose columns sum to zero, , can be post-multiplied by any matrix and still have zero-sum columns: ; the same goes for pre-multiplication of a matrix with zero-sum rows. jac_zero_mean has zero-sum columns, therefore so does , and its pseudo-inverse, has zero-sum rows. Therefore, for Gauss-Newton optimization,
Both and can be computed during the first pass over patch samples, using the identity for the latter. The division and subtraction of one_pass is then applied to N coefficients only. The result is that patch samples need not be stored and revisited in a second pass to apply normalization. This one-pass approach improves computational efficiency, especially in the case of large patches, which might otherwise incur cache misses on a second pass.
Note that this approach cannot be applied to sparse_framework, since each cost must be normalized separately.
B. D. Lucas and T. Kanade, “An iterative image registration technique with an
application to stereo vision,” in
Proceedings of the International Joint Conference on Artificial Intelligence, 1981, pp. 674–679.
G. D. Hager and P. N. Belhumeur, “Efficient region tracking with parametric models of geometry and illumination,”IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 20, no. 10, pp. 1025–1039, 1998.
S. Baker and I. Matthews, “Lucas-Kanade 20 years on: A unifying framework,”
International Journal of Computer Vision, vol. 56, no. 3, pp. 221–255, 2004.
-  S. Benhimane and E. Malis, “Real-time image-based tracking of planes using efficient second-order minimization,” in Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, vol. 1, 2004, pp. 943–948.
-  J. Engel, V. Koltun, and D. Cremers, “Direct sparse odometry,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 4, 2017.
-  I. Matthews and S. Baker, “Active appearance models revisited,” International Journal of Computer Vision, vol. 60, no. 2, pp. 135–164, 2004.
-  D. I. Barnea and H. F. Silverman, “A class of algorithms for fast digital image registration,” IEEE Transactions on Computers, vol. C-21, no. 2, pp. 179–186, Feb 1972.
-  Y. Furukawa and J. Ponce, “Accurate, dense, and robust multiview stereopsis,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 32, no. 8, pp. 1362–1376, 2010.
-  J. L. Schönberger, E. Zheng, J.-M. Frahm, and M. Pollefeys, “Pixelwise view selection for unstructured multi-view stereo,” in Proceedings of the European Conference on Computer Vision. Springer, 2016, pp. 501–518.
-  D. Wagner, D. Schmalstieg, and H. Bischof, “Multiple target detection and tracking with guaranteed framerates on mobile phones,” in Proceedings of the IEEE International Symposium on Mixed and Augmented Reality, 2009, pp. 57–64.
-  J. Nocedal and S. J. Wright, Numerical Optimization (second edition). Springer, 2006.
-  M. Irani and P. Anandan, “Robust multi-sensor image alignment,” in Proceedings of the IEEE International Conference on Computer Vision. IEEE, 1998, pp. 959–966.
R. Brooks and T. Arbel, “Generalizing inverse compositional image alignment,”
Proceedings of the International Conference on Pattern Recognition, vol. 2. IEEE, 2006, pp. 1200–1203.
-  G. D. Evangelidis and E. Z. Psarakis, “Parametric image alignment using enhanced correlation coefficient maximization,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 30, no. 10, pp. 1858–1865, 2008.
-  M. Brown and D. G. Lowe, “Automatic panoramic image stitching using invariant features,” International Journal of Computer Vision, vol. 74, no. 1, pp. 59–73, 2007.
-  J. F. Henriques, R. Caseiro, P. Martins, and J. Batista, “High-speed tracking with kernelized correlation filters,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 37, no. 3, pp. 583–596, 2015.
-  Kristan, M. , “The visual object tracking VOT2017 challenge results,” in Proceedings of the IEEE International Conference on Computer Vision Workshops, Oct 2017, pp. 1949–1972.
-  M. Kuse and S. Shen, “Robust camera motion estimation using direct edge alignment and sub-gradient method,” in Proceedings of the IEEE International Conference on Robotics and Automation, 2016, pp. 573–579.
-  X. Wang, W. Dong, M. Zhou, R. Li, and H. Zha, “Edge enhanced direct visual odometry.” in Proceedings of the British Machine Vision Conference, 2016.
-  A. Crivellaro and V. Lepetit, “Robust 3d tracking with descriptor fields,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2014, pp. 3414–3421.
-  E. Antonakos, J. Alabort-i Medina, G. Tzimiropoulos, and S. P. Zafeiriou, “Feature-based lucas–kanade and active appearance models,” IEEE Transactions on Image Processing, vol. 24, no. 9, pp. 2617–2632, 2015.
-  H. Alismail, B. Browning, and S. Lucey, “Robust tracking in low light and sudden illumination changes,” in Proceedings of the International Conference on 3D Vision. IEEE, 2016, pp. 389–398.
-  C.-H. Chang, C.-N. Chou, and E. Y. Chang, “CLKN: Cascaded Lucas-Kanade networks for image alignment,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017.
-  G. Silveira and E. Malis, “Real-time visual tracking under arbitrary illumination changes,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. IEEE, 2007, pp. 1–6.
-  P. Viola and W. M. Wells III, “Alignment by maximization of mutual information,” International Journal of Computer Vision, vol. 24, no. 2, pp. 137–154, 1997.
-  F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and P. Suetens, “Multimodality image registration by maximization of mutual information,” IEEE Transactions on Medical Imaging, vol. 16, no. 2, pp. 187–198, 1997.
-  B. Triggs, P. F. McLauchlan, R. I. Hartley, and A. W. Fitzgibbon, “Bundle adjustment—a modern synthesis,” in International workshop on vision algorithms. Springer, 1999, pp. 298–372.
-  M. J. Black and A. Rangarajan, “On the unification of line processes, outlier rejection, and robust statistics with applications in early vision,” International Journal of Computer Vision, vol. 19, no. 1, pp. 57–91, 1996.
-  E. Eade, “Lie groups for computer vision,” http://ethaneade.com/lie_groups.pdf.
-  P. W. Holland and R. E. Welsch, “Robust regression using iteratively reweighted least-squares,” Communications in Statistics-theory and Methods, vol. 6, no. 9, pp. 813–827, 1977.
-  O. Woodford, “Computer vision datasets,” https://sites.google.com/site/oliverwoodford/datasets.
-  R. Hartley and A. Zisserman, Multiple view geometry in computer vision, 2nd ed. Cambridge University Press, 2004.
-  R. Szeliski, Computer vision: algorithms and applications. Springer Science & Business Media, 2010.