US6597816B1 - Correcting distortion in an imaging system using parametric motion estimation - Google Patents

Abstract

Distortion correction in an imaging system using parametric motion estimation. Parametric motion estimation is used to correct for geometric lens distortion and illumination changes. Models for geometric lens distortion and different illumination models are presented, as is a method for an iterative linearising method of solving the resulting sets of nonlinear equations.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

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.

2. Art Background

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.

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.

SUMMARY OF THE INVENTION

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.

BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

FIG. 3 is a flowchart of the linearized solution.

DETAILED DESCRIPTION

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.

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.

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.

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.

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.

Optical Flow

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

c(t;x,y,t ₁)=(x(t;x,y,t ₁),y(t;x,y,t ₁),t),

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.

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{I\left(c\left(t\right)\right)}{t}=0.& \left(1\right)\end{array}$

Using the chain rule of vector calculus, Equation (1) can be written as $\begin{array}{cc}\frac{I\left(c\left(t\right)\right)}{t}=\frac{\partial I\left(x,y,t\right)}{\partial x}\frac{x\left(t\right)}{t}+\frac{\partial I\left(x,y,t\right)}{\partial y}\frac{y\left(t\right)}{t}+\frac{\partial I\left(x,y,t\right)}{\partial t}\frac{t}{t}=0.& \left(2\right)\end{array}$

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}$

results in

I _x(x,y,t)x′(t)+I _y(x,y,t)y′(t)+I _t(x,y,t)=0,  (4)

where I_x and I_y are the spatial image gradients, I_t is 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:

I _x x′+I _y y′+I _t=0.  (5)

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:

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)

For convenience, define:

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

and

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

to obtain

I _x(x,y)u(x,y)+I _y(x,y)v(x,y)+I ₂₁(x,y)=0.  (9)

As a result of assuming I_x and I_y are 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_x and I_y can 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_x and I_y using either frame.

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_x and I_y using either frame, and the interpretation of this amounts to whether we employ a forward or backward approximation to estimate the derivatives.

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):

I ₁(x ₁ ,y ₂)=I ₂(x ₂ ,I _{y 2}),  (10)

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

x ₂ =x ₁ +u(x ₁ ,y ₁)

y ₂ =y ₁ +v(x ₁ ,y ₁)  (11)

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,

I ₁(x ₁ , y ₁)=I ₂(x ₁ +u(x ₁ , y ₁), y _1+v(x ₁ , y ₁))

$\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}$

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}$

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

I _x(x,y)u(x,y)+I _y(x,y)v(x,y)+I ₂₁(x,y)=0,  (14)

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

Parametric Motion Models

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.

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}$

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.

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.

Translational Motion Estimation

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.

The motion model constraints for this case are

u(x,y)=u _g v(x,y)=v _g  (16)

where u_g and v_g refer 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}$

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}$

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

Affine Motion

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

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

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

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.

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$

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]$

Bilinear Motion

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.

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

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

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

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]$

Perspective Motion

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.

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}-xv\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}$

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

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

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

xy−a₈ y ²  (22)

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$

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}$

where

η=x ² I _x +xyI _y

μ=xyI _x +y ² I _y

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.

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}$

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

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}$

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

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

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}$

The solution of this simple 2 by 2 matrix equation yields the optimal translational motion model parameters u_g and v_g under the given quadratic constraint cost function. Of particular importance is the cost function

φ(u,v)=(u−u _c)²+(v−v _c)²  (28)

where we know the expected translational parameters.

Intensity Changes

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

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

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.

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

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

From which Taylor Series expansion results in:

I _x(x,y)(xa ₁ +ya ₂ +a ₃)+I _y(x,y)(xa ₄ +ya ₅ +a ₆)+I ₂(x,y)−κ₁ I ₁(x,y)−κ₂≅0  (31)

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}$

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}-1I_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]$

Lens Distortion

Inexpensive lens systems introduce geometric distortion into the imaging system. This is shown in FIG. 2. Geometric distortion causes the point shown in undistorted first frame 210 to appear in a slightly different location in geometrically distorted image 220. This same distortion causes the apparent motion of this point in undistorted second frame 230 to appear instead at the position shown in distorted second frame 240.

In this section, we will develop methods of motion estimation between images captured with an imperfect lens that causes geometric distortions. The gamma factor and external illumination changes are also taken into account. A nonlinear optimization framework to solve for distortion parameters is introduced. Linearization of this model results in OFE-like equations.

The geometric lens distortions can be modeled as $\begin{array}{cc}{x}_g=\frac{x}{1+{\kappa }_7\left(x^2+y^2\right)}& \left(33\right)\\ y_g=\frac{y}{1+{\kappa }_7\left(x^2+y^2\right)}& \left(34\right)\end{array}$

where κ₇ is the lens geometric distortion parameter, (x_g,y_g) are the distorted pixel positions and (x,y) are the corrected pixel positions.

Using Equation (33) and (34), we write $\begin{array}{cc}\frac{x_g}{y_g}=\frac{x}{y},& \left(35\right)\end{array}$

Putting Equation (35) into Equation (33) yields $\begin{array}{cc}{y}^2\left[{\kappa }_7\left(\frac{x_g^2}{y_g^2}+1\right)\right]+y\left(-\frac{1}{y_g}\right)+1=0& \left(36\right)\end{array}$

There are two solutions to this quadratic equation: $\begin{array}{cc}y=\frac{1\pm \sqrt{1-4{\kappa }_7\left(x_g^2+y_g^2\right)}}{2{\kappa }_7\left(x_g^2+y_g^2\right)}y_g,x=y\frac{x_g}{y_g}& \left(37\right)\end{array}$

As κ₇ goes to 0, y should approach y_g. Since, $\underset{{\kappa }_7\to 0}{\lim }\{\frac{1+\sqrt{1-4{\kappa }_7\left(x_g^2+y_g^2\right)}}{2{\kappa }_7\left(x_g^2+y_g^2\right)}\}y_g=\infty ,$

this can not be a solution. On the other hand, the second solution satisfies this requirement, hence it is the true solution. $\underset{{\kappa }_7\to 0}{\lim }\frac{1-\sqrt{1-4{\kappa }_7\left(x_g^2+y_g^2\right)}}{2{\kappa }_7\left(x_g^2+y_g^2\right)}y_g=y_g$

There are many parametric motion models that map pixels from one frame to another. The affine model is used for subsequent derivation in this section. The translational, bilinear, and perspective motion models follow exactly the same procedure, and are therefore not treated for the sake of space.

Using the geometric distortion model given by Equation (33), we write $\begin{array}{cc}{x}_g^{\prime }=\left({\kappa }_1,{\kappa }_2,{\kappa }_3,{\kappa }_4,{\kappa }_5,{\kappa }_6,{\kappa }_7,x_g,y_g\right)=\frac{x^{\prime }}{1+{\kappa }_7\left(x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}\right)}& \left(38\right)\end{array}$

$\begin{array}{cc}{y}_g^{\prime }=\left({\kappa }_1,{\kappa }_2,{\kappa }_3,{\kappa }_4,{\kappa }_5,{\kappa }_6,{\kappa }_7,x_g,y_g\right)=\frac{y^{\prime }}{1+{\kappa }_7\left(x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}\right)}& \left(39\right)\end{array}$

where

y′=κ ₄ x(x _g ,y _g)+κ₂ y(x _g ,y _g)+κ₂ y(x _g ,y _g)+κ₃ y=κ ₄ x(x _g ,y _g)+κ₅ y(x _g ,y _g)+κ₆  (40)

where x(x_g,y_g) and y(x_g,y_g) are as given by Equation (37).

Therefore, we have a nonlinear mapping from one frame to another characterized by 7 parameters. With the assumption that intensity does not change along the motion trajectory, we write:

I ₁(x _g ,y _g)=I ₂(x′ _g ,y′ _g)  (41).

Nonlinear Geometric Correction

We can define a cost function given by $\begin{array}{cc}\Psi =\sum _{x_g,y_g\in I_1}{\left(I_1\left(x_g,y_g\right)-I_2\left(x_g^{\prime },y_g^{\prime }\right)\right)}^2.& \left(42\right)\end{array}$

This is a function of the affine and lens-distortion parameters,

k=(κ₁,κ₂,κ₃,κ₄, κ₅,κ₆,κ₇)

A variant of Powell's method for multivariable nonlinear optimization is used to find the parameters that minimize the cost function. Other multivariable nonlinear optimization methods known in the art may also be used. To speed up the calculation, rather than summing over all x_g,y_g ε I₁, summation over several small rectangular regions in the image may be performed. The algorithm's performance depends on the size, location, and number of regions used. The accuracy of the estimate may be improved by choosing regions with features such as corners, rather than uniform regions, and by choosing regions towards the periphery of the image (where warping is more pronounced). Nonlinear optimization methods also benefit from a good initial guess. The affine-motion estimation method is used on the middle portions of the images (where the warping is relatively mild) for initial estimates of the affine parameters.

The widely used lens distortion model given by x_g=x(1+κ₇(x²+y²)) is not only a cruder approximation of the physics of the lens, but also give rise to more complicated equations, which requires the calculation of roots of a cubic polynomial.

A hybrid system that alternates between the affine-motion estimation method to estimate the six affine motion parameters, and a one-dimensional nonlinear optimization technique to estimate the lens-distortion parameter is used for optimization. For the one-dimensional optimization, the “Golden Search” method is used to minimize the cost function Ψ, here considered as a function of the single parameter κ₇ with the affine parameters (κ₁,κ₂,κ₃,κ₄,κ₅,κ₆) clamped.

Linearization in Geometric Correction

We may also use a Taylor Series approximation to get a set of linear equations with a closed-form solution. This process is shown in FIG. 3. Once the nonlinear model is selected 310, giving a system of nonlinear equations, an initial set of parameters is selected 320. Given the set of parameters, the system of nonlinear equations is linearized around the set of parameters 330, resulting in a linear system. This linear system of equations is solved 340 using standard techniques to reach a new set of parameters. The process of linearizing the system of nonlinear equations around the set of parameters and resolving the resulting linear solution is iterated 350 until our convergence criteria is met, or for a set number of steps. The linearization 330 is done as shown previously by Taylor Series expansion. If we were to center the Taylor Series around the parameter set k_def=(1,0,0,0,1,0,0) as we have done in previous sections, we would get no information about the lens-distortion parameter κ₇, since. ${\frac{\partial }{\partial {\kappa }_7}}_{k_{\mathrm{def}}}=0\mathrm{and}{\frac{\partial }{\partial {\kappa }_7}}_{k_{\mathrm{def}}}=0.$

This is expected, since k_def=(1,0,0,0,1,0,0) induces zero motion, and with no motion, both images are the same, so we have no hope of recovering the lens-distortion parameter. Likewise, a Taylor Series approximation around a parameter-set corresponding to rotation around the center of the image provides no information on the lens distortion, as the distortion is radial. But other affine-parameter sets (κ₁,κ₂,κ₃,κ₄,κ₅,κ₆) yield an observable system for the lens-distortion parameter. More precisely, a Taylor Series approximation of I₂(x′_g,y′_g) around the parameters K_w=(κ₁,κ₂,κ₃,κ₄,κ₅,κ₆, 0) would yield $\begin{array}{cc}{I}_2\left(x_g^{\prime },y_g^{\prime }\right)=I_2\left(x_g,y_g\right)+\sum _{i=1}^7{\frac{\partial I_2}{\partial x}\frac{\partial }{\partial {\kappa }_i}}_{k_w}+\sum _{i=1}^7{\frac{\partial I_2}{\partial y}\frac{\partial }{\partial {\kappa }_i}}_{k_w}& \left(43\right)\end{array}$

where $\begin{array}{ccc}{\frac{\partial }{\partial {\kappa }_1}}_{k_w}=x_g& {\frac{\partial }{\partial {\kappa }_2}}_{k_w}=y_g& {\frac{\partial }{\partial {\kappa }_2}}_{k_w}=1\\ {\frac{\partial }{\partial {\kappa }_4}}_{k_w}=0& {\frac{\partial }{\partial {\kappa }_5}}_{k_w}=0& {\frac{\partial }{\partial {\kappa }_6}}_{k_w}=0\end{array}$ $\begin{array}{c}{\frac{\partial }{\partial {\kappa }_7}}_{k_w}=\left(-k_1^3+k_1-k_1k_4^2\right)x_g^3+\\ \left(-2k_1k_4k_5-3k_1^2k_2-k_2k_4^2+k_2\right)x_g^2y_g+\\ \left(-3k_1k_2^2+k_1-2k_2k_4k_5-k_1k_5^2\right)x_gy_g^2+\\ \left(-k_2^3+k_2-k_2k_5^2\right)y_g^3+\left(-3k_1^2k_3-2k_1k_4k_6-k_3k_4^2\right)x_g^2+\\ \left(-6k_1k_2k_3-2k_2k_4k_6-2k_3k_4k_5-2k_1k_5k_6\right)x_gy_g+\\ \left(-2k_2k_5k_6-k_3k_5^2-3k_2^2k_3\right)y_g^2+\\ \left(-3k_1k_3^2-k_1k_6^2-2k_3k_4k_6\right)x_g+\\ \left(-k_2k_6^2-3k_2k_3^2-2k_3k_5k_6\right)y_g+\left(-k_3^3-k_3k_6^2\right)\end{array}$ $\begin{array}{ccc}{\frac{\partial }{\partial {\kappa }_1}}_{k_w}=0& {\frac{\partial }{\partial {\kappa }_2}}_{k_w}=0& {\frac{\partial }{\partial {\kappa }_2}}_{k_w}=0\\ {\frac{\partial }{\partial {\kappa }_4}}_{k_w}=x_g& {\frac{\partial }{\partial {\kappa }_5}}_{k_w}=y_g& {\frac{\partial }{\partial {\kappa }_6}}_{k_w}=1\end{array}$ $\begin{array}{c}{\frac{\partial }{\partial {\kappa }_7}}_{k_w}=\left(-k_4k_1^2-k_4^3+k_4\right)x_g^3+\\ \left(-k_5k_1^2-3k_4^2k_5-2k_4k_1k_2+k_5\right)x_g^2y_g+\\ \left(-2k_5k_1k_2+k_4-k_4k_2^2-3k_4k_5^2\right)x_gy_g^2+\\ \left(-k_5k_2^2-k_5^3+k_5\right)y_g^3+\left(-k_6k_1^2-3k_4^2k_6-2k_4k_1k_3\right)x_g^2+\\ \left(-6k_4k_6k_5-2k_4k_3k_2-2k_5k_1k_3-2k_6k_1k_2\right)x_gy_g+\\ \left(-2k_5k_3k_2-k_6k_2^2-3k_6k_5^2\right)y_g^2+\\ \left(-2k_6k_1k_3-3k_4k_6^2-k_4k_3^2\right)x_g+\\ \left(-3k_6^2k_5-2k_6k_3k_2-k_5k_3^2\right)y_g+\left(-k_6^3-k_6k_3^2\right)\end{array}$

The method described earlier is used for an initial estimate of the affine parameters, which are then used as the center of the Taylor series expansion. If the initial estimate is almost entirely translational, one may expand the Taylor series around just the translational component of the initial estimate to simplify the computations, taking k_w=(1,0,κ₃, 0,1,κ₆,0). This gives ${\frac{\partial }{\partial {\kappa }_7}}_{k_w}={\mathrm{–3k}}_3x_g^2-2k_6x_gy_g-k_3y_g^2-\left(k_6^2+3k_3^2\right)x_g-2k_3k_6y_g-\left(-k_3^3+k_3k_6^2\right)\mathrm{and}{\frac{\partial }{\partial {\kappa }_7}}_{k_w}=-k_6x_g^2-2k_3x_gy_g-3k_6y_g^2-2k_6k_3x_g-\left(3k_6^2+k_3^2\right)y_g-\left(k_6^3+k_6k_3^2\right)$

but may give inaccurate results if there is a substantial non-translational aspect to the motion.

Substituting Equation (43) into Equation (42), and taking the derivatives and setting them equal to zero yields the following linear matrix equation: $[\begin{array}{ccccccc}\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}\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}\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}\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}\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}\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}\mu I_y\\ \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{\mu \mu }\end{array}][\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\Delta k_1\\ \Delta k_2\end{array}\\ \Delta k_3\end{array}\\ \Delta k_4\end{array}\\ \Delta k_5\end{array}\\ \Delta k_6\end{array}\\ \Delta k_7\end{array}]=\left[\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_w\\ -\sum _{x,y\in I}{\mathrm{yI}}_xI_w\end{array}\\ -\sum _{x,y\in I}I_xI_w\end{array}\\ -\sum _{x,y\in I}{\mathrm{xI}}_yI_w\end{array}\\ -\sum _{x,y\in I}{\mathrm{yI}}_yI_w\end{array}\\ -\sum _{x,y\in I}I_yI_w\end{array}\\ -\sum _{x,y\in I}\mu I_w\end{array}\right]$

where ${\mu ={\frac{\partial }{\partial {\kappa }_7}}_{k_w}I_x+\frac{\partial }{\partial {\kappa }_7}}_{k_w}$

I_y and I_w=I₂(x,y)−I₁(x_w,y_w), with x_w and y_w the images of x and y under the transformation given by k_w=(κ₁,κ₂,κ₃,κ₄,κ₅,κ₆,0).

Illumination Correction

With the incorporation of external illumination change and offset variation in capture devices as explained earlier, we modify Equation (41) as

κ₈ I ₁(x _g ,y _g)+κ₉ =I ₂(x′ _g y′ _g),  (44)

and we redefine the cost function in terms of motion, lens distortion, and illumination change parameters as: $\begin{array}{cc}\Psi =\sum _{x_g,y_g\in I_1}{\left({\kappa }_8I_1\left(x_g,y_g\right)+{\kappa }_9-I_2\left(x_g^{\prime },y_g^{\prime }\right)\right)}^2& \left(45\right)\end{array}$

To account for vignetting effect, the vignetting correction function is introduced, and Equation (44) modified as $\begin{array}{cc}\frac{{\kappa }_8I_1\left(x_g,y_g\right)+{\kappa }_9}{W\left(x_g,y_g\right)}=\frac{I_2\left(x_g^{\prime },y_g^{\prime }\right)}{W\left(x_g^{\prime },y_g^{\prime }\right)}& \left(46\right)\end{array}$

where $W\left(x_g,y_g\right)={\kappa }_{10}{\cos }^4\left(\arctan \left(\frac{\sqrt{x_w^2+y_w^2}}{f}\right)\right)+\left(1-{\kappa }_{10}\right)$

where f is the camera focal length, and (x_w,y_w) is the camera imaging plane coordinates of the pixel (x_g,y_g).

Then we can define a cost function in terms of (κ₁,κ₂, κ₃,κ₄, κ₅,κ₆,κ₇,κ₈,κ₉,κ₁₀) as $\begin{array}{cc}\Psi =\sum _{x_g,y_g\in I_1}{\left(\frac{{\kappa }_{10}I_1\left(x_g,y_g\right)+{\kappa }_{11}}{W\left(x_g,y_g\right)}-\frac{I_2\left(x_g^{\prime },y_g^{\prime }\right)}{W\left(x_g^{\prime },y_g^{\prime }\right)}\right)}^2& \left(47\right)\end{array}$

This cost function can again be minimized via non-linear optimization techniques, or we can perform Taylor Series approximation to arrive at a set of linear equations, which can easily be solved.

The foregoing detailed description of the present invention is provided for the purpose of illustration and is not intended to be exhaustive or to limit the invention to the precise embodiments disclosed. Accordingly the scope of the present invention is defined by the appended claims.

Claims (16)

What is claimed is:

1. A method for correcting distortion in an imaging system, the method comprising:

a) applying a motion estimation model to a first and a second image to generate a set of nonlinear equations representing distortion, wherein the motion estimation model accommodates a change in illumination along at least one motion trajectory between the first and second images;

b) solving the set of nonlinear equations to produce a set of correction parameters; and

c) correcting at least one of the first and second images using the correction parameters to remove distortion.

2. The method of claim 1 where the distortion includes geometric lens distortion.

3. The method of claim 2 wherein geometric lens distortion is modeled by the transform: $x_g=\frac{x}{1+{\kappa }_7\left(x^2+y^2\right)}y_g=\frac{y}{1+{\kappa }_7\left(x^2+y^2\right)}$

where (x,y) are the corrected pixel locations, (x_g,y_g) are distorted pixel locations, and κ₇ is a geometric lens distortion factor.

4. The method of claim 1 wherein the correction parameters correct a change in illumination along the at least one motion trajectory.

5. The method of claim 4 wherein illumination changes are modeled with a linear model.

6. The method of claim 4 wherein illumination changes are modeled with a Vignetting model.

7. The method of claim 1 where the method of solving the set of nonlinear equations is a multivariable nonlinear optimization method.

8. The method of claim 1 wherein step b) further comprises:

a) selecting an initial set of solution parameters;

b) linearizing the set of nonlinear equations around the set of solution parameters;

c) solving the resulting linearized set of equations to reach an updated set of solution parameters; and

d) repeating steps b) and c) until a condition defined by a selected one of a finite number of iterations and a predetermined convergence criteria is met.

9. The method of claim 8 where the set of nonlinear equations includes correction for geometric lens distortion.

10. The method of claim 8 where the set of nonlinear equations includes correction for illumination changes.

11. The method of claim 8 where the set of nonlinear equations includes correction for both geometric lens distortion and illumination changes.

12. The method of claim 11 where the model for illumination changes is a linear model.

13. The method of claim 11 where the model for illumination changes is a Vignetting model.

14. The method of claim 11 wherein geometric lens distortion is modeled by the transform: $x_g=\frac{x}{1+{\kappa }_7\left(x^2+y^2\right)}y_g=\frac{y}{1+{\kappa }_7\left(x^2+y^2\right)}$

where (x,y) are the corrected pixel locations, (x_g,y_g) are distorted pixel locations, and κ₇ is the distortion factor.

15. A computer storage medium having instructions stored therein for programming a computer to perform the method of claim 11.

16. The method of claim 1 wherein step a) further comprises the step of performing motion estimation modeling using at least one of a translational motion model, an affine motion model, a bilinear motion model, and a perspective motion model.

Images