US 6,791,616 B2Grant

Image lens distortion correcting method

Issue Date:2004-09-14

•16 Claims

•6 Drawing Sheets

Abstract

An image lens distortion correcting method having step (A) for printing an arbitrary image I₁in a computer; step (B) for photographing the printed image I₁with a camera having a lens distortion and obtaining a photographed image I₂in the computer; step (C) for obtaining a parameter θ such that the image I₁is artificially distorted with the parameter θ and an obtained distorted image I₁^udequals with the photographed image ₂; and step (D) for using the obtained parameter θ to correct the image photographed by the camera. Thereby, (1) no correspondence point needs to be given, (2) no special tool or object is used, and (3) all points on the image are used, so that an internal parameter can automatically be obtained with high precision, and image lens distortion can be corrected with high precision.

Metadata

Assignee

Riken

Inventors

Toru Tamaki
Tsuyoshi Yamamura
Noboru Ohnishi

Application Information

Application Number:US 09/903,619

Filing Date:2001-07-13

Priority Date:2000-09-05

Art Unit:7

Classifications

IPC:

H04N 5225H04N 1700H04N 100G06K 940G01B 900G01B 1100

Field of Search:

34833534834034817534819234820799396 91382275358406356124356394

Patent Drawings (6 sheets)

Description

Background of the Invention

[0002] 1. Field of the Invention

[0003] The present invention relates to a correcting method of lens distortion of an image.

[0004] 2. Description of Related Art

[0005] An image taken by a video camera or a digital camera causes distortion (distortion aberration) depending upon a camera lens. The distortion aberration is generally extreme with a wide angle, increases toward a periphery of the image, and increases to 3% or more particularly with a zoom lens. Moreover, since the image is deformed in a so-called barrel or bobbin shape, this distortion aberration disadvantageously adversely affects estimation in computer vision, measurement using the image, and the like.

[0006] Parameters such as a distortion coefficient for determining a camera lens distortion are called a camera internal parameter, and cameral calibration means that the camera internal parameter is obtained. When the camera internal parameter can be obtained through camera calibration, the image lens distortion can be corrected by an image processing using the internal parameter.

[0007] A conventional camera calibrating method is disclosed in many research papers.

[0008] In the aforementioned conventional method, a relation between a certain point in a three-dimensional space and a point signal to the certain point (correspondence point) in the image is derived, and the internal parameter is estimated. There are two methods of obtaining the correspondence point as follows:

[0009] (1) A person designates the correspondence point with a mouse or the like.

[0010] (2) A grid intersection point, a polygon or polyhedron apex, or another characteristic point is automatically detected.

[0011] Since either one of these methods includes a measurement error, a large number of correspondence points have to be given in order to inhibit an error influence. That is, the more the number of correspondence points is, the more precise the obtained internal parameter becomes. Therefore, correspondence of an enormous number of (several tens to several hundreds) points has to be designated.

[0012] Moreover, these methods have the following problems.

[0013] In the method (1) in which the person designates the correspondence point with the mouse or the like, the error is large, operation is monotonous, and much labor is necessary. Furthermore, fatigue increases the error during designation.

[0014] In the method (2) (characteristic point detecting method) of automatically detecting the characteristic point by the image processing, the number of points detectable by the processing is limited. For example, when the grid intersection point is detected, the number of intersection points is only 88 in total even with 8 longitudinal lines and 11 lateral lines. It is said that at least 200 points or more are necessary for obtaining the parameter with good precision, and it is necessary to prepare the corresponding grid, polyhedron, and the like.

Summary of the Invention

[0015] The present invention has been developed to solve the aforementioned various problems. That is, an object of the present invention is to provide an image lens distortion correcting method in which (1) no correspondence point needs to be given, (2) any special tool or object is not used, and (3) all points on an image are used, so that an internal parameter can automatically be obtained with high precision, and an image lens distortion can be corrected with high precision.

[0016] According to the present invention, there is provided an image lens distortion correcting method comprising: an image printing step (A) for printing an arbitrary image I₁in a computer; an image pickup step (B) for taking the printed image I₁with a camera having a lens distortion and obtaining a photographed image I₂in the computer; a parameter estimating step (C) for obtaining a parameter θ such that the image I₁is artificially distorted with the parameter θ and an obtained distorted image I₁^udagrees with the photographed image I₂; and an image correcting step (D) for using the obtained parameter θ to correct the image taken by the camera.

[0017] According to the aforementioned method of the present invention, first the arbitrary image I₁(Computer graphics, a picture extracted with a scanner, or an image taken by a digital camera) is prepared as a calibration pattern in the computer. Subsequently, the image I₁is printed with a printer having no distortion. The printed image I₁is taken by a camera to be corrected, and the distorted photographed image I₂is obtained. When the certain internal parameter θ is given, the image I₁can artificially be distorted. When the artificially distorted image I₁^udis compared with the distorted image I₂and the images accurately equal with each other, the given parameter θ is desirable, and the obtained parameter θ can be used to correct the camera lens distortion.

[0018] According to a preferred embodiment of the present invention, the parameter θ includes a position correction parameter θ^ufor correcting a position, and a distortion correction parameter θ^dfor correcting the distortion.

[0019] The position correction parameter θ^uis a parameter for conversion to an image I₁^uwith the corrected position from the image I₁. The conversion parameter is obtained in a least-squares method such that a difference r=I₁(p)−I₁^u(p+u) between a luminance value I₁(p) of a point p in the image I₁and a luminance value I₁^u(p+u) of a point p+u in the image I₁^ucorresponding to the point p is minimized entirely on the image. The obtained parameter is used as the position correction parameter θ^u.

[0020] Moreover, the distortion correction parameter θ^dis a parameter for conversion to the image I_i^ufrom the distorted image I₁^ud. The conversion parameter is obtained in the least-squares method such that a difference r=I₂(p)−I₁^u(f(p)) between a luminance value I₂(p) of the point p in the image I₂and a luminance value I₁^u(f(p)) of a point f(p) in the image I₁^ucorresponding to the point p is minimized entirely on the image. The obtained parameter is used as the distortion correction parameter θ^d.

[0021] Furthermore, in the image correcting step (D), the obtained distortion correction parameter θ^dis used to correct the image taken by the camera.

[0022] That is, as shown in FIG. 1, a position in the photographed image I₂in which the original image I₁is projected is not known. Therefore, two stages are performed in order to obtain the distorted image I₁^ud. That is, first the position corrected image I₁^uis obtained from the image I₁, subsequently the image I₁^uis distorted and the distorted image I₁^udis obtained.

[0023] First Stage

[0024] First, the parameter θ^ufor conversion to the image I₁^ufrom the image I₁is obtained. This does not include any internal parameter. The point in I₁^ucorresponding to the point p in the image I₁deviates by u and is represented by p+u. Here, u changes with p and θ^u. When the image I₁exactly equals with the photographed image I₂, the luminance value I₁(p) in p has to be equal to the luminance value I₂(p+u) in p+u. That is, the luminance difference r=I₁(p)−I₂(p+u) in each point has to be 0. An evaluation function Σr²is obtained by summing squared differences for all the points, and the parameter θ^uis automatically obtained by repeating calculations such that the evaluation function is minimized.

[0025] Second Stage

[0026] Subsequently, the parameter θ^dfor conversion to the image I₁^ufrom the distorted image I¹^udis obtained. This is the internal parameter such as a distortion coefficient. This stage is similar to the first stage, but instead of the conversion to the distorted image I₁^udfrom the image I₁^u, a reverse conversion to the image I₁^ufrom the distorted image I₁^udis considered.

[0027] The point in image I₁^ucorresponding to the point p in the photographed image I₂is represented by f(p). Here, f( ) changes with p and θ^d.

[0028] When the photographed image I₂exactly equals with the image I₁^u, the luminance value I₂(p) in p has to be equal to the luminance value I₁^u(f(p)) in f(p). That is, the luminance difference r=I₂(p)−I₁^u(f(p)) in each point has to be 0. The evaluation function Σr²is obtained by summing squared differences for all the points, and the parameter θ^uis automatically obtained by repeating calculations such that the evaluation function is minimized.

[0029] Correction

[0030] The internal parameter θ^dused for obtaining the distorted image I₁^udfrom the image I₁can be used to conversely obtain I₁from I₁^ud. That is, a processing for correcting the applied distortion can be performed.

[0031] When the image I₂taken by the camera is subjected to the same processing as described above, the image with the corrected distortion can be obtained.

[0032] In the parameter estimating step (C), it is preferable to alternately repeatedly use the position correction parameter θ^uand the distortion correction parameter θ^dand obtain the parameter θ.

[0033] Moreover, the parameter θ is a parameter for conversion to the image I₁^udfrom the image I₁. The conversion parameter is obtained in the least-squares method such that a difference r=I₁(p)−I₂(d(p+u(p)) between a luminance value I₁(p) of the point p in the image I₁and a luminance value I₂(d(p+u(p)) of a point d(p+u(p)) in the image I₂corresponding to the point p is minimized entirely on the image. The obtained parameter is used as the parameter θ.

[0034] According to the method, the parameter θ can be estimated by a single step.

[0035] Other objects and advantageous characteristics of the present invention will be apparent from the following description with reference to the accompanying drawings.

Brief Description of the Drawings

[0036] FIG. 1 is a diagram schematically showing a method of the present invention.

[0037] FIG. 2 shows a halftone image displayed on a display in a calibration pattern used in an embodiment.

[0038] FIG. 3 shows the halftone image displayed on the display before and after correction according to the method of the present invention.

[0039] FIG. 4 shows another halftone image displayed on the display before and after the correction according to the method of the present invention.

[0040] FIG. 5 is a schematic diagram of a method of repeating two steps.

[0041] FIG. 6 is a schematic diagram of an estimating method with one step.

Description of the Preferred Embodiments

[0042] A preferred embodiment of the present invention will be described hereinafter with reference to the drawings. Additionally, common portions in the respective drawings are denoted with the same reference numerals, and redundant description is omitted.

[0043] Camera calibration and lens distortion correction are important processes in computer vision. In recent years, self calibration has been studied. In many researches, to facilitate the calibration, distortion is not considered and respective problems are formulated. Therefore, it is necessary to calibrate an internal parameter of a camera beforehand.

[0044] Some pieces of calibration software have been utilized via Internet (e.g., Tsai's method). However, in such ordinary method, it is necessary to obtain a correspondence between a point on an image and a point on a known three-dimensional coordinate (a plane, or a cube, a model or another structure), and subsequently point conversion is estimated.

[0045] In this method, the point correspondence must manually be estimated. Therefore, an error is generated in a manual operation and this lacks in reliability. Additionally, the method requires much time and endurance. Therefore, it is almost impossible to measure a change of a distortion parameter during zoom-in/zoom-out of the camera.

[0046] One of alternative methods is detection of a marker. This can be performed by a template matching method to an extend of sub-pixels. However, there is another problem that the marker on the image corresponds to a coordinate in a known space. To improve estimation precision, it is necessary to increase the number of markers, and this problem cannot be ignored. Even if this problem can be avoided, the number of points for the correspondence is still limited.

[0047] In the present invention, there is proposed a new correcting method for correcting an image distortion by a lens zoom operation. In the proposed method, the correspondence is set, and the estimation is expected to be more precise than the marker detection. This is because all points of the image are used in the method, instead of using some points of the marker. This method is based on an image registration (superposition) method applied to a motion analysis region, and includes the following three steps, that is, affine transformation, plane projection transformation, and lens distortion restoration.

[0048] Three Estimating Steps by Registration

[0049] A basic idea is that a correspondence between points is necessary for calibration, and registration can satisfy this requirement. In the proposed method, an ideal calibration pattern is associated with a distorted image obtained by taking the printed pattern with the camera. Since a size of the pattern printed on a sheet can easily be measured, a three-dimensional coordinate of each pixel in the pattern can easily be determined. As other characteristics of the method, any image can be used as the calibration pattern. Moreover, once the parameter is estimated, the parameter can be used as long as the lens zoom operation is unchanged.

[0050] The method includes the following three steps. In a first step, the pattern is coarsely transformed into an image represented by an affine parameter and translatory parameter. In a second step, an accurate parameter of plane transformation under through-view projection is estimated by a nonlinear optimizing method. In a third step, a distortion parameter is estimated in order to minimize a residual which remains after application of the second step.

[0051] Affine Transformation by Detected Marker

[0052] First, a pattern image I₁is prepared. The pattern image has three-color (red, green, blue) markers having coordinates m_r, m_g, m_bin respective corners. Subsequently, the pattern image I₁is printed in color in order to facilitate detection of the marker.

[0053] Therefore, an image I₂of the printed pattern I₁is taken by a camera to be calibrated. Coordinates of markers m_r′, m_g′, m_b′ in the image I₂are detected by thresholding or template matching.

[0054] Here, a parameter for transforming the pattern image I₁to the pattern projected in the image I₂is calculated. The transformation is represented by six affine parameters θ^α=(θ^α₁, . . . θ^α₆)^T.

[0055] While p=(x, y)^Tis set as one point on the pattern I₁, p+a(p;θ^α) is set as one point on I₂corresponding to p. Here, equation (1) of [Equation 1] is established. Subsequently, a linear equation of equation (2) is solved, and parameter θ^α is obtained.

[0056] Equation 1

\begin{matrix} a (p; θ^{a}) = (\begin{matrix} x & y & 0 & 0 & 1 & 0 \\ 0 & 0 & x & y & 0 & 1 \end{matrix}) θ^{a} & (1) \end{matrix}

m′_i=m_i+α(m_i; θ^α),i=r, g, b (2)

[0057] Image Registration; Perspective Fitting

[0058] In the first step described above, only three correspondence points are employed with respect to the affine transformation. Subsequently, plane projection transformation is utilized, and accurate correspondence among individual points is set as image registration.

[0059] Model and Formulation

[0060] Here, a luminance residual between p_iin the image I₁and p₁+u(p_i; θ^u) in the image I₂corresponding to p_iis minimized. Here, when the residual is defined as in equations (3), (4) of [Equation 2], a function to be minimized is obtained as in equation (5).

[0061] Equation 2

[0062] r_i=I₁(p_i)−I₂(p_i+u(p; θ^u)) (3)

\begin{matrix} \begin{matrix} u (p; θ^{u}) = (\begin{matrix} x & y & 0 & 0 & 1 & 0 & x^{2} & x y \\ 0 & 0 & x & y & 0 & 1 & x y & y^{2} \end{matrix}) θ^{u} \\ = M^{u} θ^{u} \end{matrix} & (4) \\ \begin{matrix} \min_{θ^{u}} \sum_{i} {ρ (r_{i})}_{1} & ρ (r_{i}) = r_{i}^{2} \end{matrix} & (5) \end{matrix}

[0063] Here, u( ) denotes displacement of a point on a plane, when a point of view is changed from a certain point to another point under perspective projection. The displacement is a motion model on the plane under the perspective projection, and is represented by eight parameters θ^u=(θ₁^u, . . . , θ₈^u)^Tfrequently used in motion analysis.

[0064] Minimizing Method

[0065] The function (5) is minimized by Gauss-Newton method as a famous nonlinear optimization technique and the parameter θ^uis estimated. The parameter is updated by equation (6) of [Equation 3].

[0066] As an initial value of first six elements of the parameter θ^u, the affine parameter θ^α obtained in the first stage is used. Last two elements of the parameter θ^uare initialized/set to 0.

[0067] A direction δθ^uin which the function (5) is decreased is calculated by equations (7) to (10).

[0068] Equation 3

[0069] θ^u←θ^u+δθ^u (6)

[0070] δθ^u=−(J{tilde over (D)}J^T)⁻¹J{tilde over (D)}r (7)

\begin{matrix} J = J (θ^{u}) = \frac{\partial r}{\partial θ^{u}} = [\frac{\partial r_{i}}{\partial θ_{j}^{u}}] & (8) \\ \overline{D} = d i a g [\frac{\dot{ρ} (r_{i})}{r_{i}}] & (9) \\ \dot{ρ} (r_{i}) = \frac{\partial ρ (r)}{\partial r} |_{r = r_{i}} & (10) \end{matrix}

[0071] This is the same as formulation of the least-squares method. That is, the equation is the same as a simultaneous linear equation represented by equation (11) of [Equation 4] for k=1, . . . , 8. A differentiating chained law is used to convert a partial derivation as represented by equation (12). The calculation of δθ^uin the equation (6) is repeated until δθ^uconverges. In each repetition, the parameter estimated in the previous repetition is used in calculation of u(p). When the repeated calculation stops, the estimated parameter is θ^u′,

[0072] Equation 4

\begin{matrix} \sum_{l, i} \frac{\dot{ρ} (r_{i})}{r_{i}} \frac{\partial r_{i}}{\partial θ_{k}^{u}} \frac{\partial r_{i}}{\partial θ_{l}^{u}} {δθ}_{l}^{u} = - \sum_{i} \frac{\dot{ρ} (r_{i})}{r_{i}} r_{i} \frac{\partial r_{i}}{\partial θ_{k}^{u}} for k = 1, \dots, 8. & (11) \\ \frac{\partial r}{\partial θ^{u}} = \frac{\partial u}{\partial θ^{u}} \frac{\partial r}{\partial u} = - M^{u T} \nabla I_{2} (p + u (p)) & (12) \end{matrix}

[0073] Image Registration; Distortion Fitting

[0074] In the previous step, except the influence by the lens distortion, the image registration between the pattern image and the distorted image pattern ends.

[0075] Distortion Model

[0076] A relation between a coordinate having no distortion is and coordinate having the distortion in a certain image is usually modeled by the following five camera internal parameters, that is, distortion parameters k₁and k₂, image center coordinate (c_x, c_y)^Tand a ratio s_xof a pixel width to height. These parameters are set to θ^d=(k₁, k₂, c_x, c_y, s_x)^T.

[0077] The distortion is represented by a coordinate system having an origin in (c_x, c_y)^T. On the other hand, the previously used coordinate system has an origin in an uppermost left corner. Therefore, another notation is introduced as follows.

[0078] That is, P_u=(x_u, y_u)^Tand p_d=(x_d, u_d)^Tare set as points in an image having no distortion and image having the distortion, respectively, and both images have an origin in the left upper corner of the image. Moreover, (ξ_u, η_u)^Tand (ξ_d, η_d)^Tare set as points in the image having no distortion and image having the distortion, respectively, in which respective origins are in the image center (c_x, c_y)^T. These points can be represented by equations (13) to (15) of [Equation 5].

[0079] As described above, ξ_uis represented to be positive as the function of ξ_dbut this does not apply to ξ_d. In order to obtain ξ_dfrom ξ_u, equation (16) is repeatedly solved.

[0080] Equation 5

[0081] (ξ_u, η_u)^T=(x_u, y_u)^T−(c_x, c_y)^T (13)

[0082] (ξ_d, η_d)^T=(ξ_u, η_u)^T−(κ₁R²+κ₂R⁴)(ξ_d, η_d)^T (14)

[0083] (x_d, y_d)^T=(s_xξ_d, η_d)^T+(c_x, c_y)^T (15)

[0084] R={square root over (ξ

[0085]

\begin{matrix} {(ξ_{d k}, η_{d k})}^{T} = \frac{{(ξ_{u}, η_{u})}^{T}}{1 + κ_{1} R_{k - 1}^{2} + κ_{2} R_{k - 1}^{4}} & (16) \end{matrix}

R_k={square root over (ξ_dk²+η_dk²)},

[0086] R₀={square root over (ξ_u²+η_u²)}

[0087] The aforementioned relation is used to obtain two functions between p_uand p_dfor the coordinate system having the origin in the uppermost left corner as represented by equations (17) and (18) of [Equation 6]. Here, f and d are mutually inverse functions. However, since d corresponds to the equation (16), d is not a explicit function of p_u.

[0088] In any case, the equations (17) and (18) can be used to represent transformation between two images. While I₁^uis the image of the pattern transformed by applying θ^u′ to I₁, I₁^udis the image of the pattern transformed by applying θ^dto I₁^u. That is, the relation is represented by equations (20) to (22). These relations are shown in FIG. 1.

[0089] Equation 6

[0090] P_d=d(p_u; θ^d) (17)

[0091]

\begin{matrix} \begin{matrix} p_{u} = f (p_{d}; θ^{d}) \\ = (\begin{matrix} \frac{x_{d} - c_{r}}{s_{r}} (1 + κ_{1} R^{′2} + κ_{2} R^{′4}) + c_{x} \\ (y_{d} - c_{y}) (1 + κ_{1} R^{′2} + κ_{2} R^{′4}) + c_{y} \end{matrix}) \end{matrix} & (18) \\ R^{'} = \sqrt{{(\frac{x_{d} - c_{r}}{s_{r}})}^{2} + {(y_{d} - c_{y})}^{2}} & (19) \end{matrix}

I_i(p)=I_i^u(p+u(p; θ^u)) (20)

[0092] I₁^u(p)=I₁^ud(d(p)) (21)

[0093] I_i^u(f(p))=i₁^ud(p) (22)

[0094] Minimization by Reverse Registration

[0095] The transformation to the image I₁^udfrom I₁is represented by u and d. However, since d is not the explicit function, the parameter between I₁and I₁^udcannot directly be estimated by registration.

[0096] Here, reverse registration is considered. As shown in FIG. 1, the image I₂is matched with I₁^ud. That is, the residual of luminance between two images represented by equation (23) of [Equation 7] is minimized.

[0097] However, this equation can be represented by equation (24) using equation (22). Then, an estimating method is the same as that of the previous step. The minimization is performed on the function of equation (25) by the Gauss-Newton method.

[0098] Here, Ω={i;p_iεI₂, ∃pεI₁, f(p_i)=p+u(p)} means that the minimization has to be applied to the point inside the region corresponding to the pattern I₁among points in I₂.

[0099] Equation 7

[0100] r_i=I₁^ud(p_i)−I₂(p_i) (23)

[0101] r₁=I₂(p_i)−I₁^u(f(p_i; θ^d) (24)

[0102]

\begin{matrix} \min_{θ^{d}} \sum_{Ω} ρ (r_{i}) Ω = {i; p_{i} ε I_{2}, \exists p ε I_{1}, f (p_{i}) = p + u (p)} & (25) \end{matrix}

[0103] The simultaneous linear equation in which solution needs to be found has the same form as that of the equation (11), and is represented by equation (26) of [Equation 8]. The derivation in the equation (26) is represented by equation (27).

[0104] Equation 8

\begin{matrix} \sum_{l, i εΩ} \frac{\dot{ρ} (r_{i})}{r_{i}} \frac{\partial r_{i}}{\partial θ_{k}^{d}} \frac{\partial r_{i}}{\partial θ_{l}^{d}} {δθ}_{l}^{d} = - \sum_{i εΩ} \frac{\dot{ρ} (r_{i})}{r_{i}} r_{i} \frac{\partial r_{i}}{\partial θ_{k}^{d}} & (26) \\ \frac{\partial r}{\partial θ^{d}} = \frac{\partial f}{\partial θ^{d}} \frac{\partial r}{\partial f} = \frac{\partial f}{\partial θ^{d}} (- \nabla I_{1}^{u} (f (p))) & (27) \end{matrix}

[0105] Jacobian matrix of the equation (18) is represented by equations (28) to (34) of [Equation 9].

[0106] Initial parameters c^x, c_y, k₂and s_xfor solving the equation (26) are set to ½ of width and height of I₂, 0 and 1, respectively. On the other hand, in all of equations (29) to (32), k₁is arbitrarily initially set to avoid 0 by initial setting of k₁=k₂=0. Here, k₁∈ [−10⁻⁷, 10⁻⁷] is empirically selected.

[0107] Equation 9

\begin{matrix} \frac{\partial f (p_{d})}{\partial θ^{d}} = (\begin{matrix} R^{′2} x_{d} - \frac{c_{r}}{s_{r}} & R^{′2} (y_{d} - c_{y}) \\ R^{′4} x_{d} - \frac{c_{r}}{s_{r}} & R^{′4} (y_{d} - c_{y}) \\ \frac{\partial x_{y}}{\partial c_{r}} & \frac{\partial y_{u}}{\partial c_{r}} \\ \frac{\partial x_{u}}{\partial c_{y}} & \frac{\partial y_{u}}{\partial c_{y}} \\ \frac{\partial x_{u}}{\partial s_{r}} & \frac{\partial y_{u}}{\partial s_{r}} \end{matrix}) & (28) \\ \begin{matrix} \frac{\partial x_{u}}{\partial c_{x}} = 1 - \frac{1}{s_{x}} (1 + κ_{1} R^{′2} + κ_{2} R^{′4}) \\ - 2 (κ_{1} + 2 κ_{2} R^{′2}) \frac{{(x_{d} - c_{x})}^{2}}{s_{x}^{3}} \end{matrix} & (29) \\ \frac{\partial y_{u}}{\partial c_{x}} = - 2 (κ_{1} + 2 κ_{2} R^{′2}) \frac{x_{d} - c_{x}}{s_{x}^{2}} (y_{d} - c_{y}) & (30) \\ \frac{\partial x_{u}}{\partial c_{y}} = - 2 (κ_{1} + 2 κ_{2} R^{′2}) \frac{x_{d} - c_{x}}{s_{x}} (y_{d} - c_{y}) & (31) \\ \begin{matrix} \frac{\partial y_{u}}{\partial c_{y}} = 1 - (1 + κ_{1} R^{′2} + κ_{2} R^{′4}) \\ - 2 {(y_{d} - c_{y})}^{2} (κ_{1} + 2 κ_{2} R^{′2}) \end{matrix} & (32) \\ \begin{matrix} \frac{\partial x_{u}}{\partial s_{x}} = \frac{- (x_{d} - c_{x})}{s_{x}^{2}} (1 + κ_{1} R^{′2} + κ_{2} R^{′4}) \\ - 2 (κ_{1} + 2 κ_{2} R^{′2}) \frac{{(x_{d} - c_{x})}^{3}}{s_{x}^{4}} \end{matrix} & (33) \\ \frac{\partial y_{u}}{\partial s_{x}} = - 2 (y_{d} - c_{y}) (κ_{1} + 2 κ_{2} R^{′2}) \frac{{(x_{d} - c_{x})}^{2}}{s_{x}^{3}} & (34) \end{matrix}

[0108] Repetition of Two Steps

[0109] FIG. 5 is a schematic diagram showing a method of repeating two steps. Even when there is a relatively large residual in an initial state, a strategy for shifting to a fine resolution from a coarse resolution (coarse-to-fine) is employed in order to shorten a calculation time and perform precise estimation. That is, first a filtered and blurred image is processed, a filter is gradually changed to a filter having a smaller blur, and a cleared image is processed. Therefore, while the resolutions of I₁, I₁^u, I₂are changed to the fine resolution from the coarse resolution, second and third steps are alternately repeated (see FIG. 5).

[0110] Therefore, the estimation of the second step has to be corrected to include an estimated value θ^dof the third step.

[0111] Here, I₂^fis assumed to be an image obtained by correcting distortion of the photographed image I₂through transformation f using the parameter θ^d. That is, the image can be represented by equations (35), (36).

[0112] This is used to correct the equation (3) of the residual as represented by equation (37).

[0113] In this case, Jacobian equation (12) turns to equation (38).

[0114] The corrected second step, and the third step are alternately repeated until the estimated value converges.

[0115] The image having no distortion can be obtained from the photographed image by the transformation f using the parameter θ^dout of the obtained parameters.

[0116] Equation 10

[0117] I₂^𝑓(𝑓(p^{; θd}); θ^d)=I₂(p) (35)

[0118] I₂^f(p^{; θd})=I₂(d(p^{; θd})) (36)

[0119] r_i=I₁(p_i)−I₂^𝑓(p_i+u(p_i; θu); θ^d) (37)

[0120]

\begin{matrix} \frac{\partial r}{\partial θ^{u}} = - M^{u T} \nabla I_{2}^{f} (p + u (p; θ^{u}); θ^{d}) & (38) \end{matrix}

[0121] Estimation with One Step

[0122] FIG. 6 is a schematic diagram showing an estimating method with one step. The estimation by the repetition of two steps has been described above, but these two steps can be unified into one estimating step.

[0123] Since d is not the explicit function, that is, d is an implicit function unable to be represented by any equation, a function obtained by synthesizing u and d also becomes implicit function. This synthesized function u·d transforms the image I₁to I₁^ud. Therefore, the parameter θ^u, θ^dmay be obtained such that I^udaccurately equals with I₂. In this case, Jacobian of the synthesized function needs to be explicit, and the explicit function is obtained from Jacobian of u and d.

[0124] In actual, since d is an implicit function, Jacobian cannot be obtained by directly differentiating the equation. However, when the implicit function theorem known in a mathematical field is used, Jacobian of d can be obtained as the explicit function.

[0125] Here, it is assumed that p_uand θ^dare unified into one, and represented by q=(p_u, θ^d) Then, function F of equation (39) as the function of q is considered.

[0126] In this case, d is an implicit function defined by p_d=d(q). Therefore, Jacobian of d can be calculated by the implicit function theorem as represented by equations (40) and (41).

[0127] Here, E₂is a 2×2 unit matrix, and represented by equation (41a). Moreover, ∂f/∂θ^dis the same as the equation (28). Moreover, ∂f/∂p_dcan be calculated from the equation (18), and is actually represented by equations (42) and (43) using the equations (29) to (32).

[0128] Equation 11

[0129] F(q,p_d)=p_u−𝑓(p_d, θ^d)=0 (39)

[0130]

\begin{matrix} \frac{\partial d}{\partial q} = - \frac{\partial F^{- 1}}{\partial p_{d}} \frac{\partial F}{\partial q} & (40) \\ = - \frac{\partial f^{- 1}}{\partial p_{d}} (- E_{2} \frac{\partial f}{\partial θ^{d}}) & (41) \\ E_{2} = (\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}) & (41a) \\ \frac{\partial f}{\partial p_{d}} = (\begin{matrix} \frac{\partial x_{u}}{\partial x_{d}} & \frac{\partial y_{u}}{\partial x_{d}} \\ \frac{\partial x_{u}}{\partial y_{d}} & \frac{\partial y_{u}}{\partial y_{d}} \end{matrix}) & (42) \\ = (\begin{matrix} 1 - \frac{\partial x_{u}}{\partial c_{x}} & - \frac{\partial y_{u}}{\partial c_{x}} \\ - \frac{\partial x_{u}}{\partial c_{y}} & 1 - \frac{\partial y_{u}}{\partial c_{y}} \end{matrix}) & (43) \end{matrix}

[0131] As described above, Jacobian of the function obtained by synthesizing transformations u, d can be calculated. A squared sum of luminance differences from the point in the photographed image I₂corresponding to the point p_iin the pattern I₁, and is represented by equations (44), (45).

[0132] Parameter θ={θ^u, θ^d} for minimizing the squared sum is estimated (see FIG. 6). The optimizing method for the estimation is performed similarly as the estimation of the second and third steps. That is, the simultaneous linear equation to be solved is represented by equation (46).

[0133] Here, partial differential ∂r_i/∂θ_kis derived as an element of Jacobian of equations (47) to (52).

[0134] Here, p_d=d(p_u), p_u=p+u (p) Moreover, E₅is a 5×5 unit matrix.

[0135] The image having no distortion can be obtained from the photographed image by the transformation f using the parameter θ^dout of the obtained parameters.

[0136] Equation 12

[0137] r_i=I₁(p_i)−I₂(d(p_i+u(p_i^{, θu})), θ^d) (44)

[0138]

\begin{matrix} \min_{θ^{u}} \sum_{i} ρ (r_{i}), ρ (r_{i}) = r_{i}^{2} & (45) \\ \sum_{l} \sum_{i} \frac{\dot{ρ} (r_{i})}{r_{i}} \frac{\partial r_{i}}{\partial θ_{k}} \frac{\partial r_{i}}{\partial θ_{l}} {δθ}_{l} = - \sum_{i} \frac{\dot{ρ} (r_{i})}{r_{i}} r_{i} \frac{\partial r_{i}}{\partial θ_{k}} & (46) \\ \frac{\partial r}{\partial θ} = \frac{\partial r}{\partial q} \frac{\partial q}{\partial θ} & (47) \\ = \frac{\partial r}{\partial d} \frac{\partial d}{\partial q} (\frac{\partial q}{\partial θ^{u}} \frac{\partial q}{\partial θ^{d}}) & (48) \\ = \frac{\partial r}{\partial d} \frac{\partial d}{\partial q} (\begin{matrix} \frac{\partial p_{u}}{\partial θ^{u}} & \frac{\partial p_{u}}{\partial θ^{d}} \\ \frac{\partial θ^{d}}{\partial θ^{u}} & \frac{\partial θ^{d}}{\partial θ^{d}} \end{matrix}) & (49) \\ = \frac{\partial r}{\partial d} \frac{\partial d}{\partial q} (\begin{matrix} \frac{\partial p_{u}}{\partial θ^{u}} & 0 \\ 0 & E_{5} \end{matrix}) & (50) \\ = - \nabla I_{2} (d (p_{u})) \frac{\partial d}{\partial q} (\begin{matrix} M^{u} & 0 \\ 0 & E_{5} \end{matrix}) & (51) \\ = \nabla I_{2} (d (p_{u})) \frac{\partial f^{- 1}}{\partial p_{d}} (\begin{matrix} - M^{u} & \frac{\partial f}{\partial θ^{d}} \end{matrix}) & (52) \end{matrix}

[0139] Several Strategies

[0140] (1) Interpolation of Pixel Value

[0141] When the luminance of the point having no coordinate on an integer grid is necessary, bilinear interpolation is applied between adjacent pixel values.

[0142] (2) Histogram Matching

[0143] For the image taken by the camera, since pattern luminance often changes, image histogram is converted so that the histogram of the image becomes the same as that of the pattern.

Examples

[0144] While the camera zoom parameter was changed, an experiment was performed for the actual image taken by the camera in the proposed method. A picture shown in FIG. 2 was used as the calibration pattern printed by a laser printer (EPSON LP-9200PS2), and an image of the picture was taken by a CCD camera (Sony EVI-D30) fixed on a support with image taking software (SGI O2). While the camera zoom was changed from a widest angle, two images of grid and calibration patterns were photographed. The wider the angle is, the more the distortion influence increases. Then, a grid line is curved.

[0145] An estimation result of the internal parameter according to the present invention is shown in Table 1. Moreover, this internal parameter is used to show images before and after the correction in FIG. 3 and FIG. 4. It is seen that the curved line in the grid pattern is transformed into a straight line by the correction.

	TABLE 1

	(a)	(e)

[0146] A gradation change of an image illuminance cannot be removed by a simple histogram conversion which is to be replaced with the estimation of the illuminance change by some methods such as linear luminance restriction. Therefore, the grid line is still slightly curved particularly in a peripheral portion of the image, and as shown on the right side of FIG. 3, there is a slight error in the experiment result.

[0147] In an experiment simulation in which a noise was added to each pixel of the transformed pattern and the pattern was used as the distorted image, the proposed method satisfactorily functioned even with a uniform amplitude of added noise of ±50 or more.

[0148] As described above, in the present invention, there is proposed an inventive method for obtaining the internal parameter of the camera to correct the image distortion. The proposed method is based on the image registration, and includes two nonlinear optimizing steps of perspective fitting and distortion fitting by geometric transformation. The experiment result has proved an efficiency of the present proposed method which can reduce an operator's labor. The nonlinear optimization requires time, but the technique is sufficiently effective as a batch process.

[0149] As described above, there is proposed an inventive image lens distortion correcting method for obtaining the camera internal parameter to correct the distortion of the image taken by the camera with the zoom lens attached thereto. This method is based on image registration. First, a certain calibration pattern is transformed into the distorted image pattern using affine transformation. Subsequently, the nonlinear optimization of Gauss-Newton method for minimizing the luminance difference of pixel values of two images is applied, and the two images are registered. Finally, the distortion parameter is estimated to minimize the luminance difference remaining after the second step. The experiment result shows effectiveness of the proposed method.

[0150] Therefore, in the image lens distortion correcting method of the present invention, (1) no correspondence point needs to be given, (2) any special tool or object is not used, and (3) all points on the image are used, so that the internal parameter can automatically be obtained with high precision, the image lens distortion can be corrected with high precision, and other superior effects are produced.

[0151] Additionally, the present invention is not limited to the aforementioned embodiment, and can of course variously be changed within the scope of the present invention.

Claims

What is claimed is:

1. An image lens distortion correcting method comprising:

an image printing step (A) for printing a prepared calibration pattern image I₁in a computer;

an image pickup step (B) for photographing the printed image I₁with a camera having a lens distortion and obtaining a photographed image I₂in the computer;

a parameter estimating step (C) for obtaining a parameter θ such that the image I₁is artificially distorted with the parameter θ and so as to obtain a distorted image I₁^udequal with the photographed image I₂; and

an image correcting step (D) for using the parameter θ to correct the photographed image I₂.

2. The image lens distortion correcting method according to claim 1 wherein said parameter θ comprises a position correction parameter θ^ufor correcting a position, and a distortion correction parameter θ^dfor correcting the distortion.

3. The image lens distortion correcting method according to claim 2, wherein said position correction parameter θ^uis a first parameter for conversion to an image I₁^uwith the corrected position from the image I₁, the first conversion parameter is calculate in a least-squares method such that a difference r=I₁(p)−I₁^u(p+u) between a luminance value I_i(p) of a point p in the image I₁and a luminance value I₁^u(p+u) of a point p+u in the image I₁^ucorresponding to the point p is minimized entirely on the image I₁, and the first parameter thus calculated is used as the position correction parameter θ^u.

4. The image lens distortion correcting method according to claim 3, wherein said distortion correction parameter θ^dis a second parameter for conversion to the image I₁^ufrom the distorted image I₁^ud, the second conversion parameter is calculated in the least-squares method such that a difference r=I₂(p)−I₁^u(f(p)) between a luminance value I₂(p) of the point p in the image I₂and a luminance value I₁^u(f(p)) of a point f(p) in the image I₁^ucorresponding to the point p is minimized entirely on the image I₂, and the second parameter thus calculated is used as the distortion correction parameter θ^d.

5. The image lens distortion correcting method according to claim 4, wherein said parameter estimating step (C) comprises a step of alternately repeatedly using said position correction parameter θ_uand said distortion correction parameter θ^dto calculate the parameter θ.

6. The image lens distortion correcting method according to claim 1, wherein said parameter θ is a parameter for conversion to the image I₁^udfrom the image I₁, the conversion parameter is calculated in a least-squares method by using an implicit function theorem such that a difference r=I₁(p)−I₂(d(p+u(p)) between a luminance value I₁(p) of a point p in the image I₁and a luminance value I₂(d(p+u(p)) of a point d(p+u(p)) in the image I₂corresponding to the point p is minimized entirely on the image I₁, and the conversion parameter thus calculated is used as the parameter θ.

7. The image lens distortion correcting method according to claim 4, wherein said image correcting step (D) comprises a step of using the calculated distortion correction parameter θ^dto correct the photographed image I₂.

8. An image lens distortion correcting method comprising the steps of:

(a) printing a prepared calibration pattern image I₁in a computer;

(b) photographing the printed image I₁with a camera having a lens distortion to obtain a distorted photographed image I₂in the computer;

(c) calculating a parameter θ such that the image I₁is artificially distorted with the parameter θ so as to obtain a distorted image I₁^udequal with the photographed image I₂; and

(d) using the calculated parameter θ to correct the photographed image I₂.

9. An image lens distortion correcting method according to claim 8, wherein said parameter θ comprises a position correction parameter θ^ufor correcting a position, and a distortion correction parameter θ_dfor correcting distortion.

10. An image lens distortion correcting method according to claim 9, wherein said position correction parameter θ^uis a first parameter for conversion to an image I₁^uwith the corrected position from the image I₁, the first conversion parameter is calculated in a least-squares method such that a difference r=I₁(p)−I₁^u(p+u) between a luminance value I₁(p) of a point p in the image I₁and a luminance value I₁^u(p+u) of a point p+u in the image I₁^u, corresponding to the point p, is minimized entirely on the image I₁, and the first conversion parameter thus calculated is used as the position correction parameter θ^u.

11. An image lens distortion correcting method according to claim 10, wherein said distortion correction parameter θ^dis a second parameter for conversion to the image I₁^ufrom the distorted image I₁^ud, the second conversion parameter is calculated in the least-squares method such that a difference r=I₂(p)−I₁^u(f(p)) between a luminance value I₂(p) of the point p in the image I₂and a luminance value I₁^u(f(p)) of a point f(p) in the image I₁^u, corresponding to the point p, is minimized entirely on the image I₂, and the second conversion parameter thus calculated is used as the distortion correction parameter θ^d.

12. An image lens distortion correcting method according to claim 11, wherein calculation of the parameter θ comprises alternately repeatedly using said position correction parameter θ^uand said distortion correction parameter θ^dto calculate the parameter θ.

13. An image lens distortion correcting method according to claim 8, wherein said parameter θ is a parameter for conversion to the image I₁^udfrom the image I₁that is calculated in a least-squares method by using an implicit function theorem such that a difference r=I₁(p)−₂(d(p+u(p)) between a luminance value I₁(p) of a point p in the image I₁and a luminance value I₂(d(p+u(p)) of a point d(p+u(p)) in the image I₂, corresponding to the point p, is minimized entirely on the image I₁, and the parameter for conversion to the image I₁^udfrom the image I₁thus calculated is used as the parameter θ.

14. The image lens distortion correcting method according to claim 11, wherein the calculated distortion correction parameter θ^dis used to correct the photographed image I₂.

15. The image lens distortion correcting method according to claim 9, wherein said distortion correction parameter θ^dis a five parameter distortion comprising k₁, k₂, c_x, c_y, and s_x, where k₁and k₂are distortion parameters, c_xand c_yrepresent an image center coordinate, and s_xis a ratio of pixel width to pixel height.

16. The image lens distortion correcting method according to claim 15, wherein said position correction parameter θ^uis an eight parameter displacement comprising θ₁^u, θ₂^u, θ₃^u, θ₄^u, θ₅^u, θ₆^u, θ₇^u, and θ₈^u.

Patent Citations (12)

Patent	Date	Inventor	Cited By
US4299482(A)	1981-11-01	Task	Examiner
US5406375(A)	1995-04-01	Brandstetter	Examiner
US5555313(A)	1996-09-01	Zheng et al.	Examiner
US6281931(B1)	2001-08-01	Tsao et al.	Examiner
US6437823(B1)	2002-08-01	Zhang	Examiner
US6573986(B2)	2003-06-01	Smith et al.	Examiner
US6618494(B1)	2003-09-01	Nonay et al.	Examiner
JP2-236433	1990-09-01		Examiner
JP8-201021	1996-08-01		Examiner
JP10-253499	1998-09-01		Examiner
JP11-237216	1999-08-01		Examiner
JPP2000-3437(A)	2000-01-01		Examiner

Non-Patent Literature (5)

Pedrick, G., A First Course in Analysis, pp. 160-168, 1994, Springer-Verlag New York, Inc.Applicant
Protter, M. H. et al., A First Course in Real Analysis, pp. 332-364, 1977, Springer-Verlag New York, Inc.Applicant
Tamaki, T. et al., “Unified Approach to Image Distortion.”, Not a Prior Art, 2002.Applicant
Microsoft Research, Microsoft Corporation, Technical Report, MSR-TR-97-71, Dec. 2, 1998 (Last Update: Mar. 25, 1999), “A Flexible New Technique for Camera Calibration”, Z. Zhang.Applicant
IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 18, No. 8, Aug. 1996, “Compact Representations of Videos Through Dominant and Multiple Motion Estimation”, S. Sawhney et al., pp. 814-830.Applicant