CN104374338A

CN104374338A - Single-axis rotation angle vision measurement method based on fixed camera and single target

Info

Publication number: CN104374338A
Application number: CN201410506279.0A
Authority: CN
Inventors: 全权; 董洪信; 张瑞峰; 付强; 蔡开元
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2014-09-28
Filing date: 2014-09-28
Publication date: 2015-02-25
Anticipated expiration: 2034-09-28
Also published as: CN104374338B

Abstract

The invention provides a single-axis rotation angle vision measurement method based on a fixed camera and a single target. The method includes the four steps of firstly, calibrating internal parameters and external parameters of the camera; secondly, obtaining the initial values of the external parameters related to pivoting movement; thirdly, conducting nonlinear optimization on the internal parameter pin and the external parameter pout; fourthly, measuring the rotation angle through the calibrated internal parameters and the calibrated external parameters. The method is a method for measuring the rotation angle of the pivoting movement through the monocular camera, and solves the problem that the target is hard to install in an existing angle vision measurement method; meanwhile, the method is suitable for the conditions with blockage as well and is quite high in practicability.

Description

Visual measurement method of single-axis rotation angle based on fixed camera and single target

Technical Field

The invention relates to a vision measurement method of a single-axis rotation angle based on a fixed camera and a single target, and belongs to the technical field of vision measurement.

Background

The angular displacement measurement is widely applied to the fields of industry, aerospace and the like. Because the electrical and mechanical parameters of the angular displacement sensor change in a severe environment, the angular displacement sensor needs to be calibrated and corrected periodically. The vision measurement technology is favored because of the advantages of simple structure, non-contact real-time measurement, convenient field installation and debugging, etc. The existing visual angle measuring methods are divided into binocular vision and monocular vision. The binocular vision method is greatly limited due to small common visual field and complex installation. The existing monocular vision method has certain requirements on the installation of a target and has high installation difficulty. The present patent application therefore proposes a visual measurement method based on a single axis rotation angle of a fixed camera and a single target, which does not have any requirements for the mounting of the target. In addition, the invention is also true for the case of shielding, and is more practical.

Disclosure of Invention

The invention provides a visual measurement method of a single-axis rotation angle based on a fixed camera and a single target, which is a method for measuring the rotation angle of the rotation motion around an axis by using a monocular camera. The method solves the problem of difficult target installation in the existing visual angle measurement method. Meanwhile, the invention is also applicable to the shielding condition and has good practicability.

The method calibrates the camera internal parameters and the external parameters relative to the axis through the two-dimensional target, and further measures the rotation angle.

The camera model adopted by the invention is a nonlinear perspective projection model, and the model is described as follows:

1. linear camera model

As shown in fig. 1, the imaging position of any point P in space in the image can be approximately represented by a pinhole imaging model, i.e. the projection position P of the point P in the image is the intersection point of the line connecting the optical center O and the point P and the image plane. Thus the coordinate (X) of point P in the world coordinate system_w,Y_w,Z_w)^TWith the pixel coordinate (u, v) of the projection point p^TThe relationship between them is as follows:

wherein alpha is_xf/dX is a scale factor on the u-axis, alpha_yf/dY is the scale factor on the v-axis. Alpha is alpha_x,α_y,u₀,v₀Only the camera intrinsic parameters are related, which are called the intrinsic parameters of the camera. R_cw,T_cwThe rotation matrix and the translation vector of the camera coordinate system and the world coordinate system are respectively called external parameters of the video camera.

1. Non-linear camera model

In practice, the lens is not ideal perspective imaging, but rather with some degree of distortion. The invention uses a non-linear distortion model, consisting of world coordinates (X)_w,Y_w,Z_w)^TFinding the coordinates of a pixel (u, v)^TThe process of (2) is as follows:

[\begin{matrix} X_{c} \\ Y_{c} \\ Z_{c} \end{matrix}] = R_{cw} [\begin{matrix} X_{w} \\ Y_{w} \\ Z_{w} \end{matrix}] + T_{cw}, - - - (1.2)

x＝X_c/Z_c,y＝Y_c/Z_c,r²＝x²+y², (1.3)

x_{d} = [\begin{matrix} (1 + k_{1} r^{2} + k_{2} r^{4} + k_{5} r^{6}) x + {2 k}_{3} xy + k_{4} (r^{2} + {2 x}^{2}) \\ (1 + k_{1} r^{2} + k_{2} r^{4} + k_{5} r^{6}) y + {2 k}_{4} xy + k_{3} (r_{2} + {2 y}^{2}) \end{matrix}], - - - (1.4)

wherein k is₁,k₂,k₃,k₄,k₅Is the distortion coefficient. Alpha is alpha_x,α_y,u₀,v₀,k₁,k₂,k₃,k₄,k₅Forming an intrinsic parameter, R, of the camera_cw,T_cwBeing an external parameter of the camera.

For convenience of description, we remember u and v as p_in＝(α_x,α_y,u₀,v₀,k₁,k₂,k₃,k₄,k₅)^T,R_cw,T_cw,X_w,Y_w,Z_wA function of, i.e.

u＝f_u(p_in,R_cw,T_cw,X_w,Y_w,Z_w),v＝f_v(p_in,R_cw,T_cw,X_w,Y_w,Z_w) (1.6)

The world coordinate system established by the invention is as shown in figure (2), and the rotating shaft is taken as Z_wAxis and such that the origin of the target is located at X in the reference position_wOn the shaft. Assuming that the rotation matrix and the translation vector between the reference position target coordinate system and the world coordinate system areAnd setting the rotation matrix of the camera coordinate system and the world coordinate system as R_cwTranslation vector is T_cw.

The method provided by the invention comprises the following two steps: and calibrating internal and external parameter information and solving at any angle. The calibration of the internal and external parameter information is a process of solving the required external parameter information by utilizing the special geometric relationship between the rotating coordinate systems around the axis on the basis of the calibration of the common camera.t,R_cw,T_cwAnd the included angle between each position and the reference position is an external parameter to be estimated by calibrating the internal and external parameter information. The arbitrary angle solving means that the rotation angle is obtained according to the obtained internal and external parameter information. At any angleIn the degree solving process, a point matching algorithm is provided, so that the angle measuring method can solve the shielding problem.

In summary, the invention provides a visual measurement method based on a single-axis rotation angle of a fixed camera and a single target, which comprises the following specific steps:

the method comprises the following steps: and calibrating internal and external parameters of the camera.

In a specific implementation, we use a two-dimensional checkerboard (as in fig. 5) as the target. The sequence of pictures we use includes one picture at a reference position and m pictures at different angles of rotation of the rotating member. Suppose that the m angles are θ₁,θ₂,...,θ_m. Then, a checkerboard corner extraction function detectChecker-boardpoings carried by MATLAB2014 can be used for corner extraction, and an MATLAB calibration tool is combined for internal reference calibration to obtain an internal parameter p_inAnd external parameters for m +1 poses. The calibration tool kit can be downloaded on http:// www.vision.caltech.edu/bouguetj/caliib _ doc/. These extrinsic parameters may be expressed as

[R_{{cb}_{i}}, T_{{cb}_{i}}], i = 0,1,2, . . ., m .

Expressed at a rotation angle of theta_iWhen i is 0, the reference position is represented, the rotation matrix and the translation vector between the camera and the target.

Step two: initial value of external parameter related to pivoting motion

1. From the external parameters obtained in step one, the angle of rotation θ is obtained from (1.7)_iThe rotational and translational vectors of the target relative to the reference

[R_{b_{i} b_{0}}, T_{b_{i} b_{0}}], i = 1,2, . . ., m .

R_{b_{i} b_{0}} = R_{{cb}_{i}}^{T} R_{{cb}_{0}}, T_{b_{i} b_{0}} = R_{{cb}_{i}}^{T} (T_{{cb}_{0}} - T_{{cb}_{i}}) . - - - (1.7)

2. Obtainingt and theta_iInitial value:

definition of r₁,r₂,r₃Are respectively as1,2,3, and vector [ t ]]_×Is expressed as the following matrix

Then we can find out

Wherein,S_i,jrepresenting the elements of the ith row and j columns of the matrix S. Then, r₁Can be obtained from the following formula

The above formula can be passed throughObtained by SVD decomposition if

Then x is the third column of V. Then we can get

θ_i＝2arctan((A^TA)^-1(A^TB))， (1.12)

Wherein

A = vec (r_{2} r_{3}^{T} + r_{3} r_{2}^{T}),

Vec (-) here represents the straightening of the matrix. Therefore, the first and second electrodes are formed on the substrate,

R_{{wb}_{0}} = R_{b_{0} w}^{T} - - - (1.13)

3. according to the obtainedt and theta_iSolving the initial value, and solving R according to (1.16)_cw,T_cwValue of (A)

R_{cw} = R_{{cb}_{0}} R_{{wb}_{0}}^{T}, T_{cw} = {(t, 0,0)}^{T} - R_{cw} T_{{wb}_{0}} . - - - (1.15)

Step three: intra-pair parameter p_inAnd an extrinsic parameter p_outAnd carrying out nonlinear optimization. Where p is_in∈□⁹,p_out∈□¹⁰Is shown below

Wherein,are respectively asR_cwThe corresponding Rodrigues vector. The conversion relationship between a Rodrigues vector v and its corresponding rotation matrix R is as follows

Wherein, I₃An identity matrix of order 3 is represented.

We label the coordinates of each image point, on target 1, at a rotation angle θ_iThen, the corner point coordinate corresponding to the jth target feature pointIs marked as

(u_{j}^{i}, v_{j}^{i}), i = 0,1, . . ., m, j = 1,2, . . ., n,

Points on the target are noted

X_b＝{(x_j，y_j，0)^T，j＝1，2，...，n}.

Because errors exist in the image extraction process, the reprojection error sum of the images is selected as an optimization target. From the target point and p_outObtaining the reprojection pointsThe process of (2) is as follows:

1. angle of rotation theta_iTime, target pointConversion to points in world coordinate system

2. Finding the corresponding reprojection point

The reprojection error is defined as

The optimization problem can be described as

The optimization problem can be solved by utilizing a sparse Levenberg-Marquardt algorithm, and an optimization tool box of the sparse LM algorithm can be downloaded on http:// users. ics. forth. With the above optimization we can obtaint,R_cw,T_cwAnd (3) the information is used for the rotation angle calculation in the step four.

Step four: and measuring the rotation angle by utilizing the calibrated internal and external parameters.

Firstly, matching the corner points according to the calibrated information. Let the corner coordinate of the picture taken at θ be X^c＝{(u_i,v_i)^T,i＝1,2,...,N_cAt this time, occlusion may exist in the picture, so that matching cannot be automatically completed, and the following algorithm may be used:

1. a set of points that may match is found.

For each corner point and target point pairCalculating a score S_i,jThe calculation method is as follows:

(1) obtaining (x)_j,y_j,0)^TCorresponding world coordinate

[\begin{matrix} x_{j}^{w} \\ y_{j}^{w} \\ z_{j}^{w} \end{matrix}] = R_{{wb}_{0}} [\begin{matrix} x_{j} \\ y_{j} \\ 0 \end{matrix}] + [\begin{matrix} t \\ 0 \\ 0 \end{matrix}] .

(2) Finding an image point (u)_i,v_i)^TIn thatPrimary image on planeThe obtaining method is as follows

Z_{w}^{j} = z_{w}^{j},

Wherein,

is (u)_i,v_i)^TThe corresponding normalized coordinates can be obtained iteratively.

(3) We define S_i,jThe following formula

S_{i, j} = \sqrt{{(X_{w}^{j})}^{2} + {(Y_{w}^{j})}^{2}} - \sqrt{{(x_{w}^{j})}^{2} + {(y_{w}^{j})}^{2}} . - - - (1.21)

(4) Selecting a threshold τ_mWhere d is the length of the checkerboard grid, then define the set of points that may match

<math> <mrow> <mi>Ω</mi> <mo>=</mo> <mo>{</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&LeftRightArrow;</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1,2</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <mi>n</mi> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>1,2</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msub> <mi>N</mi> <mi>c</mi> </msub> <mo>|</mo> <msub> <mi>S</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo><</mo> <msub> <mi>τ</mi> <mi>m</mi> </msub> <mo>}</mo> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1.22</mn> <mo>)</mo> </mrow> </mrow> </math>

2. Using the histogram, the angle is calculated.

For the elements in Ω, θ is calculated by the following equation_i,j

Then draw outA histogram of the distribution, the histogram having 360 squares, the number of elements in the kth square being (k-181) °, (k-179) °) # Ω. Our estimate of θ is

3. And selecting a correct matching pair by using the obtained angle.

i corresponds to the point M (j) with the smallest score in the correct angle range, i.e. the point M (j)

M (i) = \min_{j} S_{i, j},

Satisfies theta_i,j∈[θ-1°,θ+1°].

If theta is not present_i,j∈[θ-1°,θ+1°]Then, m (i) is made 0.

4. Non-linear optimization

The optimization problem is described as

Wherein,

the optimization problem can be solved by the lsqnolin function of the MATLAB itself.

The advantages and the effects are as follows:

the invention relates to a visual measurement method of a single-axis rotation angle based on a fixed camera and a single target, which is a method for measuring the single-axis rotation angle by using a monocular camera and a two-dimensional checkerboard. The invention adopts a two-dimensional target as a calibration feature, and can rotate the angle in a single axis through an optimization algorithm. The method has the advantages of simple installation process, high calibration precision, low cost and the like, and can solve the problem of shielding after calibration.

Drawings

FIG. 1: pinhole imaging model schematic

FIG. 2: schematic diagram of experimental equipment installation

FIG. 3 a: graph of mean value of measurement error as a function of noise

FIG. 3 b: measurement error variance versus noise plot

FIG. 4 a: angular point number graph detected when angle changes

FIG. 4 b: angle error map when angle changes

FIG. 5: two-dimensional checkerboard target picture used for real experiment

FIG. 6: flow chart of the invention

The symbols in the figures are as follows:

symbolic illustration in fig. 1: o denotes the camera optical center, O₁Representing the origin of the image coordinate system. X, Y denote coordinate axes of the image coordinate system and X, Y, z denote coordinate axes of the camera coordinate system. P (X)_w,Y_w,Z_w) Representing the coordinates of the three-dimensional point in the camera coordinate system and P representing the projection of the point P on the image.

Symbolic illustration in fig. 2: o is_w-X_wY_wZ_wAs a world coordinate system, Z_wThe axis is a rotation axis, and the origin of the target at the reference position is located at X_wOn the shaft; target coordinate system of target in reference positionThe coordinate system surrounds Z_wRotation theta_iAfter that, the target coordinate system becomesCoordinate system of camera O_c-X_cY_cZ_c。

Detailed Description

Referring to fig. 1-6, the invention provides a visual measurement method based on a single-axis rotation angle of a fixed camera and a single target, and simulation experiments and real experiment verification are carried out.

Simulation experiment

The simulation experiment is carried out on a computer with a main frequency of 3.07GHz and a memory of 4.00GB and MATLAB2014a in a Windows XP environment. The camera model used for simulation is a second-order distortion model, and the parameters are as follows:

α_x＝782.5109,α_y＝782.9155,u₀＝357.3909,v₀＝264.5240

k₁＝-0.4232,k₂＝0.2664,k₃＝-0.000555,k₄＝-0.0031,k₅＝0

the target was a 6X 8 checkerboard, with individual checkerboard lengths of 39.5mm by 39.5 mm. The simulation experiment mainly verifies the stability and high precision of the calibration algorithm. The simulation steps are as follows:

setting a simulation environment: and acquiring the coordinates of the image points according to the set internal parameters and the pose. Wherein the selected angle during rotation is (unit is degree)

θ_i＝2.25×i

The method comprises the following steps: using MATLAB to calibrate the tool box, obtaining the initial value of the internal parameter, and the result is as follows (calibrating the noise of the internal reference image to be 2 pixels):

α_x＝774.34147,α_y＝773.91064,u₀＝349.58372,v₀＝269.85450

k₁＝-0.43943,k₂＝0.33274,k₃＝-0.00222,k₄＝-0.00206,k₅＝0

step two: and calculating the initial value of the external parameter according to the calibration result of the common camera.

Step three: and binding and adjusting according to the initial values of the internal and external parameters.

Step four: and selecting an angle and measuring.

In order to verify the robustness of the algorithm, in step four, 100 angles are selected, measurement is performed, and the maximum value of the angle errors is taken as the measurement error. White gaussian noise with a mean value μ of 0 and a standard deviation σ of 0 □ 4 pixels was added to the image points, and the mean and variance of the measurement errors of 100 experiments were plotted, and the results are shown in fig. 3(a) and (b). The result shows that even if Gaussian white noise with the standard deviation of 4 pixels is added, the result with the angle error mean value smaller than 0.5 degrees can be obtained, and the result proves that the angle measurement method is robust to image noise. To verify that our method works for occlusion cases, we simulated the case when the target was rotated at (-180 °,180 °). Fig. 4(a) and (b) show the simulation results, where (a) indicates the number of detected corner points, and (b) indicates the measured angle error. Experimental results show that the method can effectively solve the problem of shielding.

Real experiment

To further verify the feasibility of the present invention, we performed real experimental verification. The target we used was a printed checkerboard (as shown in figure 5), with a 6 x 9 gauge, each grid being 29.5mm in length; the adopted camera is a Basler camera, the resolution is 658 x 492 pixels, and the focal length is about 800 pixels. The angle verification adopts a photoelectric tracking system, and the angle measurement precision of the system is 0.01 degree. The results are shown in table 1, with errors indicating that the maximum value of the angular error is less than 0.15 degrees under the experimental conditions.

TABLE 1 true Experimental rotation Angle estimation results

	1	2	3	4	5	6	7	8	9	10
											Rotary table value	-54.00	-36.00	-27.00	-10.80	-1.8	5.40	16.20	28.80	36.00	54.00
Truth value	-53.985	-35.978	-26.966	-10.749	-0.036	5.477	0.070	28.894	36.135	54.147
											Error of the measurement	-0.015	-0.022	-0.034	-0.051	-1.764	0.077	16.270	0.094	0.135	0.147

Claims

1. A vision measurement method based on a single-axis rotation angle of a fixed camera and a single target is characterized in that: the method comprises the following specific steps:

the method comprises the following steps: calibrating internal and external parameters of a camera;

in the specific implementation process, a two-dimensional checkerboard is used as a target; the adopted picture sequence comprises a picture at a reference position and pictures when the m rotating components rotate to different angles; suppose that the m angles are θ₁,θ₂,...,θ_mNext, the checkerboard corner points of MATLAB2014 are utilized to extractTaking a function detectCheckerboardPoints for corner extraction, and combining an MATLAB calibration tool for internal parameter calibration to obtain an internal parameter p_inAnd external parameters of m +1 poses; the calibration tool box is downloaded on http:// www.vision.caltech.edu/bouguetj/caliib _ doc/; these extrinsic parameters are expressed as

[R_{{cb}_{i}}, T_{{cb}_{i}}], i = 0,1,2, . . ., m .

Expressed at a rotation angle of theta_iWhen the i is 0, the rotation matrix and the translation vector between the camera and the target are represented as the reference position;

step two: initial value of external parameter related to pivoting motion

(1) Determining the angle of rotation theta from (1.7) based on the external parameters determined in step one_iThe rotational and translational vectors of the target relative to the reference

[R_{b_{i} b_{0}}, T_{b_{i} b_{0}}], i = 1,2, . . ., m .

R_{b_{i} b_{0}} = R_{{cb}_{i}}^{T} R_{{cb}_{0}}, T_{b_{i} b_{0}} = R_{{cb}_{i}}^{T} (T_{{cb}_{0}} - T_{{cb}_{i}}) . - - - (0.1)

(2) And obtainingt and theta_iInitial value:

Then, obtain

Wherein,S_i,jelements representing the ith row and j columns of the matrix S, r₁Obtained from the following formula

The above type passing pairObtained by SVD decomposition if

Then x is the third column of V, then, get

θ_i＝2arctan((A^TA)^-1(A^TB)), (0.6)

Wherein

A = vec (r_{2} r_{3}^{T} + r_{3} r_{2}^{T}),

Vec (-) here represents the straightening of the matrix, so,

R_{{wb}_{0}} = R_{b_{0} w}^{T} - - - (0.7)

(3) according to the obtainedt and theta_iSolving the initial value, and solving R according to (1.16)_cw,T_cwValue of (A)

R_{cw} = R_{{cb}_{0}} R_{{wb}_{0}}^{T}, T_{cw} = {(t, 0,0)}^{T} - R_{cw} T_{{wb}_{0}} . - - - (0.9)

Step three: intra-pair parameter p_inAnd an extrinsic parameter p_outCarrying out nonlinear optimization; where p is_in∈□⁹,p_out∈□¹⁰Is shown below

Wherein,are respectively asR_cwThe corresponding Rodrigues vector; the conversion relationship between a Rodrigues vector v and its corresponding rotation matrix R is as follows

Wherein, I₃An identity matrix of order 3;

the coordinates of the individual image points are marked on the target 1 at a rotation angle θ_iThen, the corner point coordinate corresponding to the jth target feature point is recorded as

(u_{j}^{i}, v_{j}^{i}), i = 0,1, . . ., m, j = 1,2, . . ., n,

Points on the target are noted

X_b＝{(x_j,y_j,0)^T,j＝1,2,...,n}.

Because errors exist in the image extraction process, the reprojection error sum of the images is selected as an optimization target, and the target point and the p_outObtaining the reprojection pointsThe process of (2) is as follows:

<1>angle of rotation theta_iTime, target point (x)_j,y_j,0)^TConversion to points in world coordinate system

<2>And calculating the corresponding reprojection point

The reprojection error is defined as

The optimization problem is described as

The optimization problem is solved by utilizing a sparse Levenberg-Marquardt algorithm, and an optimization tool box of the sparse LM algorithm is downloaded in http:// users. ics.forth.gr/. Loursakis/spark LM; with the above optimization, obtaint,R_cw,T_cwThe information of (1), which is used for the rotation angle calculation of the step four;

step four: measuring the rotation angle by utilizing the calibrated internal and external parameters;

firstly, matching angle points according to calibrated information; let the corner coordinate of the picture taken at θ be X^c＝{(u_i,v_i)^T,i＝1,2,...,N_cAt this time, occlusion may exist in the picture, so matching cannot be automatically completed according to the following algorithm:

scattering and searching a point set which is possibly matched;

(1) obtaining (x)_j,y_j,0)^TCorresponding world coordinate

[\begin{matrix} x_{j}^{w} \\ y_{j}^{w} \\ z_{j}^{w} \end{matrix}] = R_{{wb}_{0}} [\begin{matrix} x_{j} \\ y_{j} \\ 0 \end{matrix}] + [\begin{matrix} t \\ 0 \\ 0 \end{matrix}] .

Z_{w}^{j} = z_{w}^{j},

Wherein,

is (u)_i，v_i)^TThe corresponding normalized coordinates are obtained by iteration;

(3) definition of S_i,jThe following formula

S_{i, j} = \sqrt{{(X_{w}^{j})}^{2} + {(Y_{w}^{j})}^{2}} - \sqrt{{(x_{w}^{j})}^{2} + {(y_{w}^{j})}^{2}} . - - - (0.15)

<math> <mrow> <mi>Ω</mi> <mo>=</mo> <mo>{</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>v</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&LeftRightArrow;</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1,2</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <mi>n</mi> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>1,2</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <msub> <mi>N</mi> <mi>c</mi> </msub> <mo>|</mo> <msub> <mi>S</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo><</mo> <msub> <mi>τ</mi> <mi>m</mi> </msub> <mo>}</mo> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>0.16</mn> <mo>)</mo> </mrow> </mrow> </math>

Utilizing the histogram to calculate an angle;

for the elements in Ω, θ is calculated by the following equation_i,j

Then draw outA histogram of distribution having 360 squares, the number of elements in the kth square being (k-181) °, (k-179) ° n Ω; an estimate of theta of

Selecting a correct matching pair by using the obtained angle;

M (i) = \min_{j} S_{i, j},

Satisfies theta_i,j∈[θ-1°,θ+1°].

If theta is not present_i,j∈[θ-1°,θ+1°]Then, let m (i) be 0;

good and non-linear optimization

The optimization problem is described as

Wherein,

the optimization problem is solved by the lsqnolin function of the MATLAB itself.