A calibration method for non-coaxial cameras

Abstract

A calibration method for a non-coaxial camera includes: obtaining a calibration plate image taken by the non-coaxial camera; obtaining feature points in the calibration plate image and their image coordinates and world coordinates; calculating a homography matrix; decomposing the homography matrix according to a preset transformation model from world coordinates to image coordinates to obtain internal and external parameters of the non-coaxial camera, wherein the internal parameters include a tilt matrix representing the transformation from a tilted image plane coordinate system to a non-tilted image plane coordinate system, the tilted image plane is an image plane perpendicular to the optical axis of the lens, and the non-tilted image plane is an image plane of the non-coaxial camera; performing nonlinear optimization on the distortion coefficient and the decomposed internal and external parameters to obtain the final internal parameters, external parameters and distortion coefficient. Since the tilted image plane and the non-tilted image plane are introduced, the tilt matrix is ​​added to describe the transformation from the tilted image plane coordinate system to the non-tilted image plane coordinate system, which can effectively improve the calibration accuracy of the non-coaxial camera.

Description

[0001] This application is a divisional application of the original application. The application number of the original application is: 202111526560.7. The name of the invention is: A calibration method for a non-coaxial camera. The application date is: December 15, 2021. Technical Field

[0002] The present invention relates to the technical field of camera calibration, and in particular to a calibration method for a non-coaxial camera. Background Art

[0003] In the image measurement process and machine vision applications, in order to determine the three-dimensional geometric position of a point on the surface of a spatial object, it is necessary to establish a geometric model of camera imaging, that is, to determine the correspondence between the three-dimensional geometric position of a point on the surface of a spatial object and its corresponding point on the image. In this way, after obtaining the image coordinates of the image taken by the camera, the corresponding three-dimensional spatial coordinates can be inferred according to the geometric model of camera imaging. The parameters in the geometric model are the parameters of the camera, and the process of determining the parameters of the camera is called camera calibration. The calibration of camera parameters is a very critical link. The accuracy of the calibration results and the stability of the calibration algorithm directly affect the accuracy of the results produced by the camera. Therefore, doing a good job of camera calibration is a prerequisite for doing a good job of subsequent work. Camera calibration is often carried out using calibration plates. Calibration plates are widely used in machine vision, image measurement, photogrammetry, three-dimensional reconstruction, etc. The camera captures an image of a calibration plate with a fixed-pitch pattern array. After calculation by the calibration algorithm, the geometric model of camera imaging can be obtained, thereby obtaining high-precision measurement and reconstruction results. At present, calibration plates with checkerboard patterns or solid circle array patterns are usually used for camera calibration. The checkerboard calibration plate obtains feature points by locating the checkerboard corners, and the circle array calibration plate obtains feature points by locating the centers of the circles. After determining the coordinates of the feature points and their correspondence with the world coordinates, subsequent calibration work can be carried out. Summary of the invention

[0004] The present application provides a non-coaxial camera calibration method, which can be used to calibrate a non-coaxial camera.

[0005] According to a first aspect, an embodiment provides a calibration method for a non-coaxial camera, wherein the non-coaxial camera includes an image plane and a lens, wherein a normal vector of the image plane and an optical axis of the lens are not coaxial, and the lens is an image-space telecentric lens or a bilateral telecentric lens, and the calibration method includes:

[0006] Obtain a calibration plate image taken by a non-coaxial camera;

[0007] Acquire feature points in the calibration plate image, as well as image coordinates and corresponding world coordinates of the feature points;

[0008] Calculate the homography matrix H according to the image coordinates of the feature points and the corresponding world coordinates;

[0009] According to the preset transformation model from world coordinates to image coordinates, the homography matrix is ​​decomposed and calculated to obtain the intrinsic parameters and extrinsic parameters of the non-coaxial camera. The transformation model is:

[0010]

[0011] Among them, the homography matrix (r,c) T is the image coordinate of the feature point, (x w ,y w ,z w ) T is the world coordinate of the feature point, is the transformation matrix from the world coordinate system to the camera coordinate system, R is the rotation matrix, t is the displacement matrix, z c is the z coordinate of the feature point in the camera coordinate system, is the transformation matrix from the camera coordinate system to the tilted image plane coordinate system, f is the focal length of the non-coaxial camera, H tilt is a tilt matrix, representing the transformation from the tilted image plane coordinate system to the non-tilted image plane coordinate system, wherein the tilted image plane is the image plane perpendicular to the optical axis of the lens, and the non-tilted image plane is the image plane of the non-coaxial camera, is the transformation matrix from the non-tilted image plane coordinate system to the image coordinate system, s x and y are the pixel sizes of the non-coaxial camera in the horizontal and vertical directions, respectively. x ,c y ) is the principal optical axis point, For the internal reference part, is the external parameter part, the tilt matrix H tilt Specifically

[0012]

[0013] Among them, q 11 ,q 12 ,q 21 ,q 22 is an element in the rotation matrix Q, which represents the rotation transformation of the tilted image plane relative to the original coordinate system, and

[0014]

[0015] Wherein ρ represents the angle of rotation around the Z axis, τ represents the angle of rotation around the X axis, the X axis of the original coordinate system is the horizontal direction of the non-tilted image plane, the Y axis is the vertical direction of the non-tilted image plane, and the Z axis is the vertical line of the non-tilted image plane;

[0016] The distortion coefficient of the non-coaxial camera and the decomposed internal and external parameters are nonlinearly optimized to obtain the final internal and external parameters and distortion coefficient of the non-coaxial camera.

[0017] In one embodiment, decomposing and calculating the homography matrix to obtain the intrinsic parameters and extrinsic parameters of the non-coaxial camera includes:

[0018] Calculate the parameter matrix A according to the following constraints

[0019]

[0020] in

[0021] H=[h1 h2 h3]=A[r1 r2 t],

[0022]

[0023] [r1 r2 t]=[R|t];

[0024] Where h1 is the first column vector of the homography matrix H, h2 is the second column vector of the homography matrix H, h3 is the third column vector of the homography matrix H, r1 is the first column vector of the rotation matrix R, and r2 is the second column vector of the rotation matrix R;

[0025] According to r1=A -1 h1, r2 = A -1 h2,t=A -1 h3 calculates the matrix [r1 r2 t], according to Calculate the tilt matrix H tilt .

[0026] In one embodiment, the external parameters of the non-coaxial camera further include an equivalent rotation axis k and an equivalent axis angle θ, and the internal parameters and external parameters of the non-coaxial camera obtained by decomposing and calculating the homography matrix include:

[0027] Calculate the parameter matrix A according to the following constraints

[0028]

[0029] in,

[0030] H=[h1 h2 h3]=A[r1 r2 t],

[0031]

[0032] [r1 r2 t] = [R|t];

[0033] According to r1=A -1 h1, r2 = A -1 h2,t=A -1 h3 calculates the matrix [r1 r2 t];

[0034] Where h1 is the first column vector of the homography matrix H, h2 is the second column vector of the homography matrix H, h3 is the third column vector of the homography matrix H, r1 is the first column vector of the rotation matrix R, and r2 is the second column vector of the rotation matrix R;

[0035] According to the calculated rotation matrix R, the equivalent rotation axis k and the equivalent axis angle θ are obtained, where the transformation relationship between the rotation matrix R and the equivalent rotation axis k and the equivalent axis angle θ is as follows:

[0036]

[0037] k x , k y , k z are the three components of the equivalent rotation axis k;

[0038] according to Calculate the tilt matrix H tilt .

[0039] In one embodiment, the nonlinear optimization of the distortion coefficient of the non-coaxial camera and the decomposed internal parameters and external parameters to obtain the final internal parameters, external parameters and distortion coefficient of the non-coaxial camera includes:

[0040] The initial value of the distortion coefficient is set in advance, and the decomposed internal and external parameters are used as the initial values ​​of the internal and external parameters. The optimal solution is iteratively solved according to the following loss function to obtain the final internal and external parameters and distortion coefficient of the non-coaxial camera:

[0041]

[0042] Among them, n m is the number of feature points in the calibration plate image, n c is the number of cameras, n0 is the number of calibration plate images taken by the camera, and p j is the coordinate of the feature point in the world coordinate system, e l (l=1,…,n0) represents the external parameters of the calibration plate image in the reference camera, r k (k=1,…,n c ) represents the transformation of the kth camera relative to the reference camera, i k represents the transformation under the kth camera, pjkl is the image coordinate of the jth feature point in the lth calibration plate image taken by the kth camera, v jkl The value is 0 or 1. When the jth feature point is visible in the lth calibration plate image taken by the kth camera, it is 1, otherwise it is 0; function represents the transformation from the image coordinate system to the tilted image plane coordinate system, which includes using the intrinsic parameters to transform the image coordinate system to the tilted image plane coordinate system and using the distortion coefficient to dedistort the tilted image plane coordinate system; the function φ u (p j ,e l ,r k ,i k ) represents the transformation from the world coordinate system to the tilted image plane coordinate system, which includes the use of extrinsic parameters to transform the world coordinate system into the camera coordinate system.

[0043] In one embodiment, the calibration method of the non-coaxial camera further includes: before each iteration, using the calculated distortion coefficient to perform distortion correction on the tilted image plane coordinates of the feature point.

[0044] In one embodiment, according to the formula q k+1 =q k +δ iterates to find the optimal solution, where q k represents the vector composed of the intrinsic parameters, extrinsic parameters and distortion coefficients of the non-coaxial camera at the kth iteration. δ is determined by the formula Jδ = ε, where ε is the number of feature points corresponding to each calibration plate image captured by each camera at the current iteration. With φ u (p j ,e l ,r k ,i k ), J is the Jacobian matrix, which is composed of the Jacobian matrices of each camera, and the Jacobian matrix of the i-th camera is

[0045]

[0046] Among them J i,j and J i ' ,j They respectively represent the partial derivatives of the intrinsic and extrinsic parameters corresponding to the j-th calibration plate image taken by the ith camera.

[0047] In one embodiment, the calibration plate image is a circular array calibration plate image; the feature points in the calibration plate image and the image coordinates of the feature points are obtained by:

[0048] Performing image processing on the calibration plate image to obtain circular feature points therein;

[0049] Perform edge extraction on the circular feature points to obtain edge points of the circular feature points, and use the edge points to perform ellipse fitting to obtain image coordinates of the circular feature points, wherein the image coordinates of the circular feature points refer to the image coordinates of the center of the circular feature points;

[0050] Determine the correspondence between the image coordinates and the world coordinates of the circular feature points;

[0051] The image coordinates of the circular feature points are corrected using the ellipse equation to obtain the final image coordinates of the circular feature points.

[0052] In one embodiment, the error correction of the image coordinates of the circular feature points using the ellipse equation includes:

[0053] Calculate the non-distorted ellipse equation matrix D to the distorted ellipse equation matrix The transformation matrix H D , where the undistorted elliptic equation matrix D satisfies

[0054] p T Dp=0,

[0055] Distorted ellipse equation matrix satisfy

[0056]

[0057] Transformation matrix H D satisfy

[0058]

[0059] Where p is the image coordinate of the circular feature point after error correction, are the image coordinates of the circular feature points before error correction, and

[0060]

[0061] in Λ0 and U satisfy

[0062] U T DU=Λ,

[0063] Λ=diag(λ1,λ2,λ3),

[0064]

[0065] Where λ1, λ2 and λ3 are the eigenvalues ​​of the elliptic equation matrix D, and U is the matrix consisting of the corresponding eigenvectors. and is the elliptic equation matrix The characteristic value of is the matrix consisting of the corresponding eigenvectors;

[0066] According to the following objective function, the image coordinates p of the circular feature point after error correction are obtained. i :

[0067]

[0068] The subscript i represents the i-th point.

[0069] In one embodiment, the error correction of the image coordinates of the circular feature points using the ellipse equation includes:

[0070] According to the following objective function, the edge points of the circular feature points are used to fit the ellipse and obtain the coefficients in the ellipse equation:

[0071]

[0072] st4ac-b 2 =1

[0073] Where a, b, c, d, e, and f are the coefficients of the ellipse equation, (x i ,y i ) is the image coordinate of the edge point, w i is the weight, n is the number of edge points;

[0074] The image coordinates of the center point of the ellipse (r) are calculated according to the following formula d ,c d ):

[0075]

[0076] The image coordinates (r′) of the center point of the unfitted ellipse are calculated according to the following formula: d ,c′ d ):

[0077]

[0078] Where F is the set of points in the area where the circular feature point is located, p i is a point in the set, I(p i ) is point p i The gray value, (x i ,y i ) is point p i The image coordinates of , where the subscript i represents the i-th point;

[0079] According to the image coordinates (r d ,c d ) and (r′ d ,c′d ) and perform error correction on the image coordinates of the circular feature points.

[0080] According to the second aspect, an embodiment provides a computer-readable storage medium, on which a program is stored. The program can be executed by a processor to implement the calibration method of the non-coaxial camera described in the first aspect.

[0081] According to the calibration method and computer-readable storage medium of the non-coaxial camera of the above embodiment, for the non-coaxial camera, when establishing the transformation model from the world coordinate to the image coordinate, the situation that the optical axis of the lens of the non-coaxial camera is not coaxial with the optical axis of the imaging plane is taken into consideration, and the concepts of tilted image plane and non-tilted image plane are introduced, wherein the tilted image plane is an image plane perpendicular to the optical axis of the lens of the non-coaxial camera, and the non-tilted image plane is the imaging plane of the non-coaxial camera, and the tilt matrix H is added. tilt This parameter is used to describe the transformation from the tilted image plane coordinate system to the non-tilted image plane coordinate system. The tilt matrix H tilt is regarded as part of the camera internal parameters; at the same time, the rotation angle of the tilted image plane around the coordinate axis of the original coordinate system is introduced to represent the tilt matrix H tilt , where the X-axis of the original coordinate system is the horizontal direction of the non-tilted image plane, the Y-axis is the vertical direction of the non-tilted image plane, and the Z-axis is the vertical line of the non-tilted image plane. The mathematical model in the existing method is corrected, which can well describe the process of coordinate system conversion in the non-coaxial camera, improve the calibration accuracy of the non-coaxial camera, and reduce the error of the working result generated in the process of applying the non-coaxial camera. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 Schematic diagram of the optical structure of the coaxial camera;

[0083] Figure 2 Schematic diagram of the optical structure of a non-coaxial camera;

[0084] Figure 3 Schematic diagram of the transformation of each coordinate system in the pinhole camera model;

[0085] Figure 4 is a flow chart of a calibration method for a non-coaxial camera in an embodiment;

[0086] Figure 5 is a schematic diagram of tilt transformation;

[0087] Figure 6 A flowchart of a method for extracting high-precision coordinates of feature points of a circular array calibration plate according to an embodiment;

[0088] Figure 7A flowchart of performing image processing on a calibration plate image to obtain circular feature points therein in an embodiment;

[0089] Figure 8 A flowchart of performing image processing on a calibration plate image to obtain circular feature points therein in another embodiment;

[0090] Fig. 9 is a schematic diagram of a circular array calibration plate with triangular markers;

[0091] Fig.10 is a schematic diagram of a circular array calibration plate with hollow points;

[0092] Fig.11 A flowchart for determining the correspondence between the image coordinates and the world coordinates of circular feature points in a circular array calibration plate with triangular markers;

[0093] Fig.12 The present invention is a flowchart for determining the correspondence between the image coordinates and the world coordinates of circular feature points in a circular array calibration plate with hollow points. DETAILED DESCRIPTION

[0094] The present invention is further described in detail below by specific embodiments in conjunction with the accompanying drawings. Wherein similar elements in different embodiments adopt associated similar element numbers. In the following embodiments, many detailed descriptions are for making the present application better understood. However, those skilled in the art can easily recognize that some features can be omitted in different situations, or can be replaced by other elements, materials, methods. In some cases, some operations related to the present application are not shown or described in the specification, this is to avoid the core part of the present application being overwhelmed by too much description, and for those skilled in the art, it is not necessary to describe these related operations in detail, and they can fully understand the related operations according to the description in the specification and the general technical knowledge in the art.

[0095] In addition, the features, operations or characteristics described in the specification can be combined in any appropriate manner to form various implementations. At the same time, the steps or actions in the method description can also be interchanged or adjusted in a manner that is obvious to those skilled in the art. Therefore, the various sequences in the specification and the drawings are only for the purpose of clearly describing a certain embodiment and are not meant to be a required sequence, unless otherwise specified that a certain sequence must be followed.

[0096] The serial numbers assigned to the components in this document, such as "first", "second", etc., are only used to distinguish the objects described and do not have any order or technical meaning. The "connection" and "coupling" mentioned in this application, unless otherwise specified, include direct and indirect connections (couplings). The "image plane" and "imaging plane" mentioned in this document are the same concept.

[0097] The optical axis of the lens in currently used cameras is coaxial with the optical axis of the imaging plane (that is, the normal vector of the imaging plane). However, in the field of machine vision, due to manufacturing or design reasons, the lens of some cameras is not necessarily parallel to the imaging plane, that is, the optical axis of the camera lens is not coaxial with the optical axis of the imaging plane. If this is not taken into account and calibration is performed using traditional methods, errors that cannot be ignored will be introduced.

[0098] The main methods of current camera calibration are designed and calculated based on Zhang Zhengyou's calibration method, which mainly includes the following calculation steps:

[0099] (1) Obtain the homography matrix based on the correspondence between the world coordinates and image coordinates of the feature points in the calibration plate;

[0100] (2) Decomposing the homography matrix and calculating the initial parameters of the internal or external parameters;

[0101] (3) The LM (Levenberg-Marquardt) algorithm is used to perform nonlinear optimization on the initial parameters, and the internal parameters, external parameters and distortion coefficients are iteratively calculated to obtain the final calibration results.

[0102] However, Zhang Zhengyou's calibration method mainly considers the situation that the optical axis of the lens and the optical axis of the imaging plane are coaxial, and does not provide a processing method and mathematical model design for the non-coaxial situation. Therefore, it is very necessary to provide a calibration method that can be used for non-coaxial cameras.

[0103] In the case of no eccentricity or coaxiality, the optical structure of the camera is as follows Figure 1 As shown, the optical axis of the image plane (i.e., the vertical axis of the image plane or the vector) and the optical axis of the lens coincide with each other. In the case of off-axis or non-coaxial, the optical structure of the camera is as follows: Figure 2 As shown, the optical axis of the image plane and the optical axis of the lens do not coincide with each other, but there is a certain angle θ.

[0104] In the case of coaxial, the projection transformation relationship of each coordinate system in the imaging process of the camera can be expressed by the pinhole camera model, such as Figure 3 Point P in the World Coordinate System (WCS) wProject the lens projection center to point P on the imaging plane to obtain point P w The image coordinates q projected onto the imaging plane i , it needs to be converted to the camera coordinate system (CCS) first. The x-axis and y-axis of the camera coordinate system are parallel to the c-axis and r-axis of the image respectively, the z-axis is perpendicular to the imaging plane where the image is located, and the direction of the z-axis is set so that the z coordinates of all points in front of the camera are positive numbers, where the c-axis direction of the image is the horizontal direction of the image, and the r-axis direction is the vertical direction of the image. Figure 3 Medium c Axis, y c Axis and z c The axes represent the x-axis, y-axis, and z-axis of the camera coordinate system. The transformation from the world coordinate system to the camera coordinate system can be expressed as c = c H w p w To represent, where p c =(x c ,y c ,z c ) T is the coordinate in the camera coordinate system, p w =(x w ,y w ,z w ) T is the coordinate in the world coordinate system, c H w It can be represented by the rotation matrix R and the translation matrix t.

[0105] After converting the world coordinate system to the camera coordinate system, it is necessary to convert it to the image plane coordinate system. This is a process of converting 3D coordinates to 2D coordinates. For non-telecentric lenses such as CCTV (Closed Circuit Television) lenses, this transformation can be expressed as:

[0106]

[0107] Where f represents the focal length of the camera lens, (u,v) T Represents the coordinates in the image plane coordinate system.

[0108] For a telecentric lens, this transformation can be expressed as:

[0109]

[0110] Where m represents the magnification of the lens.

[0111] After projection onto the imaging plane, the lens distortion will result in coordinates q c =(u,v) T changes, so that the coordinates formed on the imaging plane are distorted This change can be modeled on the imaging plane alone, which means that no three-dimensional information is needed. For most lenses, their distortion can be adequately approximated as radial distortion. There are usually two models that can be used to describe distortion, one is the division model and the other is the polynomial model. The division model is as follows:

[0112]

[0113] The parameter κ represents the magnitude of radial distortion. If κ is negative, it becomes barrel distortion. If κ is positive, it becomes pincushion distortion. The distortion can be corrected by the following formula:

[0114]

[0115] The polynomial model is as follows:

[0116]

[0117] in k1, k2, k3, p1, p2 are model coefficients. According to the above model, u and v can be solved by Newton's method, and the initial value of the iteration is the undistorted initial value itself.

[0118] Finally, the image plane coordinate system is converted to the image coordinate system (ICS), which can be expressed as follows:

[0119]

[0120] where s x and y are the pixel sizes of the camera in the horizontal and vertical directions, respectively, (c x ,c y ) is the principal optical axis point, generally the center of the image.

[0121] Therefore, the above transformation can be expressed as follows if the distortion is not considered:

[0122]

[0123] This is the mathematical model on which camera calibration is based.

[0124] In a non-coaxial camera, since the optical axis of the image plane and the optical axis of the lens do not coincide with each other but there is a deviation angle θ, the above model is not applicable in the process of converting the camera coordinate system to the image coordinate system. Therefore, the present application introduces the concepts of tilted image plane and non-tilted image plane, wherein the tilted image plane is the image plane perpendicular to the optical axis of the lens of the non-coaxial camera, and the non-tilted image plane is the imaging plane of the non-coaxial camera. In the process of converting the camera coordinate system to the image coordinate system, the camera coordinate system is first converted to the tilted image plane coordinate system, and then the tilted image plane coordinate system is converted to the non-tilted image plane coordinate system, and finally the non-tilted image plane coordinate system is converted to the image coordinate system. The tilt matrix H is used to transform the tilted image plane coordinate system to the non-tilted image plane coordinate system. tilt Therefore, if the distortion is not considered, the entire transformation process can be expressed as:

[0125]

[0126] Based on the above conversion model, this application provides a calibration method for a non-coaxial camera, please refer to Figure 4 In one embodiment, the method includes steps 110 to 150, which are described in detail below.

[0127] Step 110: Acquire a calibration plate image taken by a non-coaxial camera.

[0128] The calibration plate can be a checkerboard calibration plate, a circular array calibration plate, etc. When calibrating, the non-coaxial camera can be placed in multiple positions (positions are also the positions and angles of the non-coaxial camera relative to the calibration plate) based on experience, and the calibration plate is photographed at each position, thereby obtaining multiple different calibration plate images for calibration.

[0129] Step 120: Acquire feature points in the calibration plate image, as well as the image coordinates of the feature points and the corresponding world coordinates.

[0130] For the checkerboard calibration plate, the feature point is the corner point of the checkerboard, and for the circular array calibration plate, the feature point is the center of the circular feature point in the circular array, and the circular feature point is the circular pattern on the circular array calibration plate.

[0131] The world coordinate system can be constructed based on the parameter information of the calibration plate to obtain the world coordinates corresponding to the feature points. The parameter information of the calibration plate includes the size of the calibration plate, the size of the chessboard, the radius of the circular feature points, the spacing between the feature points, etc. The image coordinates of the feature points can be obtained by image processing the calibration plate image. The accurate extraction of the image coordinates of the feature points is crucial for camera calibration. It is necessary to take measures to improve the accuracy of feature point positioning. This application proposes a method for extracting high-precision coordinates of feature points for circular array calibration plates. The ellipse equation can be used to perform error correction on the image coordinates of the feature points, which effectively improves the accuracy of the image coordinates of the feature points, thereby improving the accuracy of camera calibration. This method will be described in detail below.

[0132] Step 130: Calculate the homography matrix H according to the image coordinates of the feature points and the corresponding world coordinates. The homography matrix H can be calculated using the image coordinates of multiple feature points and the corresponding world coordinates.

[0133] Step 140: According to the preset transformation model from world coordinates to image coordinates, the homography matrix H is decomposed and calculated to obtain the intrinsic parameters and extrinsic parameters of the non-coaxial camera, where the intrinsic parameters include the tilt matrix H tilt .

[0134] It can be understood that the homography matrix Therefore, the internal and external parameters of the non-coaxial camera can be obtained by decomposing the homography matrix H. tilt is regarded as part of the internal reference, then the internal reference part of the non-coaxial camera is The external reference part is

[0135] Please refer to Figure 5 , this application introduces three parameters to represent the tilt matrix H tilt , namely the image plane distance d, the rotation angle τ and ρ, the tilted image plane can be regarded as the translation distance d relative to the non-tilted image plane (i.e. the image plane of the non-coaxial camera), the rotation angle τ around the X-axis of the original coordinate system, and the rotation angle ρ around the Z-axis, where the X-axis of the original coordinate system is the horizontal direction of the non-tilted image plane, the Y-axis is the vertical direction of the non-tilted image plane, and the Z-axis is the vertical line of the non-tilted image plane. The rotation matrix Q is used to represent the rotation transformation of the tilted image plane relative to the original coordinate system, then the geometric transformation relationship can be calculated to get

[0136]

[0137] For image-side telecentric lenses and bilateral telecentric lenses, the tilt matrix H tilt for

[0138]

[0139] When the homography matrix H is decomposed to calculate the intrinsic and extrinsic parameters of the non-coaxial camera,

[0140]

[0141] As a whole, let H = A[R|t]. According to orthogonality, we can obtain:

[0142] H=[h1 h2 h3]=A[r1 r2 t],

[0143] Where [r1 r2 t] = [R|t], h1 is the first column vector of the homography matrix H, h2 is the second column vector of the homography matrix H, h3 is the third column vector of the homography matrix H, r1 is the first column vector of the rotation matrix R, and r2 is the second column vector of the rotation matrix R. The parameter matrix A can be calculated according to the following constraints:

[0144]

[0145] Then according to r1=A -1 h1, r2 = A -1 h2,t=A -1 h3 calculates the matrix [r1 r2 t] to obtain the external parameters. The internal parameters are the camera focal length f and the pixel size s. x and y , the principal optical axis point (c x ,c y ) can be known in advance and then based on The tilt matrix H can be calculated tilt .

[0146] In one embodiment, the rotation matrix R can be represented by an equivalent rotation axis k and an equivalent axis angle θ, and the equivalent rotation axis k and the equivalent axis angle θ can be regarded as part of the external parameters. The transformation relationship between the rotation matrix R and the equivalent rotation axis k and the equivalent axis angle θ is as follows:

[0147]

[0148] k x , k y , k z are the three components of the equivalent rotation axis k.

[0149] According to r1=A -1 h1, r2 = A -1 h2,t=A -1 After h3 calculates the matrix [r1 r2 t], we get the rotation matrix R. According to the above formula, we can calculate the equivalent rotation axis k and the equivalent axis angle θ.

[0150] Step 150: nonlinearly optimize the distortion coefficient of the non-coaxial camera and the decomposed internal parameters and external parameters to obtain the final internal parameters, external parameters and distortion coefficient of the non-coaxial camera.

[0151] In this step, the distortion coefficient of the non-coaxial camera and the decomposed internal and external parameters are nonlinearly optimized through the set loss function. The initial values ​​of the internal and external parameters in the iterative process can be the internal and external parameters obtained by the above decomposition, and the initial value of the distortion coefficient can be pre-set based on experience. Since the distortion occurs in the process of projecting the point through the lens to the tilted image plane, and the distortion is a nonlinear change, it is possible to choose to establish a loss function on the tilted image plane, and divide the entire transformation into two parts: one part is the transformation from the image coordinate system to the tilted image plane coordinate system, and the other part is the transformation from the world coordinate system to the tilted image plane. The closer the results of the two are, the better the calibration results are. Therefore, the loss function can be constructed as follows:

[0152]

[0153] This loss function can be adapted to the calibration of multiple cameras, where one camera is set as the reference camera, and the other cameras can be transformed into the coordinate system of the reference camera for unified calculation. m is the number of feature points in the calibration plate image, n c is the number of cameras, n0 is the number of calibration plate images taken by the camera, and p j is the coordinate of the feature point in the world coordinate system, e l (l=1,…,n0) represents the external parameters of the calibration plate image in the reference camera, r k (k=1,…,n c ) represents the transformation of the kth camera relative to the reference camera, i k Indicates that the transformation is the transformation under the kth camera, p jkl is the image coordinate of the jth feature point in the lth calibration plate image taken by the kth camera, v jkl The value is 0 or 1. It is 1 when the jth feature point is visible in the lth calibration plate image taken by the kth camera, otherwise it is 0. Function Represents the transformation from the image coordinate system to the tilted image plane coordinate system. From the above, we can see that this includes the process of using the intrinsic parameters to transform the image coordinate system to the non-tilted image plane coordinate system, transforming from the non-tilted image plane coordinate system to the tilted image plane coordinate system, and using the distortion coefficient to perform anti-distortion on the tilted image plane coordinate system. Function φ u (p j ,e l ,r k ,i k) represents the transformation from the world coordinate system to the tilted image plane coordinate system. From the above, we can see that this includes the process of using external parameters to transform the world coordinate system to the camera coordinate system and transform the camera coordinate system to the tilted image plane coordinate system.

[0154] The LM algorithm can be used for iterative calculation, and the update of parameters during the iteration process can be expressed as q k+1 =q k +δ, where q k represents the vector composed of the intrinsic parameters, extrinsic parameters and distortion coefficients of the non-coaxial camera at the kth iteration. δ is determined by the formula Jδ = ε, where ε is the number of feature points corresponding to each calibration plate image captured by each camera at the current iteration. With φ u (p j ,e l ,r k ,i k ), J is the Jacobian matrix, which is composed of the Jacobian matrices of each camera, and the Jacobian matrix of the i-th camera is

[0155]

[0156] Among them J i,j and J i ' ,j They represent the partial derivatives of the intrinsic and extrinsic parameters corresponding to the j-th calibration plate image captured by the ith camera. The solution of the partial derivatives is explained below.

[0157] For the external parameter part, the rotation matrix R is represented by the equivalent rotation axis and the equivalent axis angle. Let the rotation axis be [r x ,r y ,r z ] T , then the equivalent axis angle is equal to: Therefore, the unit rotation vector can be obtained as Then the rotation matrix R can be expressed as:

[0158] R=cosθ·I+(1-cosθ)rr T +sinθR z ,

[0159] Where I is the identity matrix, R z is an antisymmetric matrix, and

[0160] definition Available

[0161] a0=-sinθl i ,a1=[sinθ-2(1-cosθ)θ′]l i ,

[0162] a2=2(1-cosθ)θ′, a3=[cosθ-θ′sinθ]l i ,a4=θ′sinθ,

[0163] When i=0,l i = l x ; i=1,l i = l y ; i=2,l i = l z .

[0164] Define a vector:

[0165] i=[1,0,0,0,1,0,0,0,1],

[0166] dr0=[2l x ,l y ,l z ,l y ,0,0,l z ,0,0],dr1=[0,l x ,0,l x ,2l y ,l z ,0,l z ,0],dr2=[0,0,l x ,0,0,l y ,l x ,l y ,2l z ],

[0167] q x =[0,-r z ,r y ,r z ,0,-r x ,-r y ,r x ,0],

[0168] dq0=[0,0,0,0,0,-1,0,1,0], dq1=[0,0,1,0,0,0,-1,0,0], dq2=[0,-1,0,1,0,0,0,0,0],

[0169] Then the partial derivative of the extrinsic parameter of the i-th camera can be expressed as:

[0170] J′ i =[J″0,J″1,J″2] T ,

[0171] Where J″ i=a0i+a1r t +a2dr i +a3q x +a4dq i (i=0,1,2).

[0172] As for the internal parameter part, it can be seen from formula (2) that it can be divided into three parts for solution. The first part is Since formula (3) involves the transformation from the image coordinate system to the non-tilted image plane coordinate system, it is Taking the derivative, we can get

[0173]

[0174] According to formula (1), without considering the distortion,

[0175] u=(cc x )s x ,v=(rc y )s y ,

[0176] Therefore, the partial derivative can be obtained:

[0177]

[0178] The second part is the tilt matrix H tilt , the same thing needs to be done here To perform the derivation, H can be obtained by taking the partial derivative of the rotation matrix Q tilt For the rotation matrix Q, if we directly calculate it based on τ and ρ, due to the ambiguity of the rotation angle, the obtained τ and ρ are not the only solution, but two solutions that meet the conditions will appear. In order to eliminate the ambiguity, the constraint condition is added here:

[0179]

[0180] make Then the parameters in the rotation matrix Q can be expressed as:

[0181]

[0182] Therefore, we can get

[0183]

[0184] According to the formula

[0185]

[0186] The partial derivatives of t2 and c2 with respect to S and C can be obtained. Substituting them into the rotation matrix Q, the partial derivatives of the rotation matrix Q with respect to S and C can be obtained, and thus the tilt matrix H can be obtained. tilt The partial derivative of the inverse matrix is ​​obtained by using the differential matrix formula d(X -1 )=-X -1 (dX)X -1 , we can obtain The partial derivative of .

[0187] The third part is The same thing needs to be done here For CCTV lens, The partial transformation can be expressed as:

[0188]

[0189] Then we can get:

[0190] For a telecentric lens, The partial transformation can be expressed as:

[0191]

[0192] Then we can get:

[0193] For the distortion coefficient part, for the division model, then

[0194]

[0195] For the polynomial model, Therefore, we can get:

[0196]

[0197] Write the derivative variables as a vector representation ν = [k1, k2, k3, k4, k5, k6, p1, p2] T , then the above formula can be expressed as:

[0198]

[0199] Combining the above parts, the partial derivative of the entire internal parameter part can be obtained according to the chain rule.

[0200] Since distortion is not considered when calculating the intrinsic parameters and extrinsic parameters in step 140, distortion correction can be added during the nonlinear optimization process. In one embodiment, during the nonlinear optimization process, the calculated distortion coefficients are used to perform distortion correction on the tilted image plane coordinates of the feature points before each iteration. Specifically, for the division model, since a reversible solution can be directly obtained, it can be directly calculated according to the following formula:

[0201]

[0202] For the polynomial model, it is assumed that the distortion model can be expressed as in is the distorted tilted image plane coordinate, p c is the coordinate in the camera coordinate system, the vector f d By u and f v It consists of two parts, and

[0203]

[0204] Therefore, the corrected coordinates can be expressed as In order to calculate f d It can be Taylor expanded to consider only the linear part:

[0205]

[0206] Therefore, we can get

[0207]

[0208] Therefore, the process of distortion correction for the polynomial model is as follows:

[0209] (1) The image coordinates (r, c) T Transform to oblique image plane coordinates

[0210] (2) Iterative calculation to eliminate distortion:

[0211] initialization calculate and Then update x and y according to the following formula and iterate:

[0212] Δx=p1(r 2 +2u 2 )+2p2uv

[0213] Δy=2p1uv+p2(r 2 +2v 2 )

[0214]

[0215] According to the calibration method of the non-coaxial camera of the above embodiment, for the non-coaxial camera, when establishing the transformation model from the world coordinate to the image coordinate, the situation that the optical axis of the lens of the non-coaxial camera is not coaxial with the optical axis of the imaging plane is taken into consideration, and the concepts of tilted image plane and non-tilted image plane are introduced, wherein the tilted image plane is an image plane perpendicular to the optical axis of the lens of the non-coaxial camera, and the non-tilted image plane is the imaging plane of the non-coaxial camera, and the tilt matrix H is added. tilt This parameter is used to describe the transformation from the tilted image plane coordinate system to the non-tilted image plane coordinate system. The tilt matrix H tilt The method is regarded as a part of the camera internal parameters, and the mathematical model in the existing method is corrected; and the rotation angle τ around the X-axis of the original coordinate system and the rotation angle ρ around the Z-axis are introduced to represent the transformation between the tilted image plane and the non-tilted image plane, which can well describe the coordinate system transformation process in the non-coaxial camera; in one embodiment, the tilted image plane coordinates are also distorted during the nonlinear optimization process. In summary, the calibration method of the non-coaxial camera provided in the present application improves the calibration accuracy of the non-coaxial camera and reduces the error of the working results generated in the process of applying the non-coaxial camera.

[0216] The following describes the method for extracting high-precision coordinates of feature points for the circular array calibration plate mentioned in step 120. Please refer to Figure 6 In one embodiment, the method includes steps 210 to 250.

[0217] Step 210: Acquire a calibration plate image. It should be noted that the acquired calibration plate image may be taken by a coaxial camera or a non-coaxial camera.

[0218] Step 220: Perform image processing on the calibration plate image to obtain circular feature points therein. Image processing includes binarization, filtering, feature screening, etc. Please refer to Figure 7 In one embodiment, the process of obtaining circular feature points includes steps 310 to 340, which are described in detail below.

[0219] Step 310: extract edges of the calibration plate image to obtain a calibration plate boundary box, thereby obtaining the position of the calibration plate in the calibration plate image.

[0220] Step 320: construct an image pyramid for the calibration plate area to obtain pyramid images of each layer, wherein the calibration plate area is the area in the calibration plate image that is within the calibration plate boundary box. In the image pyramid, the image resolution of the upper layer is smaller, and the image resolution of the lower layer is larger. The specific number of layers of the image pyramid can be set based on experience.

[0221] Step 330: Binarize the current layer pyramid image to search for circular feature points. The initial value of the current layer pyramid image is the topmost layer pyramid image.

[0222] In one embodiment, the binarization process may be a process that is iterated according to a gray value step length. Specifically, within a preset threshold interval, the gray threshold is selected from small to large intervals of the preset interval value. Each time the gray threshold is selected, the current layer pyramid image is threshold segmented using the gray threshold to obtain a circular area. When the number of circular areas is equal to the preset number, it is determined that a circular feature point that meets the preset conditions is searched in the current layer pyramid image, and the next gray threshold is stopped from being selected. Otherwise, the next gray threshold is continuously selected to perform threshold segmentation on the current layer pyramid image until the preset threshold interval is traversed. For example, if the preset threshold interval is 50-90, if the step length is 10, that is, the preset interval value is 10, then 50, 60, 70, 80, and 90 are selected in turn as gray thresholds to perform threshold segmentation on the current layer pyramid image until the number of circular areas is equal to the preset number. The preset threshold interval can be set based on experience. After the circular area is obtained by threshold segmentation, some morphological processing and area screening can also be performed to obtain more accurate results.

[0223] Step 340: Determine whether a circular feature point that meets the preset conditions is found in the current layer pyramid image. If so, the process ends; otherwise, the next layer pyramid image is used as the current layer pyramid image and the process returns to step 330.

[0224] Please refer to Figure 8 ,In another embodiment, the process of obtaining circular feature points includes steps 410 to 430, which are described in detail below.

[0225] Step 410: construct an image pyramid for the calibration plate image to obtain pyramid images of each layer.

[0226] Step 420: Search for circular feature points on the current layer pyramid image. The initial value of the current layer pyramid image is the topmost layer pyramid image.

[0227] The circular feature point search can be performed in the following manner: binarize the current layer pyramid image to obtain a circular area, perform statistical analysis on the area of ​​the circular area to obtain the area with the highest frequency of occurrence, calculate the radius based on the area with the highest frequency of occurrence, multiply it by the magnification corresponding to the current layer pyramid image to obtain the estimated radius of the circular feature point, and search for the circular feature point in the calibration plate image based on the estimated radius.

[0228] The estimated radius of the circular feature point can be obtained through histogram statistics. After binarizing the current layer pyramid image to obtain the circular area, the circular area is first screened according to the preset roundness range and / or area range, and the area histogram statistics of the screened circular area are performed to establish a functional mapping relationship between the area and the frequency of occurrence to obtain the area with the highest frequency of occurrence; then the radius is calculated according to the area with the highest frequency of occurrence, and multiplied by the magnification corresponding to the current layer pyramid image to obtain the estimated radius of the circular feature point; finally, the calibration plate image is filtered according to the estimated radius and then threshold segmentation is performed to obtain the feature point estimation area, and the feature point estimation area with an area greater than the preset area threshold is eliminated. The number of feature point estimation areas is calculated. When the number of feature point estimation areas is equal to the preset number, it is determined that the circular feature point that meets the preset conditions has been searched in the current layer pyramid image.

[0229] Step 430: Determine whether a circular feature point that meets the preset conditions is found in the current layer pyramid image. If so, the process ends; otherwise, the next layer pyramid image is used as the current layer pyramid image and the process returns to step 420.

[0230] The method for obtaining circular feature points in the calibration plate image in the above embodiment searches for circular feature points by constructing an image pyramid, starting from the top layer of the image pyramid and searching layer by layer to the lower layers. When a circular feature point that meets the preset conditions is found in a certain layer, the search can be stopped. A coarse-to-fine search strategy is adopted, and the original image is not directly used for processing. Since the image resolution of the upper layer of the image pyramid is small and the image is small, this method is conducive to improving the search efficiency of circular feature points. In some embodiments, the binarization process is an iterative process according to a gray value step size, and a single gray threshold is not used for segmentation, which is conducive to more accurate extraction of circular feature points.

[0231] The following will continue to introduce steps 230 to 250.

[0232] Step 230: extracting edges of circular feature points to obtain edge points of circular feature points, and using the edge points to perform ellipse fitting to obtain image coordinates of the circular feature points, wherein the image coordinates of the circular feature points refer to the image coordinates of the center of the circular feature points.

[0233] Step 240: Determine the correspondence between the image coordinates and the world coordinates of the circular feature points.

[0234] Commonly used calibration plates often do not have a reference. Workers need to manually select the reference to compare the image coordinates of the feature points with the world coordinates to determine their correspondence, which is cumbersome. This application provides two circular array calibration plates with references, and gives methods for determining the correspondence between the image coordinates of circular feature points and the world coordinates. One of them is a circular array calibration plate with triangular markers, such as Fig. 9 As shown, a triangular marker is set at one corner of the calibration plate. The triangular marker is an isosceles right triangle, and its right-angled vertex is one of the vertices of the circular array calibration plate, and the other two vertices are respectively on the two sides of the circular array calibration plate adjacent to the right-angled vertex. Another type is a circular array calibration plate with hollow points, such as Fig.10 As shown, the hollow points are aggregated into several clusters, such as Fig.10 There are 5 clusters in .

[0235] Please refer to Fig.11 , in the circular array calibration plate with triangular markers, determining the correspondence between the image coordinates and the world coordinates of the circular feature points includes the following steps:

[0236] Step 510: Detect the triangular marker in the calibration plate image and determine the relative position relationship between the circular feature point and the triangular marker. The triangular marker can be detected by detecting the hypotenuse of the triangle.

[0237] Step 520: Establish a reference coordinate system based on the triangular marker, and determine the one-to-one correspondence between the image coordinates and the world coordinates of the circular feature points according to the parameter information of the circular array calibration plate. On the basis of step 510, after establishing the reference coordinate system, the position of the circular feature points in the reference coordinate system can be obtained, and the reference coordinates and the world coordinates can be matched. The one-to-one correspondence between the image coordinates and the world coordinates of the circular feature points can be determined using the parameter information of the circular array calibration plate.

[0238] Please refer to Fig.12 , in the circular array calibration plate with hollow points, determining the correspondence between the image coordinates and the world coordinates of the circular feature points includes the following steps:

[0239] Step 610: extract hollow points from the obtained circular feature points, and use a clustering algorithm to divide the hollow points into different clusters.

[0240] Step 620: Calculate the hollow point in the cluster with the shortest sum of distances from all other hollow points, use it as the center point of the cluster, and assign the non-hollow points to the cluster with the shortest distance to it.

[0241] Step 630: Determine the position of the cluster in the circular array calibration plate according to the arrangement of the hollow points in the cluster. Fig.10In the figure, we can see that the arrangement of the hollow points in the five clusters is different, which can be used to determine the position of the cluster in the circular array calibration plate.

[0242] Step 640: Taking one of the clusters as a reference cluster, the relative positional relationship between the other clusters and the reference cluster is determined. In this way, the relative positional relationship between the circular feature points in the other clusters and the reference cluster can also be determined.

[0243] Step 650: Establish a reference coordinate system with the center point of the reference cluster as the origin, and determine the one-to-one correspondence between the image coordinates and the world coordinates of the circular feature points according to the parameter information of the circular array calibration plate. On the basis of step 640, after establishing the reference coordinate system, the position of the circular feature points in the reference coordinate system can be obtained, and the reference coordinates and the world coordinates can be matched. The one-to-one correspondence between the image coordinates and the world coordinates of the circular feature points can be determined using the parameter information of the circular array calibration plate.

[0244] In one embodiment, after obtaining the correspondence between the image coordinates and the world coordinates of the circular feature points, sub-pixel edge extraction can also be performed. Specifically, first, a homography matrix is ​​calculated based on the one-to-one correspondence between the image coordinates and the world coordinates of the hollow points in the cluster, and the world coordinates of other circular landmark points are mapped onto the image using the homography matrix to obtain mapping points; then, the circular feature points containing the mapping points are obtained, and sub-pixel edge extraction and ellipse fitting are performed on them to obtain new edge points and image coordinates of the circular feature points. The accuracy of the edge points and image coordinates obtained using sub-pixel edge extraction is further improved.

[0245] Step 250: Use the ellipse equation to perform error correction on the image coordinates of the circular feature points to obtain the final image coordinates of the circular feature points.

[0246] For a circle with a radius of r and a center at (X0, Y0), its equation can be expressed as: T Fx = 0, where

[0247]

[0248] F is called the ellipse equation matrix. The center of the circle can be expressed as: c conic =F -1 (0,0,1) T The transformation from the world coordinate system to the image coordinate system can be expressed as: i =H t p w , so the transformed circle can be expressed as: The transformed center can be expressed as:

[0249] The above transformation from the world coordinate system to the image coordinate system does not take distortion into account. If distortion is taken into account, the relationship between distortion and non-distortion needs to be established. In one embodiment of the present application, the above ellipse equation matrix is ​​used to establish an objective function based on the idea of ​​minimizing the gap between the observed value and the expected value, and the non-distorted image coordinates are solved to achieve error correction of the image coordinates of the circular feature points.

[0250] The ellipse before error correction is a distorted ellipse, which can be expressed as: in is the image coordinate of the circular feature point before error correction, is the distorted ellipse equation matrix. If the distortion is not considered, it is a standard quadratic curve. The ellipse equation matrix can be expressed as D, and the equation of the curve is p T Dp = 0, where p is the image coordinate of the circular feature point after error correction. The mapping can be used with the transformation matrix H D To express:

[0251]

[0252] Diagonalizing the elliptic equation matrix yields: U T DU=Λ, in

[0253]

[0254] Where λ1, λ2 and λ3 are the eigenvalues ​​of the elliptic equation matrix D, and U is the matrix consisting of the corresponding eigenvectors. and is the elliptic equation matrix The characteristic value of is the matrix consisting of the corresponding eigenvectors.

[0255] make but

[0256] Get H D After that, the image coordinates p of the circular feature point after error correction can be solved according to the following objective function: i :

[0257]

[0258] The subscript i represents the i-th point.

[0259] In another embodiment, the ellipse equation can be obtained by ellipse fitting to obtain the image coordinates of the ellipse center point, which are compared with the image coordinates of the ellipse center point without fitting, and the difference between the two is used to directly perform error correction on the image coordinates of the circular feature point. The ellipse fitting can be performed according to the following objective function:

[0260]

[0261] st4ac-b 2 =1

[0262] Where a, b, c, d, e, and f are the coefficients of the ellipse equation, (x i ,y i ) is the image coordinate of the edge point of the circular feature point, w i is the weight, and n is the number of edge points.

[0263] Then calculate the image coordinates of the ellipse center point (r d ,c d ):

[0264]

[0265] The image coordinates of the center point of the unfitted ellipse (r d ′,c′ d ):

[0266]

[0267] Where F is the set of points in the area where the circular feature point is located, that is, all the points in the entire circle, p i is a point in the set, I(p i ) is point p i The gray value, (x i ,y i ) is point p i The image coordinates of , where the subscript i represents the i-th point.

[0268] Calculate the image coordinates (r d ,c d ) and (r d ′,c′ d ) is used to compensate the image coordinates of the circular feature points and complete the error correction.

[0269] According to the method for extracting the high-precision coordinates of the feature points of the circular array calibration plate of the above embodiment, the calibration plate image is processed to obtain the circular feature points therein, and then the circular feature points are subjected to edge extraction and ellipse fitting to obtain the image coordinates of the circular feature points. After obtaining the image coordinates, the ellipse equation is used to correct the errors. The error correction can be achieved by ellipse fitting, error compensation, etc. In the process of obtaining the circular feature points, the circular feature points are searched by constructing an image pyramid, and the search is performed layer by layer from the top layer of the image pyramid to the lower layer. When a circular feature point that meets the preset conditions is found in a certain layer, it can be stopped. A coarse-to-fine search strategy is adopted, and the original image is not directly used for processing. Since the image resolution of the upper layer of the image pyramid is small and the image is small, this method is conducive to improving the search efficiency of the circular feature points. In some embodiments, the binarization process is an iterative process according to a gray value step size, and a single gray threshold is not used for segmentation, which is conducive to more accurate extraction of circular feature points. In summary, the method for extracting the high-precision coordinates of the feature points of the circular array calibration plate provided by the present application effectively improves the accuracy and efficiency of feature point coordinate extraction, thereby improving the accuracy of camera calibration.

[0270] Those skilled in the art will appreciate that all or part of the functions of the various methods in the above-mentioned embodiments can be implemented by hardware or by computer programs. When all or part of the functions in the above-mentioned embodiments are implemented by computer programs, the program can be stored in a computer-readable storage medium, and the storage medium can include: read-only memory, random access memory, disk, optical disk, hard disk, etc., and the program is executed by a computer to implement the above-mentioned functions. For example, the program is stored in the memory of the device, and when the program in the memory is executed by the processor, all or part of the above-mentioned functions can be implemented. In addition, when all or part of the functions in the above-mentioned embodiments are implemented by computer programs, the program can also be stored in a storage medium such as a server, another computer, disk, optical disk, flash disk or mobile hard disk, and can be downloaded or copied and saved in the memory of the local device, or the system of the local device is updated, and when the program in the memory is executed by the processor, all or part of the functions in the above-mentioned embodiments can be implemented.

[0271] The above specific examples are used to illustrate the present invention, which is only used to help understand the present invention and is not intended to limit the present invention. For those skilled in the art, according to the concept of the present invention, some simple deductions, modifications or substitutions can be made.

Claims

1. A calibration method for a non-coaxial camera, wherein the non-coaxial camera comprises an image plane and a lens, wherein the normal vector of the image plane is not coaxial with the optical axis of the lens, and the lens is an image-space telecentric lens or a bilateral telecentric lens, wherein: The calibration method comprises: Obtain a calibration plate image taken by a non-coaxial camera; Acquire feature points in the calibration plate image, as well as image coordinates and corresponding world coordinates of the feature points; Calculate the homography matrix H according to the image coordinates of the feature points and the corresponding world coordinates; According to the preset transformation model from world coordinates to image coordinates, the homography matrix is ​​decomposed and calculated to obtain the intrinsic parameters and extrinsic parameters of the non-coaxial camera. The transformation model is: Among them, the homography matrix (r,c) T is the image coordinate of the feature point, (x w ,y w ,z w ) T is the world coordinate of the feature point, is the transformation matrix from the world coordinate system to the camera coordinate system, R is the rotation matrix, t is the displacement matrix, z c is the z coordinate of the feature point in the camera coordinate system, is the transformation matrix from the camera coordinate system to the tilted image plane coordinate system, f is the focal length of the non-coaxial camera, H tilt is a tilt matrix, representing the transformation from the tilted image plane coordinate system to the non-tilted image plane coordinate system, wherein the tilted image plane is the image plane perpendicular to the optical axis of the lens, and the non-tilted image plane is the image plane of the non-coaxial camera, is the transformation matrix from the non-tilted image plane coordinate system to the image coordinate system, s x and y are the pixel sizes of the non-coaxial camera in the horizontal and vertical directions, respectively. x ,c y ) is the principal optical axis point, For the internal reference part, is the external parameter part, the tilt matrix H tilt Specifically Among them, q 11 ,q 12 ,q 21 ,q 22 is an element in the rotation matrix Q, which represents the rotation transformation of the tilted image plane relative to the original coordinate system, and Wherein ρ represents the angle of rotation around the Z axis, τ represents the angle of rotation around the X axis, the X axis of the original coordinate system is the horizontal direction of the non-tilted image plane, the Y axis is the vertical direction of the non-tilted image plane, and the Z axis is the vertical line of the non-tilted image plane; The distortion coefficient of the non-coaxial camera and the decomposed internal and external parameters are nonlinearly optimized to obtain the final internal and external parameters and distortion coefficient of the non-coaxial camera.

2. The calibration method according to claim 1, characterized in that: The method of decomposing the homography matrix to obtain the intrinsic parameters and extrinsic parameters of the non-coaxial camera includes: Calculate the parameter matrix A according to the following constraints in H=[h1 h2 h3]=A[r1 r2 t], [r1 r2 t] = [R|t]; Where h1 is the first column vector of the homography matrix H, h2 is the second column vector of the homography matrix H, h3 is the third column vector of the homography matrix H, r1 is the first column vector of the rotation matrix R, and r2 is the second column vector of the rotation matrix R; According to r1=A -1 h1, r2 = A -1 h2,t=A -1 h3 calculates the matrix [r1 r2 t], according to Calculate the tilt matrix H tilt .

3. The calibration method according to claim 1, characterized in that: The external parameters of the non-coaxial camera also include an equivalent rotation axis k and an equivalent axis angle θ. The internal parameters and external parameters of the non-coaxial camera obtained by decomposing and calculating the homography matrix include: Calculate the parameter matrix A according to the following constraints in, H=[h1 h2 h3]=A[r1 r2 t], [r1 r2 t] = [R|t]; According to r1=A- 1 h1, r2 = A- 1 h2,t=A- 1 h3 calculates the matrix [r1 r2 t]; Where h1 is the first column vector of the homography matrix H, h2 is the second column vector of the homography matrix H, h3 is the third column vector of the homography matrix H, r1 is the first column vector of the rotation matrix R, and r2 is the second column vector of the rotation matrix R; According to the calculated rotation matrix R, the equivalent rotation axis k and the equivalent axis angle θ are obtained, where the transformation relationship between the rotation matrix R and the equivalent rotation axis k and the equivalent axis angle θ is as follows: kx, ky, kz are the three components of the equivalent rotation axis k; according to The tilt matrix Htilt is calculated.

4. The calibration method according to claim 1, characterized in that: The nonlinear optimization of the distortion coefficient of the non-coaxial camera and the decomposed internal parameters and external parameters to obtain the final internal parameters, external parameters and distortion coefficient of the non-coaxial camera includes: The initial value of the distortion coefficient is set in advance, and the decomposed internal and external parameters are used as the initial values ​​of the internal and external parameters. The optimal solution is iteratively solved according to the following loss function to obtain the final internal and external parameters and distortion coefficient of the non-coaxial camera: Where, nm is the number of feature points in the calibration plate image, nc is the number of cameras, n0 is the number of calibration plate images taken by the camera, pj is the coordinate of the feature point in the world coordinate system, el(l=1,…,n0) represents the external parameter of the calibration plate image in the reference camera, rk(k=1,…,nc) represents the transformation of the kth camera relative to the reference camera, ik represents the transformation under the kth camera, pjkl is the image coordinate of the jth feature point in the lth calibration plate image taken by the kth camera, vjkl takes the value of 0 or 1, when the jth feature point is visible in the lth calibration plate image taken by the kth camera, it is 1, otherwise it is 0; Function It represents the transformation from the image coordinate system to the tilted image plane coordinate system, which includes the use of intrinsic parameters to transform the image coordinate system to the tilted image plane coordinate system and the use of distortion coefficients to dedistort the tilted image plane coordinates. The function φu(pj,el,rk,ik) represents the transformation from the world coordinate system to the tilted image plane coordinate system, which includes the use of extrinsic parameters to transform the world coordinate system to the camera coordinate system.

5. The calibration method according to claim 4, characterized in that: Also includes: Before each iteration, the calculated distortion coefficients are used to perform distortion correction on the tilted image plane coordinates of the feature points.

6. The calibration method according to claim 4, characterized in that: According to the formula q k+1 =q k +δ iterates to find the optimal solution, where q k represents the vector composed of the intrinsic parameters, extrinsic parameters and distortion coefficients of the non-coaxial camera at the kth iteration. δ is determined by the formula Jδ = ε, where ε is the number of feature points corresponding to each calibration plate image captured by each camera at the current iteration. With φ u (p j ,e l ,r k ,i k ), J is the Jacobian matrix, which is composed of the Jacobian matrices of each camera, and the Jacobian matrix of the i-th camera is Among them J i,j and J i ' ,j They respectively represent the partial derivatives of the intrinsic and extrinsic parameters corresponding to the j-th calibration plate image taken by the ith camera.

7. The calibration method according to claim 1, characterized in that: The calibration plate image is a circular array calibration plate image; the feature points in the calibration plate image and the image coordinates of the feature points are obtained by: Performing image processing on the calibration plate image to obtain circular feature points therein; Perform edge extraction on the circular feature points to obtain edge points of the circular feature points, and use the edge points to perform ellipse fitting to obtain image coordinates of the circular feature points, wherein the image coordinates of the circular feature points refer to the image coordinates of the center of the circular feature points; Determine the correspondence between the image coordinates and the world coordinates of the circular feature points; The image coordinates of the circular feature points are corrected using the ellipse equation to obtain the final image coordinates of the circular feature points.

8. The calibration method according to claim 7, characterized in that: The method of using the ellipse equation to perform error correction on the image coordinates of the circular feature points includes: Calculate the non-distorted ellipse equation matrix D to the distorted ellipse equation matrix The transformation matrix H D , where the undistorted elliptic equation matrix D satisfies p T Dp=0, Distorted ellipse equation matrix satisfy Transformation matrix H D satisfy Where p is the image coordinate of the circular feature point after error correction, are the image coordinates of the circular feature points before error correction, and where U, Λ0 and U satisfy Where λ1, λ2 and λ3 are the eigenvalues ​​of the elliptic equation matrix D, and U is the matrix consisting of the corresponding eigenvectors. and is the elliptic equation matrix The characteristic value of is the matrix consisting of the corresponding eigenvectors; According to the following objective function, the image coordinates p of the circular feature point after error correction are obtained. i : The subscript i represents the i-th point.

9. The calibration method according to claim 7, characterized in that: The method of using the ellipse equation to perform error correction on the image coordinates of the circular feature points includes: According to the following objective function, the edge points of the circular feature points are used to fit the ellipse and obtain the coefficients in the ellipse equation: Where a, b, c, d, e, and f are the coefficients of the ellipse equation, (x i ,y i ) is the image coordinate of the edge point, w i is the weight, n is the number of edge points; The image coordinates of the center point of the ellipse (r) are calculated according to the following formula d ,c d ): The image coordinates of the center point of the unfitted ellipse (r d ′,c′ d ): Where F is the set of points in the area where the circular feature point is located, pi is a point in the set, and I(p i ) is point p i The gray value, (x i ,y i ) is point p i The image coordinates of , where the subscript i represents the i-th point; According to the image coordinates (r d ,c d ) and (r d ′,c′ d ) and perform error correction on the image coordinates of the circular feature points.

10. A computer-readable storage medium, characterized in that: The medium stores a program, which can be executed by a processor to implement the calibration method according to any one of claims 1 to 9.

Images