Full-model lens distortion correction method using dual-feature joint measurement

Abstract

The present invention proposes a full-model lens distortion correction method using dual-characteristic joint measurement, and the steps are as follows: perform sub-pixel level point coordinate detection on the detected image, and determine the distortion center using the least square method according to the relationship equation between the basic matrix and the distortion center; calculate the inclination angle and transformation factor of the pixel unit on the detected image, constitute the inverse mapping matrix of the discrete imaging process of the digital image, and convert the non-standard pixel unit into a square pixel unit; calculate the coordinates of the vanishing point according to the significant characteristics of the perspective projection and the homography matrix, and use the vanishing point to update the fitted distortion line; on the basis of the full distortion model, establish a function constraint relationship with the distortion line according to the geometric constraints of the vanishing point, and use the optimization method to minimize the straight line residual to obtain the distortion parameter, inclination angle and transformation factor, so as to realize accurate distortion correction of the camera. The present invention can realize lens distortion correction quickly and accurately, and has good correction effect and strong robustness.

Description

Technical Field

[0001] The present invention relates to the technical field of close-range measurement and digital image processing, and in particular to the technical field of lens distortion correction. The present invention relates to a full-model lens distortion correction method using dual-characteristic joint measurement, and uses a chessboard pattern to achieve accurate lens distortion correction. Background Art

[0002] In recent years, with the rapid development of artificial intelligence, more and more researchers have devoted themselves to the research in this field. Since machine vision plays an important role in artificial intelligence, zoom lenses are widely used in aerial photography, tracking and measurement. Due to the imperfection of lenses and the non-alignment of optical elements, lens distortion is prevalent in optical systems, especially in fisheye and wide-angle lenses. When the lens has non-negligible distortion, the use of an ideal pinhole model may lead to high measurement errors, making the pinhole model no longer applicable. Therefore, lens distortion must be accurately eliminated to maintain a linear imaging process in projective geometry.

[0003] According to the principles of projective geometry, if the camera lens is distortion-free, the straight line in space is still a straight line on the image plane, the cross ratio of the straight line remains unchanged under the projective transformation, the parallelism of the straight lines remains unchanged, and the area ratio of parallel figures and the length ratio of parallel line segments remain unchanged. In fact, after obtaining an image containing distortion, the straight line features in the scene will be bent into curves after being projected into the image due to the influence of lens distortion. The greater the nonlinear distortion of the lens, the greater the degree of curvature of the curve after the spatial straight line is projected.

[0004] In recent years, the research on lenses has mainly focused on the use of point correspondence and projection invariants. Weng et al. used the correspondence between the three-dimensional world coordinates of the object and the image coordinates to correct radial, eccentric and thin prism distortions. Gao and Yin proposed a method to estimate the ideal homography matrix to calculate the ideal image point representing the actual position of the distorted image point. Zhang used point correspondence at different angles to estimate radial distortion and eccentric distortion parameters. However, all three methods require the calculation of point correspondences.

[0005] As the method of removing unnecessary restrictions has strong advantages and is easy to construct, line-based methods have been widely used in distortion correction. Prescott and McLean proposed a radial lens distortion correction method that detects line information, but encountered local minima during the optimization process and only considered the influence of radial distortion. Zhou proposed a line-based camera calibration method that corrects the distortion of a single image by using three squares of unknown length, but this method has high requirements on the straight line segments in the image and is more suitable for processing images in special scenes. They only focus on radial distortion rather than complete distortion.

[0006] At present, traditional lens distortion correction methods have many shortcomings, such as: (1) the image center or principal point coordinates are taken as the distortion center by default, and the distortion center is not accurately calculated. As a result, it is impossible to use distortion parameter compensation to correct it, nor is it possible to use external orientation parameters to correct it; (2) the parallelogram pixels existing in the projective transformation are not considered, and they are not processed, and the linear imaging process in projective geometry cannot be guaranteed; (3) the algorithm based on the straight line theory only focuses on radial distortion or eccentric distortion, and does not perform lens distortion correction under the complete distortion model. Summary of the invention

[0007] In view of the technical problems that traditional lens distortion correction methods do not accurately calculate the distortion center, do not consider non-square pixels, and do not handle lens distortion under the full distortion model, the present invention proposes a full-model lens distortion correction method using dual-characteristic joint measurement, which can quickly realize accurate correction of lens distortion with high accuracy and strong robustness, and can be applied in solving specific practical problems related to image analysis and processing.

[0008] In order to achieve the above object, the technical solution of the present invention is implemented as follows: a full model lens distortion correction method using dual-characteristic joint measurement, the steps of which are as follows:

[0009] Step 1: Perform sub-pixel point coordinate detection on the detected image, and determine the distortion center using the least squares method based on the relationship equation between the basic matrix and the distortion center;

[0010] Step 2: Calculate the tilt angle and transformation factor of the parallelogram pixel of the point coordinates on the detected image, form the inverse mapping matrix of the discrete imaging process of the digital image, and transform the non-standard pixel unit into a square pixel unit;

[0011] Step 3: Calculate the coordinates of the vanishing points based on the salient features of the perspective projection and the homography matrix, and use the vanishing points to update the fitted distortion line;

[0012] Step 4: Based on the full distortion model, a functional constraint relationship with the distortion line is established according to the geometric constraints of the vanishing point. The optimization method is used to minimize the straight line residual to obtain the distortion parameters, tilt angle and transformation factor, so as to achieve accurate distortion correction of the camera.

[0013] The detected image is obtained by collecting the chessboard pattern through a CCD camera.

[0014] In the step 1, the detected image is grayed to obtain a grayscale image, and the grayscale image is subjected to sub-pixel point coordinate detection, i.e., the coordinates of the square vertices in the chessboard image are detected to form a point set E, and the coordinates of the distortion center are calculated using the point set E and the equation of the distortion center.

[0015] The method for calculating the coordinates of the distortion center is: Let x w is the known coordinate point on the chessboard pattern corresponding to the detected image, x u is the image point after the distortion model correction, x d is the corresponding distortion point in the pixel coordinate system on the chessboard image; due to the influence of the distortion center c, the undistorted image point x u Move to distortion point x d For: x d =c+k(x u -c); left multiply [x d ] T [c] x , we get the equation for the distortion center:

[0016] k[x d ] T [c] x Hx w =[x d ] T [c] x Hx w =0;

[0017] Where k is the radial distortion coefficient, [c] x represents an antisymmetric 3×3 matrix, the basic matrix F = [c] x H; H is the transformation matrix and x u =Hx w , the distortion center is extracted as the left pole of the basic matrix F;

[0018] Using the least squares method, we get equation c T f j = 0, j = 1, ..., 3, f j is the jth column of the basic matrix F; the least squares solution of all equations is the coordinate point of the distortion center c.

[0019] The implementation method of converting the non-standard pixel unit into a square pixel unit in step 2 is:

[0020] Assume that the length and width of the parallelogram pixel are d x and d y , one side of the parallelogram pixel is parallel to the u-axis, and the other side forms an angle θ with the u-axis. The inclination angle θ is used for inverse transformation to transform the parallelogram pixel into a rectangular pixel:

[0021]

[0022] Among them, the mapping matrix (x r ,yr ) is the inverse mapping matrix R -1 The coordinates of the rectangular pixel after processing, (u, v) is the coordinates of the image point in the pixel coordinate system;

[0023] Let d x =t×d y , using the transformation factor t to convert rectangular pixels into square pixels:

[0024]

[0025] Among them, the mapping matrix (x q ,y q ) is the inverse mapping matrix Q -1 The pixel coordinates after processing.

[0026] The method for calculating the coordinates of the vanishing point in step 3 is as follows: suppose the homogeneous coordinates of the corresponding points X and X' on the world plane and the pixel plane are (x, y, 1) and (x', y', 1), respectively, and calculate the homography matrix H by X = HX'; according to the characteristics of perspective projection, in the world coordinate system, suppose the infinite points in the x-axis and y-axis directions are p x and p y , p x and p y The homogeneous coordinates are [1 0 0] T and [0 1 0] T , the projection of the infinite point on the pixel plane is the vanishing point, and the coordinates of the vanishing point in the horizontal and vertical directions are calculated using the homography matrix H (h vx ,h vy ) and (v vx ,v vy );

[0027] The homography matrix H and the infinite point p in the horizontal direction x Calculate w x =Hp x , and the three-dimensional homogeneous coordinates w x =[w x1 ,w x2 ,w x3 ] T , w x1 、w x2 、w x3 They are the three-dimensional homogeneous coordinates w x The three elements of , so the coordinates of the vanishing point in the horizontal direction are

[0028] The total distortion model in step 4 is: d ,y d) corresponds to the point coordinates (x u ,y u )for:

[0029]

[0030] Among them, (x q ,y q ) is the measurable distortion coordinate on the distortion line l after processing in step 2, δ x and δ y is the distortion displacement on the detected image, and:

[0031]

[0032] Among them, x e =x q -c x ,y e =y q -c y , (c x ,c y ) is the coordinate point of the distortion center, k1, k2 are the radial distortion coefficients, p1, p2 are the eccentric distortion coefficients, s1, s2 are the thin prism distortion coefficients.

[0033] The fitting method of the distortion line in step 3 is: using the calculated vanishing point (h vx ,h vy ) is used to constrain the distortion line, and the equation of the updated fitted distortion line l is:

[0034] Among them, x and y are the independent variable and dependent variable of the fitted line respectively, and a and b are the parameters of the fitted line equation;

[0035] like That is, b = 1, using the square pixel coordinates (x1, y1), (x2, y2), ..., (x n ,y n ), update the parameters of the distortion line l equation:

[0036]

[0037] On the contrary, a=1, calculate the parameters of the equation of the distorted line l:

[0038]

[0039] Among them, n represents the number of coordinate points on each straight line.

[0040] The objective function of the linear residual in step 4 is:

[0041]

[0042] subject to a k h vx +b k h vy +c k =0

[0043] Where m is the number of distorted lines l, n is the number of points on the distorted line l, (h vx ,h vy ) is the coordinate of the vanishing point in the horizontal direction, (x u ,y u ) is the coordinate of the straight line point after correction by the distortion model, (a k ,b k ,c k ) is the result of normalizing the parameters a and b of the distortion line l and calculating the vanishing point together, and:

[0044]

[0045] The optimization method in step 4 is to use the Optimize algorithm in the SciPy function library to minimize the objective function of the overall line residual using the Nelder-Mead method and calculate the total distortion model parameters k1, k s ,p1,p2,s1,s2, tilt angle θ and change factor t; when the distortion parameters converge or the iteration reaches the preset maximum number of iterations, the iterative process ends; otherwise, return to step 2 to continue the iteration until the iterative process ends.

[0046] Compared with the prior art, the present invention has the following beneficial effects: firstly, in order to solve the problem of camera distortion center positioning error, the sub-pixel level point coordinate detection is performed on the detected image, the coordinates of the distortion center are accurately estimated by using the relationship equation between the basic matrix and the distortion center, and the least square method is used to approximate the real distortion center; secondly, in order to solve the problem that the pixels of the camera after digital discretization are not squares but parallelograms under the projective transformation, the inclination angle and transformation factor of the pixel unit on the detected image are calculated to form the inverse mapping matrix of the digital image discrete imaging process, and the non-standard pixel unit is converted into a square pixel unit; finally, on the basis of the full distortion model, the dual characteristics are used to establish a functional constraint relationship with the straight line according to the geometric constraints of the vanishing point, the Optimize algorithm in the SciPy function library is used to minimize the line residual using the Nelder-Mead method, and the distortion parameters are iteratively optimized to achieve accurate distortion correction of the camera. The present invention minimizes the discrete errors in the horizontal and vertical directions of the point coordinates, so that the pixels processed during distortion correction are all strictly square pixels; based on the full distortion model, the dual characteristics of straightness and vanishing points are used to accurately estimate the distortion parameters, and the measurement error caused by lens distortion in the optical system is eliminated. The feasibility of the method of the present invention and the effectiveness of the correction results have been verified through a large number of experiments. The method can quickly and accurately realize lens distortion correction, has good correction effect and strong robustness, and can be applied in solving specific practical problems related to image analysis and processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0048] Figure 1 It is a flow chart of the present invention.

[0049] Figure 2 A grayscale image of the image captured by the camera.

[0050] Figure 3 The figure is a comparison diagram of the transformation factor and the tilt angle before and after correction of the present invention.

[0051] Figure 4 Schematic diagram of the average straightness error of the present invention.

[0052] Figure 5 1 is a curve showing the variation of the cross ratio at different centers of the present invention, wherein (a) is the row cross ratio and (b) is the column cross ratio.

[0053] Figure 6 The figure is a curve showing the change of the cross ratio between the present invention and the document [1], wherein (a) is the row cross ratio and (b) is the column cross ratio.

[0054] Figure 7 This is a comparison diagram of the image point coordinates before and after the distortion correction of the present invention.

[0055] Figure 8 for Figure 2 The grayscale image after correction of the original image in . DETAILED DESCRIPTION

[0056] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0057] like Figure 1 As shown in the figure, a full-model lens distortion correction method using dual-characteristic joint measurement is proposed. The chessboard pattern is used to achieve fast and accurate lens distortion correction. The main idea is: first, in order to solve the problem of camera distortion center positioning error, the basic matrix and the equation of the distortion center are used to determine the accurate distortion center; second, in order to solve the discrete error caused by the digital discretization of the CCD camera, the coordinates are distorted from square pixels to parallelogram pixels, and the inclination angle and transformation factor are calculated to form the inverse mapping matrix of the digital image discrete imaging process to avoid the situation of non-strict square pixels, and the discrete errors in the horizontal and vertical directions are minimized. Minimize; finally, on the basis of the full distortion model, a functional constraint relationship with the straight line is established according to the geometric constraints of the vanishing points, and the Optimize algorithm in the SciPy function library is used to minimize the straight line residual using the Nelder-Mead method, and the distortion parameters, tilt angles and transformation factors are iteratively optimized. When the number of iterations reaches the upper limit or the distortion parameters converge, the iterative process ends, thereby achieving accurate distortion correction of the camera; otherwise, according to the current correction result and the calculated distortion center, continue to iterate its distortion correction steps until the iterative process ends, and accurate distortion coefficients, tilt angles and transformation factors can be obtained. In view of the situation that the traditional method assumes that each pixel point is a parallelogram pixel when processing lens distortion, the present invention studies and processes it, so that the point pixels processed during distortion correction are all strictly square pixels, and the distortion center is accurately estimated based on the full distortion model, with high precision and strong robustness, which can be applied in solving specific practical problems related to image analysis and processing. The specific steps of the present invention are:

[0058] Step 1: Perform sub-pixel point coordinate detection on the chessboard image to be detected, and determine the distortion center using the least squares method based on the relationship equation between the basic matrix and the distortion center.

[0059] A1: Grayscale the image to be detected. The image to be detected is obtained by collecting the chessboard pattern with a CCD camera. The obtained grayscale image is as follows: Figure 2 As shown, the purpose of image processing is to reduce the amount of information stored in color images. In order to better detect point coordinates at the sub-pixel level, the cvtColor function in the opencv function library is used to implement image grayscale operations.

[0060] A2: Perform sub-pixel point coordinate detection on the obtained grayscale image, that is, detect the coordinates of the square vertices in the chessboard image in the detected image to improve the detection accuracy. The detected points are accurate floating-point corner points, which are more accurate than the point coordinates detected by ordinary means. Sub-pixel point coordinate detection is achieved through the cornerSubPix function in the opencv function library.

[0061] A3: The point coordinates obtained in step A2 are formed into a point set E, and the coordinates of the distortion center are calculated using the point set E and the equation of the distortion center. During the projection process, each point moves from its "initial" position to its "final" position, that is, from the undistorted point x u Move to distortion point x d , which is due to the influence of the distortion center c, also known as the "expansion center", or more generally, the pole. The specific calculation steps are as follows:

[0062] Let x w x is the known coordinate point on the chessboard pattern corresponding to the detected image. Points in the horizontal or vertical direction on the chessboard pattern can be used. u is the image point after the distortion model correction, x d is the corresponding distortion point in the pixel coordinate system on the chessboard image. Due to the influence of the distortion center c, the undistorted image point x u Move to distortion point x d , thus:

[0063] x d =c+k(x u -c);

[0064] In the above formula, multiply [x d ] T [c] x , according to the known coordinate point x w and the distortion point x d The corresponding relationship can get the transformation matrix H, and then according to the transformation matrix H and x c Back-projection can obtain an inaccurate undistorted image point xu , that is, x u =Hx w , we can get the equation of the distortion center:

[0065] k[x d ] T [c] x Hx w =[x d ] T [c] x Hx w =0;

[0066] Where k is the radial distortion coefficient, [c] x represents an antisymmetric 3×3 matrix. Let F = [c] x H, F is called the basic matrix, and the distortion center is extracted as the left pole of the basic matrix F. The calculation method of the transformation matrix H is the same as the homography matrix in step 3 C1. The distortion center is extracted as the left pole of the basic matrix F, so we need to find the left null space generator of vector c. In fact, no vector satisfies this equation, so we need to use the least squares method to solve it. Use the least squares method to find the equation c T f j = 0, where j = 1, ..., 3, f j is the jth column of the fundamental matrix F. For each solution, we can obtain T The least squares solution of the three equations for F=0 can determine the distortion center.

[0067] Traditionally, when performing distortion correction on a lens, the image center or principal point coordinates are always taken as the distortion center by default. However, in the actual projection process, the manufacturing tolerance of the lens leads to different imaging characteristics, so that the center is usually in different positions. Whether the distortion center is in the correct position determines whether the experimental results can be optimal. The consequence is that it is impossible to use distortion parameter compensation to correct it, nor is it possible to use external orientation parameters to correct it.

[0068] Step 2: Calculate the tilt angle and transformation factor of the parallelogram pixel of the point coordinates on the detected image, construct the inverse mapping matrix of the pixel point in the discretization process, convert the non-standard pixel unit into a square pixel unit, and minimize the discrete error in the horizontal and vertical directions of the coordinate point.

[0069] The calculation method for converting the coordinates of the distorted points in the point set E into square pixels is:

[0070] B1: For the case where the pixels after digital discretization of the camera are parallelograms instead of squares, minimize the discretization errors in the horizontal and vertical directions of the coordinate points. Let dx =t×d y , an inverse transformation is formed on the pixel plane to calculate the transformation factor t and tilt angle θ of the parallelogram pixel, where d x ,d y are the distances between two adjacent pixels in the x and y directions, respectively. Assuming that the initial values ​​of the transformation factor t and the tilt angle θ are 1 and π / 2, respectively, with each iteration of the present invention, the transformation factor t and the tilt angle θ also change accordingly, thereby processing non-square pixels and minimizing the discrete errors in the horizontal and vertical directions. Figure 3 As shown in FIG. 1 , the comparison diagrams of the tilt angle and the transformation factor calculated by the coordinates of all points before and after the distortion correction using the method of the present invention are shown in FIG. Figure 3 It can be seen that the changes in the tilt angle and transformation factor of the pixel unit before and after processing, the tilt angle is stable near 90 degrees, and the transformation factor fluctuates around 1.00, indicating that the pixel unit after processing gradually approaches the square pixel from the parallelogram pixel.

[0071] B2: Due to the limitation of manufacturing process, the digital discretization process of CCD camera makes the coordinate point pixel of the detected image not a rectangle but a parallelogram. The image point (x, y, 1) T The discretized coordinates are (u,v,1) T . Let the length and width of the parallelogram pixel be d x and d y , one side of the parallelogram is parallel to the u-axis, and the other side forms an angle θ with the u-axis. At this time, the inclination angle θ is used for inverse transformation to transform the parallelogram pixels into rectangular pixels. By using the inverse mapping matrix R of the image point -1 To process parallelogram pixels, as shown in the following formula:

[0072]

[0073] in, (x r ,y r ) is the inverse mapping matrix R -1 The coordinates of the rectangular pixels after processing. At this time, the pixel coordinate point (x r ,y r ) is a rectangular pixel, and its tilt angle is π / 2 (or infinitely close to π / 2).

[0074] B3: According to step B2, at this time (x r ,y r ) is the rectangular pixel coordinate point, then let d x =t×d y, using the transformation factor t to convert rectangular pixels into strictly square pixels, so that the pixels of the coordinate points are standardized into strict squares, more refined pixel coordinates can be obtained. The present invention calculates the inverse mapping matrix Q -1 Converting rectangular pixels to strictly square pixels is shown in the following formula:

[0075]

[0076] in, (x q ,y q ) is the inverse mapping matrix Q -1 After the pixel coordinates are processed in this step, the pixel coordinates are all strictly square pixels.

[0077] By operating the inverse mapping matrix, an inverse operation can be performed on the pixel plane, so that the coordinate point (x d ,y d ), the discrete errors in the horizontal and vertical directions of the pixel unit are minimized, that is, the tilt angle θ of the pixel unit is π / 2 (or infinitely close to π / 2), and the aspect ratio of the pixel unit is 1 (or infinitely close to 1). At this time, the coordinates of the distortion points at each sub-pixel level are strictly square pixels.

[0078] Step 3: Calculate the vanishing points based on the salient features of perspective projection and the homography matrix, and use the vanishing points to update the fitted straight line equation to improve the accuracy of straight line estimation. The calculation steps are as follows:

[0079] C1: The homography matrix is ​​used to describe the relationship between the observed target plane and the image imaging plane. A plane projection transformation can be represented by a 3×3 non-singular matrix. Assuming that the homogeneous coordinates of the corresponding points X and X' on the world plane and the pixel plane are (x, y, 1) and (x', y', 1) respectively, the homography matrix H can be calculated by X=HX'.

[0080] The present invention uses the detected point coordinates and their corresponding world point coordinates to calculate the initial value of the homography matrix H. After entering the iterative correction process, the point coordinates after distortion correction in the previous iteration and their corresponding world point coordinates are used to update the element values ​​of the homography matrix. In this way, the most accurate homography matrix in the distortion correction iterative process can be calculated.

[0081] C2: According to the characteristics of perspective projection (the image of an object extending to infinity may appear within a limited range), the parallel lines in the 3D scene projected onto the image plane intersect at a point, which is called a vanishing point. Two parallel lines in the scene intersect at a vanishing point in the image. In the 3D world coordinate system, let the infinite points in the x-axis and y-axis directions be p x and p y , their homogeneous coordinates are [1 0 0] T and [0 1 0] T , the projection of the infinite point on the image plane is the vanishing point, which is related to the homography matrix H. That is, the coordinates of the vanishing points in the horizontal and vertical directions can be calculated using the homography matrix H (h vx ,h vy ) and (v vx ,v vy ). The selection of vanishing point coordinates is related to the estimation of the straight line equation. Here, the present invention only describes the case of b=1, that is, the straight line and the vanishing point in the horizontal direction are used at this time. Correspondingly, the situation in the vertical direction is similar to this.

[0082] The homography matrix H and the infinite point p in the horizontal direction x , we can calculate w x =Hp x , H is a 3×3 homography matrix, that is, w x =[w x1 ,w x2 ,w x3 ] T , w x is a three-dimensional homogeneous coordinate, according to the formula w x =Hp x To calculate, w x1 、w x2 、w x3 They are w x The three elements of , so the coordinates of the vanishing point in the horizontal direction are Using the calculated vanishing point (h vx ,h vy ) is used to constrain the distorted straight line. The equation of the updated fitted distorted straight line is calculated by the following formula:

[0083]

[0084] Among them, x and y are the independent variable and dependent variable of the fitted line, respectively, a and b are the parameters of the fitted line equation, and the calculation formula of the equation parameter a of the fitted line l is described in step D2 of step 4 by updating and optimizing the line l. The fitted line is an optimization of the distorted line, which can better meet the requirements of the algorithm and pave the way for defining the overall line residual objective function.

[0085] Step 4: Based on the full distortion model, the geometric constraints of the vanishing points are used to establish a functional constraint relationship with the distortion line. The overall line residual is minimized using the optimization method to obtain the distortion parameters, thus achieving accurate distortion correction of the camera. The specific implementation steps are as follows:

[0086] D1: Under the full distortion model, calculate the radial, eccentric and thin prism distortion coefficients:

[0087] Assume (x d ,y d ) is the measurable coordinate of the distortion point on the straight line l. The distortion point is on the distortion line l. Due to the existence of lens distortion, l is changed from a straight line to a curved curve, so the ultimate goal is to change it back to a straight line. It is necessary to use step C2 to optimize the straight line l. As the iterative correction proceeds, the equation of the straight line l is gradually optimized and turned back to a straight line. (x u ,y u ) is the same as the distortion point (x q ,y q ) is the point coordinate after correction by the distortion model corresponding to , then according to the full distortion model, the corrected straight line point coordinate (x u ,y u ) can be defined as:

[0088]

[0089] Among them, δ x and δ y Represents the distortion displacement on the plane image (the detected image) and is expressed by the following formula:

[0090]

[0091] Among them, x e =x q -c x ,y e =y q -c y , (c x ,c y ) is the distortion center obtained in step 1, k1, k2 are radial distortion coefficients, p1, p2 are eccentric distortion coefficients, s1, s2 are thin prism distortion coefficients.

[0092] The advantage of this step is that all distortion coefficients are solved based on the full distortion model, which can meet the actual needs of camera measurement and achieve better lens distortion correction.

[0093] D2: Before calculating the distortion coefficient, it is necessary to use the vanishing point to update the straight line equation. When optimizing the straight line parameters, there are two cases: a=1 or b=1. Here we only describe the case of b=1, that is, the straight line and the vanishing point in the horizontal direction. Correspondingly, the situation in the vertical direction is similar.

[0094] In order to make the optimized fitting straight line as close to the ideal straight line as possible, if That is, b = 1. From step D1, using the square pixel coordinates (x1, y1), (x2, y2), ..., (x n ,y n ), and update the parameter a of the fitted straight line equation as follows:

[0095]

[0096] On the contrary, a=1, the parameter b of the straight line equation is calculated by the following formula:

[0097]

[0098] Among them, n represents the number of coordinate points on each straight line.

[0099] D3: Minimize the overall line residual objective function to obtain accurate lens distortion coefficients, tilt angles and transformation factors.

[0100] Since the measurement points inevitably carry measurement errors, only one straight line l = (a, b, c) can be estimated to make the set of points X (in order to derive the defined overall line residual objective function, X refers to the set of points formed after processing in step B3) as "close" to the straight line as possible. The degree of "closeness" is measured by algebraic distance. The algebraic method of straight line estimation is to find a straight line that minimizes the sum of the squares of the algebraic distances from the measurement points to it, that is, to obtain the objective function of the straight line residual defined by the following formula:

[0101]

[0102] subject to a k h vx +b k h vy +c k =0

[0103] Where m is the number of distorted lines l, n is the number of points on the line l, (h vx ,h vy ) is the calculated horizontal coordinate of the vanishing point, (x u ,y u ) is the coordinate of the straight line point after the distortion model correction in step D1, (a k ,bk ,c k ) is the result of normalizing the straight line equation parameters a and b and calculating the vanishing point, as shown in the following formula:

[0104]

[0105] Next, the Optimize algorithm in the SciPy library is used to minimize the line residual using the Nelder-Mead method to calculate the total distortion model parameters k1,k s ,p1,p2,s1,s2, tilt angle θ and change factor t, through multiple repeated iterations, the most accurate distortion coefficient is calculated to achieve accurate distortion correction of the camera. When the distortion parameters converge or the iteration reaches the preset maximum number of iterations (set to 5 in the present invention), the iterative process ends; otherwise, according to the distortion correction process, continue to iterate its distortion correction steps until the iterative process ends. The optimal distortion coefficient is calculated through multiple repeated iterations, thereby achieving accurate distortion correction of the camera. In step D1, the expression of the full distortion model is a function of the distortion coefficients k1, k2, p1, p2, s1, s2 and the transformation factor t and the tilt angle θ, so the values ​​of these parameters can be obtained after each correction and used as the initial value of the next iteration. Return to step 2 to update the inverse mapping matrix R -1 and the inverse mapping matrix Q -1 , thereby converting the distorted points into square pixels.

[0106] The advantage of this step is that the function of the distortion coefficient in step D1 is optimized and solved to obtain the most accurate distortion coefficient of the lens, so that the distortion correction effect is better and more robust.

[0107] According to the obtained distortion correction results, the cross ratio value and the average straightness error of the image point coordinates after each correction are compared to evaluate the distortion correction results. The evaluation factors of the distortion correction results of the present invention are composed of the cross ratio value and the average straightness error, and the concept definition is as follows:

[0108] a: The cross ratio is also called the harmonic ratio. It is an invariant of fractional linear transformation and one of the basic invariants in projective geometry. Any four co-pointed lines have a certain cross ratio. i The definition of the cross ratio is shown in the following formula.

[0109]

[0110] In the formula, Under any projective transformation of the line, the value of the cross ratio remains unchanged: if x' = H 2×2x, then Cross(x1',x'2,x'3,x'4) = Cross(x1,x2,x3,x4), which is a projective coordinate system in which the cross ratio does not depend on the straight line.

[0111] b: The present invention also uses the average straightness error to more objectively evaluate the effect of distortion correction. The average straightness error refers to the average distance of the corrected straight line from the sub-pixel edge point to the corresponding regression line. According to the definition, the straightness error is defined as the root mean square distance from a set of distorted edge points corresponding to the same physical line to its overall linear regression, such as Figure 4 shown.

[0112] Given a strictly square pixel point (x1,y1),...,(x N ,y N ), and thus the regression line is obtained:

[0113] αx l +βy l -γ=0

[0114] And α=sinθ, β=cosθ, γ=A x sinθ+A y cosθ, where

[0115]

[0116]

[0117] Since (α, β) is a unit vector, from point (x i ,y i ) to the straight line is

[0118] S li =αx li +βy li -γ

[0119] Given L straight lines, N l is the number of points on each line l, then the sum of the squares of the symbol matrix is

[0120]

[0121] Therefore, the average straightness error (ASE) is defined as

[0122]

[0123] in, is the total number of points.

[0124] According to the method of the present invention, a zoom lens is used to shoot a chessboard plane in different postures and positions, and the sub-pixel level point coordinate detection is performed on it, because it has strong robustness to noise. According to the method proposed by the present invention, based on the full distortion model, the correction results of the model at different centers are compared. The comparison experimental results of the cross ratio are as follows: Figure 5 As shown, the comparison results of the average straightness error at different centers are shown in Table 1.

[0125] Table 1 Average straightness error ASE (comparison of different centers)

[0126]

[0127] Depend on Figure 5 As can be seen from Table 1, the cross ratio value of the result after the distortion center correction is closer to the ideal cross ratio value, which is more accurate and precise than the correction result using the image center and the principal point coordinates, and the effect is better. This can be observed by comparing the average straightness error. At the same time, the correction result of the present invention is also compared with Liu's method---Liu Lianxiong, Hu Changhua, He Chuan, et al. An improved non-measurement distortion correction method based on straight line features [J]. Acta Optica Sinica, 2017, 37(09): 0915001.], and the comparison results of the cross ratio value are as follows. Figure 6 The average straightness error is shown in Table 2.

[0128] Table 2 Average straightness error ASE (comparison with Liu's method)

[0129]

[0130] Depend on Figure 6 From Table 2, we can see that after iteration, the row and column cross ratios of the present invention differ from the ideal cross ratio by only 5.163775e-05 at the minimum and 8.542710e-04 at the maximum, with the minimum and maximum standard deviations of 5.701253e-05 and 4.855813e-04, respectively. The row and column cross ratios of Liu's method differ from the ideal cross ratio by 7.745470e-04 at the minimum and 2.609315e-02 at the maximum, with the minimum and maximum standard deviations of 6.449447e-04 and 2.665085e-00, respectively. Figure 6 From Table 2, we can see that the correction accuracy of the present invention is more accurate than that of Liu's method. For the cross ratio value, it is probably more optimized by one to two precisions; according to the average straightness error, it is more optimized by 83.878% than Liu's method.

[0131] According to the point coordinates detected at the initial time (the initial point coordinates at the sub-pixel level obtained in step A2), six horizontal approximate straight lines and nine vertical approximate straight lines are outlined, such as Figure 7 According to the full distortion model, the result of each distortion correction is observed by remapping its pixel coordinates, which is generated by detecting a single feature point. Figure 7 The dashed lines (--) are the results after correction by the present invention, and they are all from Figure 2 The actual image captured in the image, the result of distortion correction is as follows Figure 8 As shown, it can be seen that the method proposed in the present invention can achieve accurate correction. After a large number of experimental verifications, the cross ratio value of the result after correction using the distortion center is closer to the ideal cross ratio value, which is better than the correction results using the image center and the principal point coordinates. The correction effect is good, efficient and robust, and the distortion correction effect is significantly better than other algorithms.

[0132] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A full-model lens distortion correction method using dual-characteristic joint measurement, characterized in that: The steps are as follows: Step 1: Perform sub-pixel point coordinate detection on the detected image, and determine the distortion center using the least squares method based on the relationship equation between the basic matrix and the distortion center; Step 2: Calculate the tilt angle and transformation factor of the parallelogram pixel of the point coordinates on the detected image, form the inverse mapping matrix of the discrete imaging process of the digital image, and transform the non-standard pixel unit into a square pixel unit; Step 3: Calculate the coordinates of the vanishing points based on the salient features of the perspective projection and the homography matrix, and use the vanishing points to update the fitted distortion line; Step 4: Based on the full distortion model, a functional constraint relationship with the distortion line is established according to the geometric constraints of the vanishing point. The optimization method is used to minimize the straight line residual to obtain the distortion parameters, tilt angle and transformation factor, so as to achieve accurate distortion correction of the camera.

2. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 1, characterized in that: The detected image is obtained by collecting the chessboard pattern through a CCD camera.

3. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 2, characterized in that: In the step 1, the detected image is grayed to obtain a grayscale image, and the grayscale image is subjected to sub-pixel point coordinate detection, i.e., the coordinates of the square vertices in the chessboard image are detected to form a point set E, and the coordinates of the distortion center are calculated using the point set E and the equation of the distortion center.

4. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 3, characterized in that: The method for calculating the coordinates of the distortion center is: Let x w is the known coordinate point on the chessboard pattern corresponding to the detected image, x u is the image point after the distortion model correction, x d is the corresponding distortion point in the pixel coordinate system on the chessboard image; due to the influence of the distortion center c, the undistorted image point x u Move to distortion point x d For: x d =c+k(x u -c); left multiply [x d ] T [c] x , we get the equation for the distortion center: k[x d ] T [c] x Hx w =[x d ] T [c] x Hx w =0; Where k is the radial distortion coefficient, [c] x represents an antisymmetric 3×3 matrix, the basic matrix F = [c] x H; H is the transformation matrix and x u =Hx w , the distortion center is extracted as the left pole of the basic matrix F; Using the least squares method, we get equation c T f j = 0, j = 1, ..., 3, f j is the jth column of the basic matrix F; the least squares solution of all equations is the coordinate point of the distortion center c.

5. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 3, characterized in that: The implementation method of converting the non-standard pixel unit into a square pixel unit in step 2 is: Assume that the length and width of the parallelogram pixel are d x and d y , one side of the parallelogram pixel is parallel to the u-axis, and the other side forms an angle θ with the u-axis. The inclination angle θ is used for inverse transformation to transform the parallelogram pixel into a rectangular pixel: Among them, the mapping matrix (x r ,y r ) is the inverse mapping matrix R -1 The coordinates of the rectangular pixel after processing, (u, v) is the coordinates of the image point in the pixel coordinate system; Let d x =t×d y , using the transformation factor t to convert rectangular pixels into square pixels: Among them, the mapping matrix (x q ,y q ) is the inverse mapping matrix Q -1 The pixel coordinates after processing.

6. The full model lens distortion correction method using dual-characteristic joint measurement according to any one of claims 1, 4 or 5, characterized in that: The method for calculating the coordinates of the vanishing point in step 3 is as follows: suppose the homogeneous coordinates of the corresponding points X and X' on the world plane and the pixel plane are (x, y, 1) and (x', y', 1), respectively, and calculate the homography matrix H by X = HX'; according to the characteristics of perspective projection, in the world coordinate system, suppose the infinite points in the x-axis and y-axis directions are p x and p y , p x and p y The homogeneous coordinates of are [100] T and [010] T , the projection of the infinite point on the pixel plane is the vanishing point, and the coordinates of the vanishing point in the horizontal and vertical directions are calculated using the homography matrix H (h vx ,h vy ) and (v vx ,v vy ); The homography matrix H and the infinite point p in the horizontal direction x Calculate w x =Hp x , and the three-dimensional homogeneous coordinates w x =[w x1 ,w x2 ,w x3 ] T , w x1 、w x2 、w x3 They are the three-dimensional homogeneous coordinates w x The three elements of , so the coordinates of the vanishing point in the horizontal direction are 7. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 1, characterized in that: The total distortion model in step 4 is: d ,y d ) corresponds to the point coordinates (x u ,y u )for: Among them, (x q ,y q ) is the measurable distortion coordinate on the distortion line l after processing in step 2, δ x and δ y is the distortion displacement on the detected image, and: Among them, x e =x q -c x ,y e =y q -c y , (c x ,c y ) is the coordinate point of the distortion center, k1, k2 are the radial distortion coefficients, p1, p2 are the eccentric distortion coefficients, s1, s2 are the thin prism distortion coefficients.

8. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 7, characterized in that: The fitting method of the distortion line in step 3 is: using the calculated vanishing point (h vx ,h vy ) is used to constrain the distortion line, and the equation of the updated fitted distortion line l is: Among them, x and y are the independent variable and dependent variable of the fitted line respectively, and a and b are the parameters of the fitted line equation; like That is, b = 1, using the square pixel coordinates (x1, y1), (x2, y2), ..., (x n ,y n ), update the parameters of the distortion line l equation: On the contrary, a=1, calculate the parameters of the equation of the distorted line l: Among them, n represents the number of coordinate points on each straight line.

9. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 8, characterized in that: The objective function of the linear residual in step 4 is: subject to a k h vx +b k h vy +c k =0 Where m is the number of distorted lines l, n is the number of points on the distorted line l, (h vx ,h vy ) is the coordinate of the vanishing point in the horizontal direction, (x u ,y u ) is the coordinate of the straight line point after correction by the distortion model, (a k ,b k ,c k ) is the result of normalizing the parameters a and b of the distortion line l and calculating the vanishing point together, and:

10. The full model lens distortion correction method using dual-characteristic joint measurement according to claim 9, characterized in that: The optimization method in step 4 is to use the Optimize algorithm in the SciPy function library to minimize the objective function of the overall line residual using the Nelder-Mead method and calculate the total distortion model parameters k1, k s ,p1,p2,s1,s2, tilt angle θ and change factor t; when the distortion parameters converge or the iteration reaches the preset maximum number of iterations, the iterative process ends; otherwise, return to step 2 to continue the iteration until the iterative process ends.

Images