Correcting distortion in an imaging system using parametric motion estimation

[0002] 1. Field of the Invention

[0003] The present invention deals with the correction of distortion in imaging systems, more particularly to systems such as scanners and in cameras used in applications such as image capture and video conferencing.

[0004] 2. Art Background

[0005] With the advent of inexpensive imaging sensors, applications such as video conferencing and remote image capture are becoming more prevalent. As manufacturers provide imaging systems at lower and lower price points, compromises are made especially in the optics provided, leading to systems having geometric lens distortion. As these systems are placed in more locations, problems with illumination arise, both as a result of variations in lighting conditions, as well as the effect of automatic brightness controls in the imaging systems themselves.

[0006] The field of motion estimation deals with processing images for example to stitch together successive overlapping images to form a unified mosaic. These motion estimation techniques deal with arrays of pixels, and operate by deriving the motion between a group of pixels in one frame and the same pixels in a new position in a second frame. For motion estimation to succeed, the illumination levels and therefore pixel intensity between the two frames must be the same, and the frames must be relatively free of geometric distortion. What is needed is a way to correct for the effects of geometric lens distortion in motion estimation, and to correct for changes in illumination that effect pixel density.

[0007] Motion estimation techniques and the Optical Flow Equation (OFE) are applied to correct for distortion in an imaging system. Geometric lens distortion is corrected using a novel optical correction model. Illumination changes are corrected using both simple and Vignetting models. The solutions involve solving of a set of nonlinear equations. Typical solutions to these equations are computationally expensive. An iterative linearising solution is presented.

[0008] The present invention is described with respect to particular exemplary embodiments thereof and reference is made to the drawings in which:

[0009] FIG. 1 is a block diagram of an imaging system,

[0010] FIG. 2 is a diagram showing geometric lens distortion, and

[0011] FIG. 3 is a flowchart of the linearized solution.

[0012] Inexpensive imaging sensors are enabling more applications of imaging technologies, such as video conferencing and image scanning. Compromises in sensor design, and in sensor placement lead to distortions in the resulting images.

[0013] Motion estimation techniques deal with processing images for example to stitch together successive overlapping images to form a unified mosaic; this technique us used for example in assembling scanned images of a large area using multiple passes of an imaging sensor. Motion estimation allows the overlaps to be recognized and the smaller images from each pass stitched together to form the larger unified mosaic. Similar techniques also arise in videoconferencing when the videoconferencing camera is moved, or a videoconferencing camera is used to scan a document.

[0014] To succeed, motion estimation presupposes the illumination and therefore pixel densities of overlapping image portions to be the same, and that the images are free from geometric distortion. Inexpensive optics introduce geometric distortion, and may also introduce illumination differences such as vignetting, the gradual fall off of illumination or image density away from the central axis of the image sensor. Additional illumination distortions are introduced by automatic brightness controls in image sensors, and by poor or uneven illumination in the environment.

[0015] Two dimensional motion estimation refers to determining the displacements or velocities of pixels from one frame to another. The motion is induced as the same physical point on the surface of an object in 3-D space is mapped to different locations within the imaging plane at different time. In the sampled image domain this results in different pixel indices. This difference in pixel indices taken relative to one of the frames is called the “motion vector” at that pixel between those two frames. Each pixel in the image has a motion vector of dimension two representing horizontal and vertical displacements. The collection of all motion vectors is sometimes called the “motion vector field”, or simply “motion field”. The term “field” should be understood in a broader sense, not in strict mathematical definition of “vector field”. Motion estimation from sequences of images is discussed for example in A. M. Tekalp, Digital Video Processing, Prentice Hall Signal Processing Series; M. I. Sezan and R. L. Lagendijk, Motion Analysis and Image Sequence Processing, Norwell, Mass.: Kluwer, 1993; and J. K. Aggarwal and N. Nandhakumar, “On the computation of motion from sequences of images- A review,” Proc. IEEE, vol. 76, pp. 917-935, August 1988.

[0016] FIG. 1 shows a system for performing motion estimation. Bus 100 interconnects image capture device 110 with central processing unit (CPU) 120. Memory hierarchy 130 provides storage for instructions and data, including storage of frames if data from information capture device 110, typically containing a mixture of RAM, ROM, and permanent disc storage. Not shown are other typical items such as displays and network connections. Image capture device 110 may be a device such as a television camera with appropriate frame capture hardware, or other imaging device such as a teleconferencing camera, CCD image sensor, scanner, or the like.

[0017] In principle, all pixels within an image can be thought to follow a given trajectory through time. The trajectory of a given pixel can be parameterized via time, given an initial position at some reference time. If we take the reference time as t=t₁, then

[0018] c(t;x,y,t₁)=(x(t;x,y,t₁),y(t;x,y,t₁),t),

[0019] where c denotes the trajectory followed by a pixel originating at (x,y,t₁). For the development that follows, we abbreviate x(t;x,y,t₁) as x(t), and apply similar abbreviations for y(t) and c(t). Even though physical conditions generally limit the variation between trajectories of neighboring pixels, conditions such as occlusion and uncovered background can lead to arbitrary (independent) trajectories amongst neighboring pixels. In this treatment we neglect occlusion and uncovered background.

[0020] Assuming further the intensity/color value of a pixel does not change with time along the motion trajectory, then, with I(x,y,t) denoting the intensity function, $\begin{array}{cc}\frac{dI\left(c\left(t\right)\right)}{dt}=0.& \left(1\right)\end{array}$

[0021] Using the chain rule of vector calculus, Equation (1) can be written as $\begin{array}{cc}\frac{dI\left(c\left(t\right)\right)}{dt}=\frac{\partial I\left(x,y,t\right)}{\partial x}\frac{dx\left(t\right)}{dt}+\frac{\partial I\left(x,y,t\right)}{\partial y}\frac{dy\left(t\right)}{dt}+\frac{\partial I\left(x,y,t\right)}{\partial t}\frac{dt}{dt}=0.& \left(2\right)\end{array}$

[0022] Letting $\begin{array}{cc}{I}_x\left(x,y,t\right)=\frac{\partial I\left(x,y,t\right)}{\partial x}{I}_y\left(x,y,t\right)=\frac{\partial I\left(x,y,t\right)}{\partial y}{I}_t\left(x,y,t\right)=\frac{\partial I\left(x,y,t\right)}{\partial t}& \left(3\right)\end{array}$

[0023] results in

[0024] I_x(x,y,t)x′(t)+I_y(x,y,t)y′(t)+I_t(x,y,t)=0,  (4)

[0025] where I_xand I_yare the spatial image gradients, I_tis the temporal image gradient, and where x′ and y′ are the horizontal and vertical velocity components respectively. Equation (4) is called the optic-flow-equation (OFE), sometimes abbreviated as:

[0026] I_xx′+I_yy′+I_t=0.  (5)

[0027] We can proceed using (4) in two equivalent ways. In the first approach integrate both sides of (4) from t₁to t₂with respect to t, and assume both I_x(x,y,t) and I_y(x,y,t) are constant with respect to time to obtain:

[0028] I_x(x,y,t)(x(t₂)−x(t₁))+I_y(x,y,t)(y(t₂)−y(t₁))+(I(x,y,t₂)−I(x,y,t₁))=0.  (6)

[0029] For convenience, define:

[0030] I₂(x,y)=I(x,y,t₂), I₁(x,y)=I(x,y,t₁), I₂₁(x,y)=I₂(x,y)−I₁(x,y)  (7)

[0031] and

[0032] u(x,y)=x(t₂)−x(t₁), v(x,y)=y(t₂)−y(t₁),  (8)

[0033] to obtain

[0034] I_x(x,y)u(x,y)+I_y(x,y)v(x,y)+I₂₁(x,y)=0.  (9)

[0035] As a result of assuming I_xand I_yare constant with respect to time between t₁and t₂, we drop time dependence from the notation. Second, this assumption taken together with the continuity of intensity assumption and (1) implies that I_t(x,y,t) is constant between t₁and t₂. This allows for trivial integration of I_t(x,y,t) on this interval. Lastly, estimates of I_xand I_ycan be obtained anywhere on [t₁,t₂]. This is important to note because if we begin with just two frames at t₁and t₂, we can discretely estimate I_xand I_yusing either frame.

[0036] The second method of proceeding from (5) is to approximate the derivatives involved by first order differences and then multiply through by (t₂−t₁). Then, once again, if we have just two frames at t₁and t₂, we can discretely estimate I_xand I_yusing either frame, and the interpretation of this amounts to whether we employ a forward or backward approximation to estimate the derivatives.

[0037] In the development to this point it has been assumed that there is an explicit temporal evolution of the image intensity distribution. To enforce this assumption the chain rule has been applied to the intensity distribution temporal derivative. We can equivalently, however, begin from a completely discrete perspective and proceed instead by applying a 2-D spatial Taylor series approximation. To that end, assume from the onset that we have two captured frames at times t₁and t₂. Now, the notion of intensity invariance is expressed as (ignoring illumination change effects):

[0038] I₁(x₁,y₂)=I₂(x₂,I_y2),  (10)

[0039] where I₁(x₁,y₁) refers to the image intensity at the pixel (x₁,y₁) in the first frame captured at time t₁, and I₂(x₂,y₂) refers to the intensity at the point (x₂,y₂) in the second frame captured at time t₂. Note that (x₁,y₁) are integer valued, but (x₂,y₂) are real-valued. Both I₁and I₂are assumed to be continuous and differentiable functions with respect to x and y. Using the definition of motion vectors, one can write

[0040] x₂=x₁+u(x₁,y₁)

[0041] y₂=y₁+v(x₁,y₁)  (11)

[0042] where u(.,.) and v(.,.) are the horizontal and vertical motion vectors, respectively. Combining Equations (10) and (11), and using a Taylor Series Expansion of differentiable functions,

[0043] I₁(x₁, y₁)=I₂(x₁+u(x₁, y₁), y_1+v(x₁, y₁))

[0044] $\begin{array}{cc}\begin{array}{c}{I}_1\left(x_1,y_1\right)=I_2\left(x_1+u\left(x_1,y_1\right),y_1+v\left(x_1,y_1\right)\right)\\ =I_2\left(x_1,y_1\right)+\frac{\partial I_2\left(x,y\right)}{\partial x}{{|}_{x_1,y_1}u\left(x_1,y_1\right)+\frac{\partial I_2\left(x,y\right)}{\partial y}|}_{x_1,y_1}v\left(x_1,y_1\right)+h.o.t\end{array}& \left(12\right)\end{array}$

[0045] The Taylor series is expanded around the pixel (x₁,y₁). Assuming either that u(x₁,y₁) and v(x₁,y₁) are small quantities or that the higher order terms are all zero, we can disregard the higher order terms in Equation (12), and arrive at $\begin{array}{cc}\frac{\partial I_2\left(x,y\right)}{\partial x}{{|}_{x_1,y_1}u\left(x_1,y_1\right)+\frac{\partial I_2\left(x,y\right)}{\partial y}|}_{x_1,y_1}v\left(x_1,y_1\right)+\left(I_2\left(x_1,y_1\right)-I_1\left(x_1,y_1\right)\right)\cong 0.& \left(13\right)\end{array}$

[0046] Note that by employing the same notation as in the previous development we write

[0047] I_x(x,y)u(x,y)+I_y(x,y)v(x,y)+I₂₁(x,y)=0,  (14)

[0048] which is the same result (9) previously arrived at.

[0049] In many realistic imaging scenarios a 2-D motion field is induced by either camera motion or objects in the 3-D world undergoing a rigid motion. Each of these effects produces a smoothly varying 2-D motion field (provided that the depth of the objects does not change swiftly in proportion to its distance to imaging device). Noting the fact that the motion field is smoothly varying, there have been attempts to describe it using a parametric model. Namely,the motion vector at any pixel needs to obey a model characterized by a few parameters. Therefore, the motion vector of a pixel is completely represented by the model parameters and its location.

[0050] Given a vector of model parameters a the motion is expressed as u(x,y)=f(x,y;a) and v(x,y)=g(x,y;a). In order to estimate the model parameters we define the following cost function, $\begin{array}{cc}\Psi =\sum _{x,y\in I}\left\{{\left(I_{x}u+I_{y}v+I_{21}\right)}^2\right\}.& \left(15\right)\end{array}$

[0051] The parameter estimation procedure is to take derivatives of ψ with respect to each element of a, set those derivatives to zero, then solve the resulting set of equations. For non-linear models f and g a linearization step is typically employed.

[0052] In the following, we will examine the parametric motion models: 1) translational motion, 2) affine motion, 3) bilinear motion, and 4) perspective motion. These parametric motion models can be shown to be the projection of the motion of simple 3-D objects, such as a plane, under different projection models onto the 2-D imaging plane. Lastly, (15) is extended to include a-priori information and apply the result to estimate constrained global translational motion.

[0053] One of the simplest models that can be imposed on the motion vector field is that all the motion vectors be the same, i.e. every pixel in the image is moving in the same direction by the same amount. This motion model is usually called “global translational motion”. However, it should be noted that translational motion does not necessarily mean that objects in the 3-D (real) world translate. 3-D translational motion could only be reduced to a 2-D translational model under the assumption that the 3-D object is a planar object, and an orthographic camera model is utilized. Additionally, global translational motion can be induced by a camera translation while a “distant” scene is being imaged. It may sound a very restrictive motion model, however there are applications to which this model perfectly fits. An example is multi-frame resolution enhancement technologies for scanner applications, where multiple scans of a flat page are taken such that a translational shift occurs between each scan.

[0054] The motion model constraints for this case are

[0055] u(x,y)=u_gv(x,y)=v_g  (16)

[0056] where u_gand v_grefer to global horizontal and vertical displacements. Then, taking the derivatives of (15)with respect to u and v, and setting them equal to zero yields $\begin{array}{cc}\frac{\partial \Psi }{\partial u}=\sum _{x,y\in I}2\left(I_x{u}_g+I_y{v}_g+I_{21}\right)I_x=0\frac{\partial \Psi }{\partial v}=\sum _{x,y\in I}2\left(I_x{u}_g+I_y{v}_g+I_{21}\right)I_y=0& \left(17\right)\end{array}$

[0057] Equation (17) can be re-written in the matrix form as follows: $\begin{array}{cc}\left[\begin{array}{cc}\sum _{x,y\in I}{I}_x^2& \sum _{x,y\in I}{I}_x{I}_y\\ \sum _{x,y\in I}{I}_x{I}_y& \sum _{x,y\in I}{I}_y^2\end{array}\right]\left[\begin{array}{c}{u}_g\\ v_g\end{array}\right]=\left[\begin{array}{c}-\sum _{x,y\in I}{I}_x{I}_{21}\\ -\sum _{x,y\in I}{I}_y{I}_{21}\end{array}\right]& \left(18\right)\end{array}$

[0058] Equation (18) can be solved to estimate global displacement parameters.

[0059] Affine motion is another commonly employed motion model. It can account for 2-D transformations such as global rotation, translation, zoom and skew. Physically, affine motion is generated between two images when the 3-D scene is planar and two images are acquired using different imaging plane positions under the orthogonal projection. The affine motion model is described by

[0060] u(x,y)=a₁x+a₂y+a₃

[0061] i v(x,y)=a₄x+a₅y+a₆  (19)

[0062] where a₁,a₂,a₃,a₄,a₅,a₆are model parameters (sometimes called “Affine Parameters”) and (x,y) is the pixel location. Note that global translational motion results when a₁, a₂, a₄and a₅are all zero.

[0063] To solve for optimal affine parameters in the Least Square sense, we expand (15) as $\Psi =\sum _{x,y\in I}{\left\{I_x\left(x,y\right)\left({\mathrm{xa}}_1+{\mathrm{ya}}_2+a_3\right)+I_y\left(x,y\right)\left({\mathrm{xa}}_4+{\mathrm{ya}}_5+a_6\right)+I_{21}\left(x,y\right)\right\}}^2$

[0064] As previously mentioned, ψ is minimized by taking the derivatives with respect to the affine parameters a₁,a₂,a₃,a₄,a₅,a₆, and setting each derivative equal to zero. After carrying out the algebra, one can arrive at the following linear matrix equation. $\left[\begin{array}{cccccc}\sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_x{I}_y\\ \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_y\\ \sum _{x,y\in I}{I}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_x{I}_x& \sum _{x,y\in I}{I}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_x{I}_y\\ \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_y\\ \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_y\\ \sum _{x,y\in I}{I}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_y{I}_x& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_y{I}_y\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{a}_1\\ a_2\end{array}\\ a_3\end{array}\\ a_4\end{array}\\ a_5\end{array}\\ a_6\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}-\sum _{x,y\in I}{\mathrm{xI}}_x{I}_{21}\\ -\sum _{x,y\in I}{\mathrm{yI}}_x{I}_{21}\end{array}\\ -\sum _{x,y\in I}{I}_x{I}_{21}\end{array}\\ -\sum _{x,y\in I}{\mathrm{xI}}_y{I}_{21}\end{array}\\ -\sum _{x,y\in I}{\mathrm{xI}}_y{I}_{21}\end{array}\\ -\sum _{x,y\in I}{I}_y{I}_{21}\end{array}\right]$

[0065] The global bilinear transform is used to improve on the capabilities of the affine transform, yet still yield a linear estimation problem. In the case of an affine transformation parallel lines always remain parallel. This is not the case for motion induced by a perspective-projection imaging model. The perspective model, however, which will be treated next, yields a non-linear problem. The bilinear transform can be a good compromise in that it transforms horizontal and vertical straight lines into straight lines, but does not produce parallel lines for any given parallel pair lines. The potential problem in using a bilinear transform is that lines that are not horizontal or vertical will be transformed to have some degree of curvature. If the curvature becomes severe the transform may not be representative of the actual motion and imaging system.

[0066] The derivation of global bilinear motion model estimation is quite similar to affine motion model estimation. In particular, the motion field is governed by:

[0067] u(x,y)=a₁+a₂x+a₃y+a₄xy

[0068] v(x,y)=a₅+a₆x+a₇y+a₈xy  (20)

[0069] By substituting (20) into (15) and setting the derivatives with respect to model parameters equal to zero, we obtain, $\left[\begin{array}{cccccccc}\sum _{x,y\in I}{I}_x{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_x& \sum _{x,y\in I}{\mathrm{xyI}}_x{I}_x& \sum _{x,y\in I}{I}_y{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_x& \sum _{x,y\in I}{\mathrm{xyI}}_y{I}_x\\ \sum _{x,y\in I}{I}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_y{\mathrm{xI}}_x\\ \sum _{x,y\in I}{I}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_y{\mathrm{yI}}_x\\ \sum _{x,y\in I}{I}_x{\mathrm{xyI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xyI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xyI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_x{\mathrm{xyI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{xyI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xyI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xyI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_y{\mathrm{xyI}}_x\\ \sum _{x,y\in I}{I}_x{I}_y& \sum _{x,y\in I}{\mathrm{xI}}_x{I}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_y& \sum _{x,y\in I}{\mathrm{xyI}}_x{I}_y& \sum _{x,y\in I}{I}_y{I}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_y& \sum _{x,y\in I}{\mathrm{xyI}}_y{I}_y\\ \sum _{x,y\in I}{I}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xyI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xyI}}_y{\mathrm{xI}}_y\\ \sum _{x,y\in I}{I}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xyI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xyI}}_y{\mathrm{yI}}_y\\ \sum _{x,y\in I}{I}_x{\mathrm{xyI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xyI}}_x& \sum _{x,y\in I}{I}_y& \sum _{x,y\in I}{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xyI}}_y\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{a}_1\\ a_2\end{array}\\ a_3\end{array}\\ a_4\end{array}\\ a_5\end{array}\\ a_6\end{array}\\ a_7\end{array}\\ a_8\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\sum _{x,y\in I}{I}_x\\ \sum _{x,y\in I}\end{array}\\ \sum _{x,y\in I}{\mathrm{yI}}_x\end{array}\\ \sum _{x,y\in I}{\mathrm{xyI}}_x\end{array}\\ \sum _{x,y\in I}\end{array}\\ \sum _{x,y\in I}\end{array}\\ \sum _{x,y\in I}{\mathrm{yI}}_y\end{array}\\ \sum _{x,y\in I}\end{array}\right]$

[0070] The perspective motion model is the most versatile of the motion models used thus far. Physically,the perspective model takes into account the 2-D mapping induced among perspective projected images of a plane in 3-D under different focal plane positions. A complication is that the perspective model is nonlinear in its model parameters.

[0071] The motion field in global perspective model is governed by the following equations. $\begin{array}{cc}u\left(x,y\right)=\frac{a_{1}x+a_{2}y+a_3}{a_{7}x+a_{8}y+1}-xv\left(x,y\right)=\frac{a_{4}x+a_{5}y+a_6}{a_{7}x+a_{8}y+1}-y& \left(21\right)\end{array}$

[0072] Before using these model parameters in (15) they are first linearized around the default perspective parameter a=(1,0,0,0,1,0,0,0) to yield

[0073] u(x,y)≅(a₁−1)x+a₂y+a₃−a₇x²−a₈xy

[0074] v(x,y)≅a₄x+(a₅−1)y+a₆−a₇

[0075] xy−a₈y²  (22)

[0076] To solve for optimal perspective motion model parameters in the Least Square sense, we apply (22) to (15) and obtain $\Psi =\sum _{x,y\in I}{\left\{\begin{array}{c}{I}_x\left(x,y\right)\left(\left(a_1-1\right)x+a_{2}y+a_3-a_7{x}^2-a_8\mathrm{xy}\right)+\\ I_y\left(x,y\right)\left(a_{4}x+\left(a_5-1\right)y+a_6-a_7\mathrm{xy}-a_8{y}^2\right)+I_{21}\left(x,y\right)\end{array}\right\}}^2$

[0077] To minimize Ψ, we take the derivatives of Ψ with respect to perspective parameters a₁, a₂, a₃, a₄, a₅, a₆, a₇, and a₈and set them equal to zero. After carrying out the algebra, one can arrive at the following linear matrix equation. $\begin{array}{cc}\left[\begin{array}{cccccccc}\sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_x& \sum _{x,y\in I}\eta {\mathrm{xI}}_x& \sum _{x,y\in I}\mu {\mathrm{xI}}_x\\ \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_x& \sum _{x,y\in I}\eta {\mathrm{yI}}_x& \sum _{x,y\in I}\mu {\mathrm{yI}}_x\\ \sum _{x,y\in I}{\mathrm{xI}}_x{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_x& \sum _{x,y\in I}{I}_x{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_x& \sum _{x,y\in I}{I}_y{I}_x& \sum _{x,y\in I}\eta I_x& \sum _{x,y\in I}\mu I_x\\ \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_y& \sum _{x,y\in I}\eta {\mathrm{xI}}_y& \sum _{x,y\in I}\mu {\mathrm{xI}}_y\\ \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_y& \sum _{x,y\in I}\eta {\mathrm{yI}}_y& \sum _{x,y\in I}\mu {\mathrm{yI}}_y\\ \sum _{x,y\in I}{\mathrm{xI}}_x{I}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_y& \sum _{x,y\in I}{I}_x{I}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_y& \sum _{x,y\in I}{I}_y{I}_y& \sum _{x,y\in I}\eta I_y& \sum _{x,y\in I}\mu I_y\\ \sum _{x,y\in I}{\mathrm{xI}}_x\eta & \sum _{x,y\in I}{\mathrm{yI}}_x\eta & \sum _{x,y\in I}{I}_x\eta & \sum _{x,y\in I}{\mathrm{xI}}_y\eta & \sum _{x,y\in I}{\mathrm{yI}}_y\eta & \sum _{x,y\in I}{I}_y\eta & \sum _{x,y\in I}\mathrm{\eta \eta }& \sum _{x,y\in I}\mathrm{\mu \eta }\\ \sum _{x,y\in I}{\mathrm{xI}}_x\mu & \sum _{x,y\in I}{\mathrm{yI}}_x\mu & \sum _{x,y\in I}{I}_x\mu & \sum _{x,y\in I}{\mathrm{xI}}_y\mu & \sum _{x,y\in I}{\mathrm{yI}}_y\mu & \sum _{x,y\in I}{I}_y\mu & \sum _{x,y\in I}\mathrm{\eta \mu }& \sum _{x,y\in I}\mathrm{\mu \mu }\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{a}_1-1\\ a_2\end{array}\\ a_3\end{array}\\ a_4\end{array}\\ a_5-1\end{array}\\ a_6\end{array}\\ a_7\end{array}\\ a_8\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}-\sum _{x,y\in I}{\mathrm{xI}}_x{I}_t\\ -\sum _{x,y\in I}{\mathrm{yI}}_x{I}_{21}\end{array}\\ -\sum _{x,y\in I}{I}_x{I}_{21}\end{array}\\ -\sum _{x,y\in I}{\mathrm{xI}}_y{I}_{21}\end{array}\\ -\sum _{x,y\in I}{\mathrm{yI}}_y{I}_{21}\end{array}\\ -\sum _{x,y\in I}{I}_y{I}_{21}\end{array}\\ -\sum _{x,y\in I}\eta I_{21}\end{array}\\ -\sum _{x,y\in I}\mu I_{21}\end{array}\right]& \left(23\right)\end{array}$

[0078] where

[0079] η=x²I_x+xyI_y

[0080] μ=xyI_x+y²I_y

[0081] Note that we have utilized the Taylor Series Approximation twice, in contrast to all previous derivations. Until this point, only linearization of the intensity function has been required. However, we had to additionally linearize the motion model to derive the Linear Matrix Equation form of Equation (23). Due to this added linearization step, perspective parameter estimation can be less accurate and stable.

[0082] In some applications, one may have a-priori information relating to the motion model parameters. Such information, for instance, could be a probability distribution of the model parameters. We take this additional information into account by modifying the cost function (15) such that $\begin{array}{cc}\Psi =\sum _{x,y\in I}\left\{{\left(I_{x}u+I_{y}v+I_t\right)}^2+\varphi \left(u,v\right)\right\}& \left(24\right)\end{array}$

[0083] Where the a-priori information is embodied using φ(u,v).

[0084] We will only derive the equations to estimate the global translational motion model parameters under such constraints. In this case we proceed as before and obtain $\begin{array}{cc}\frac{\partial \Psi }{\partial u}=\sum _{x,y\in I}\left\{2\left(I_x{u}_g+I_y{v}_g+I_{21}\right)I_x+\frac{\partial \varphi }{\partial u}\right\}=0\frac{\partial \Psi }{\partial v}=\sum _{x,y\in I}\left\{2\left(I_x{u}_g+I_y{v}_g+I_{21}\right)I_y+\frac{\partial \varphi }{\partial v}\right\}=0& \left(25\right)\end{array}$

[0085] At this point, assume that the constraint cost function is a quadratic function in u and v. This assumption is made to arrive at a linear matrix equation for optimal motion parameter estimation. Although any differentiable cost function would be acceptable, the model parameter solution may otherwise require costly non-linear optimization techniques. In particular, the constraint cost function is in the form

[0086] φ(u,v)=ξ₁+ξ₂u+ξ₃v+ξ₄uv+ξ₅u²+ξ₆v².  (26)

[0087] Then, Equation (25) can be re-written as $\begin{array}{cc}\left[\begin{array}{cc}\sum _{x,y\in I}\left(I_x^2+2{\xi }_5\right)& \sum _{x,y\in I}\left(I_x{I}_y+{\xi }_4\right)\\ \sum _{x,y\in I}\left(I_x{I}_y+{\xi }_4\right)& \sum _{x,y\in I}\left(I_y^2+2{\xi }_6\right)\end{array}\right]\left[\begin{array}{c}{u}_g\\ v_g\end{array}\right]=\left[\begin{array}{c}-\sum _{x,y\in I}\left(I_x{I}_{21}+{\xi }_2\right)\\ -\sum _{x,y\in I}\left(I_y{I}_{21}+{\xi }_3\right)\end{array}\right]& \left(27\right)\end{array}$

[0088] The solution of this simple 2 by 2 matrix equation yields the optimal translational motion model parameters u_gand v_gunder the given quadratic constraint cost function. Of particular importance is the cost function

[0089] φ(u,v)=(u−u_c)²+(v−v_c)²  (28)

[0090] where we know the expected translational parameters.

[0091] All methods above are derived under the premise that the intensity of a pixel does not change along the motion trajectory. As explained earlier, this may not always be true. The external illumination may change and the gain/offset of the image capture device may vary between captured frames. To take these effects into account, we modify Equation (10) such that

[0092] κ₁I₁(x₁,y₁)+κ₂=I₂(x₂,y₂)  (29)

[0093] From here, the optimization proceeds exactly as before, except that parameters κ₁and κ₂need to be estimated as well. We shall now derive the affine motion model parameter estimation equations under the assumption (29). Derivations for other motion models are not repeated since the outline of the derivation is very similar in each case.

[0094] Incorporating the affine motion model equations given by the Equation (19), into (29), we write:

[0095] κ₁I₁(x,y)+κ₂=I₂(x+a₁x+a₂y+a₃, y+a₄x+a₅y+a₆)  (30)

[0096] From which Taylor Series expansion results in:

[0097] I_x(x,y)(xa₁+ya₂+a₃)+I_y(x,y)(xa₄+ya₅+a₆)+I₂(x,y)−κ₁I₁(x,y)−κ₂≅0  (31)

[0098] To solve for the unknown parameters, define a cost function Ψ as below: $\begin{array}{cc}\Psi =\sum _{x_1,y_1\in I}{\left\{\begin{array}{c}{I}_x\left(x,y\right)\left({\mathrm{xa}}_1+{\mathrm{ya}}_2+a_3\right)+I_y\left(x,y\right)\left({\mathrm{xa}}_4+{\mathrm{ya}}_5+a_6\right)+\\ I_2\left(x,y\right)-{\kappa }_1{I}_1\left(x,y\right)-{\kappa }_2\end{array}\right\}}^2& \left(32\right)\end{array}$

[0099] Now take the partial derivatives of Ψ with respect to the affine parameters a₁,a₂,a₃,a₄,a₅,a₆and the illumination change model parameters κ₁,κ₂, and set them equal to zero. After carrying out the algebraic steps, we obtain, $\left[\begin{array}{cccccccc}\sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_x{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_x& \sum _{x,y\in I}-I_1{\mathrm{xI}}_x& \sum _{x,y\in I}-{\mathrm{xI}}_x\\ \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_x{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_x& \sum _{x,y\in I}-I_1{\mathrm{yI}}_x& \sum _{x,y\in I}-{\mathrm{yI}}_x\\ \sum _{x,y\in I}{\mathrm{xI}}_x{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_x& \sum _{x,y\in I}{I}_x{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_x& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_x& \sum _{x,y\in I}{I}_y{I}_x& \sum _{x,y\in I}-I_1{I}_x& \sum _{x,y\in I}-I_x\\ \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_x{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{xI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{xI}}_y& \sum _{x,y\in I}-I_1{\mathrm{xI}}_y& \sum _{x,y\in I}-{\mathrm{xI}}_y\\ \sum _{x,y\in I}{\mathrm{xI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_x{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_y{\mathrm{yI}}_y& \sum _{x,y\in I}-I_1{\mathrm{yI}}_y& \sum _{x,y\in I}-{\mathrm{yI}}_y\\ \sum _{x,y\in I}{\mathrm{xI}}_x{I}_y& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_y& \sum _{x,y\in I}{I}_x{I}_y& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_y& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_y& \sum _{x,y\in I}{I}_y{I}_y& \sum _{x,y\in I}-I_1{I}_y& \sum _{x,y\in I}-I_y\\ \sum _{x,y\in I}{\mathrm{xI}}_x{I}_1& \sum _{x,y\in I}{\mathrm{yI}}_x{I}_1& \sum _{x,y\in I}{I}_x{I}_1& \sum _{x,y\in I}{\mathrm{xI}}_y{I}_1& \sum _{x,y\in I}{\mathrm{yI}}_y{I}_1& \sum _{x,y\in I}{I}_y{I}_1& \sum _{x,y\in I}-I_1{I}_1& \sum _{x,y\in I}-1I_1\\ \sum _{x,y\in I}{\mathrm{xI}}_x& \sum _{x,y\in I}{\mathrm{yI}}_x& \sum _{x,y\in I}{I}_x& \sum _{x,y\in I}{\mathrm{xI}}_y& \sum _{x,y\in I}{\mathrm{yI}}_y& \sum _{x,y\in I}{I}_y& \sum _{x,y\in I}-I_1& \sum _{x,y\in I}-1\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{a}_1\\ a_2\end{array}\\ a_3\end{array}\\ a_4\end{array}\\ a_5\end{array}\\ a_6\end{array}\\ {\kappa }_1\end{array}\\ {\kappa }_2\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}-\sum _{x,y\in I}{\mathrm{xI}}_x{I}_2\\ -\sum _{x,y\in I}{\mathrm{yI}}_x{I}_2\end{array}\\ -\sum _{x,y\in I}{I}_x{I}_2\end{array}\\ -\sum _{x,y\in I}{\mathrm{xI}}_y{I}_2\end{array}\\ -\sum _{x,y\in I}{\mathrm{yI}}_y{I}_2\end{array}\\ -\sum _{x,y\in I}{I}_y{I}_2\end{array}\\ -\sum _{x,y\in I}{I}_1{I}_2\end{array}\\ -\sum _{x,y\in I}{I}_2\end{array}\right]$