Automatic Calibration between Multi-Lines LiDAR and Visible Light Camera Based on Edge Refinement and Virtual Mask Matching

Abstract

To assist in the implementation of a fine 3D terrain reconstruction of the scene in remote sensing applications, an automatic joint calibration method between light detection and ranging (LiDAR) and visible light camera based on edge points refinement and virtual mask matching is proposed in this paper. The proposed method is used to solve the problem of inaccurate edge estimation of LiDAR with different horizontal angle resolutions and low calibration efficiency. First, we design a novel calibration target, adding four hollow rectangles for fully automatic locating of the calibration target and increasing the number of corner points. Second, an edge refinement strategy based on background point clouds is proposed to estimate the target edge more accurately. Third, a two-step method of automatically matching between the calibration target in 3D point clouds and the 2D image is proposed. Through this method, i.e., locating firstly and then fine processing, corner points can be automatically obtained, which can greatly reduce the manual operation. Finally, a joint optimization equation is established to optimize the camera’s intrinsic and extrinsic parameters of LiDAR and camera. According to our experiments, we prove the accuracy and robustness of the proposed method through projection and data consistency verifications. The accuracy can be improved by at least 15.0% when testing on the comparable traditional methods. The final results verify that our method is applicable to LiDAR with large horizontal angle resolutions.

1. Introduction

At present, traditional large-scale digital map technology has been well developed and applied, such as Google Maps, Baidu Maps, etc. However, the work related to fine terrain reconstruction in a small scale is relatively limited. In order to obtain the fine 3D topographic map in a small scale, the simplest way is to use unmanned aerial vehicles (UAVs) to reconstruct terrain using the tilt photography technology. Tilt photography technology uses pure 2D image analysis and modeling to restore the 3D structure from 2D information, and its calculation accuracy is directly affected by image quality and imaging environment. Recently, more and more research has been conducted on the fine 3D terrain reconstruction technology using UAVs equipped with light detection and ranging (LiDAR) and visible light cameras. The visible light camera can obtain high-resolution color information, but it is particularly vulnerable to the influence of external weather, illumination, and other factors; it also lacks the 3D information of the target. LiDAR can quickly obtain the 3D information of space objects, but it cannot get the texture, color, and other information of objects. Therefore, the LiDAR and visible light camera can achieve excellent complementary effects in performance. This will greatly improve the performance of existing UAVs 3D terrain reconstruction or low-altitude remote sensing works [1,2,3,4,5,6,7]. However, since the data obtained by the two sensors are based on their respective coordinate systems, and data fusion requires that the data collected by two sensors should be expressed in LiDAR or camera coordinate system. It is necessary to determine the transformation matrix between two coordinate systems through joint calibration, that is, the extrinsic parameters between LiDAR and the camera together with the intrinsic parameters of the camera.

Currently, according to the characteristics of different methods, we divide the calibration between the LiDAR and visible light camera into a target-based method and a target-less method. For details and representative works, refer to Figure 1. The target-based method is to find the 2D feature points in an image coordinate system and the 3D feature points in a LiDAR coordinate system with the help of a standard target to establish geometric constraints and solve the extrinsic parameters of sensors through perspective-n-point (PNP) or nonlinear optimization. According to the shape of calibration target, it can be divided into 1D objects, such as line-feature-based target [8,9,10], 2D objects, such as chessboard [11,12,13] or circular holes [14], and 3D objects [15,16], such as spherical targets [17]. Refs. [18,19] detected the corner points in images and point clouds by the intersection of edge fitting lines, respectively, and solved calibration parameters by a linear equation. The calibration target they used was a monochrome board. This is mainly due to the fact that two different colors have different effects on LiDAR ranging, which will affect the accuracy of plane fitting. Refs. [20,21,22] introduced nonlinear optimization to improve the calibration accuracy. Ref. [23] adopted the method of space joint calibration combining coarse measurement and fine adjustment. The innovation of [24] was to use a known plane to estimate the 3D corners. Refs. [24,25] used a LiDARTag and intensity information to locate target, which was similar to the method in [26].

The target-less method can be divided into mutual-information-based approach [27,28], motion-based method [29,30], and learning-based technique [31,32,33,34,35]. The mutual-information-based approach can make full use of environmental information to complete online calibration. However, the disadvantage is that it needs to be carried out in the natural scene, and if the scene characteristics cannot meet the expectations, it will lead to considerable deviation of results. The motion-based method takes the calibration of LiDAR and camera as a hand–eye calibration problem, which can recover the right transformation by a series of rigid transformations. This method needs accurate motion information, so the static system is not applicable. Many learning networks use the convolutional neural network (CNN) to solve this task. These networks can calibrate the LiDAR–camera system without using calibration targets, matching information, and motion information, but the accuracy is dependent on the size of the training data set and the structure of the CNN. In the application, it requires high texture information of the environment and high computational power, and often needs to use a graphic processing unit to accelerate data processing.

The target-based method is still the mainstream because of its stability and low computational performance requirements. Therefore, sensors will undergo stable and high-precision target-based calibration methods to calibrate the extrinsic parameters before leaving the factory. The existing methods mainly have the following problems. First, the adaptability of LiDAR with different horizontal angle resolutions (θ) of calibration methods is not considered. For example, the method based on edge estimation [18,19] had a great difference in the effectiveness of LiDAR with θ = 0.1° and large θ = 0.4°. The second is that the camera intrinsic parameters, such as distortion coefficients, are not considered in the optimization process. The two-stage method, first calibrating the intrinsic parameters and then solving the extrinsic parameters, brings a lot of inconveniences to large-scale commercial applications in terms of efficiency. Third, for the locating of the calibration target, manual filtering is not conducive to the improvement of the automation level. The method of using intensity information to locate special marks will have different effects on LiDAR ranging due to two different color objects.

To solve the adaptability of edge fitting under different θs, we propose an approximation edge fitting technology based on background point clouds, which can refine the edge points under different θs by using a known geometric size of the calibration target. According to the geometric shape of the designed calibration board, a set of corner points extraction schemes of full-automatic locating and fine processing are designed, which can greatly reduce the manual operation and do not depend on the intensity information of LiDAR and special marks. Finally, camera intrinsic parameters are taken into account, and their initial values are sent to the optimization function, which avoids the inefficiency caused by the two-stage method when calibrating a large amount of equipment. The main contributions of this paper are as follows: We design a novel calibration target and propose a joint automatic calibration method based on edge refinement and virtual mask matching. An improved edge refinement scheme is introduced to refine the edge points. The maximum error of the edge estimation of the calibration board in this method does not exceed 5.0 mm, which greatly improves the accuracy of the corner points. The high-precision corner detection method makes the final calibration accuracy better than state-of-the-art techniques, especially in scenes with sparse point clouds, such as θ is 0.4°. An automatic location method, locating firstly and then fining processing, is proposed to get the corner points, which do not need any parameters set by users in the calibration process. In the experiments, we take multi-lines LiDAR as an example, such as 64 lines, and verify the accuracy and robustness of proposed method.

In the following sections, Section 2 is a brief introduction of the proposed method, Section 3 is the concrete introduction, Section 4 is the experiment and discussion, and Section 5 is the conclusion.

2. Proposed Calibration System and Computational Flow Chart

Figure 2a,b show the design of the calibration board. We need to place the calibration board about 30.0 cm in front of a flat wall to better obtain the points of the background wall for edge refinement, as shown in Figure 2c. In Section 3.1, we will use Euclidean distance to cluster all objects. In order to prevent the calibration board and the flat wall from being clustered into a category, the calibration board needs to keep a certain distance from the wall, e.g., 30.0 cm is considered in this paper. Actually, other suitable distances can also be used in practical applications. Holes dug in the calibration board are convenient for us to use geometric information for locating. Figure 2d is a diagrammatic sketch of this method, which involves the locating and matching of 3D feature points, the locating and refinement of 2D feature points, the solution of optimization equations, and the 3D–2D projection results. If we unify the data to the LiDAR coordinate system, we can get the point clouds with color, or unify the data to the camera coordinate system, and we can also get an image with sparse depth information. The symbols and function we will use in the following sections are given in

Table 1. Symbols and functions definitions.

Symbols and Functions	Description
$p_i^L\left(x_i^L,y_i^L,z_i^L\right)$	a 3D point in LiDAR coordinate system
$p_i^C\left(x_i^C,y_i^C,z_i^C\right)$	a 3D point in camera coordinate system
$q_i\left(u_i,v_i\right)$	a 2D point in image coordinate system
$P/Q$	a 3D/2D points set
$P{C}_i$	a point cloud belonging to a same class of objects
$V_i$	a 3D space vector
scale	a scale factor about $V_i$
ss	the step size of scale
th	the expanded threshold of edges
$rin{g}_i$	a scanning line of LiDAR
θ/θs	θ is the horizontal angle resolutions, θs is the plural form of θ
$R\left(r_x,r_y,r_z\right)$	a rotation matrix around the X, Y, and Z axes
$T\left(t_x,t_y,t_z\right)$	a translation matrix in the direction of X, Y, and Z axes
$u_x/u_y$	the coordinates of the image principal point
$f_x/ f_y$	the scale factors in image $x$ and $y$ axes
$k_1, k_2$	the distortion coefficient of the image
$Box$	a bounding box in 2D plane
$Score\left(P{C}_i\right)$	a matching score between $P{C}_i$ and the preset point cloud
$N$($P{C}_i$)	the number of points in $P{C}_i$
$S$($Box\left(P{C}_i\right)$)	the area of $Box$ enclosing $P{C}_i$
c	the number of cells for image segmentation
ε	the gray threshold of image
$E\left(q,q^{\prime }\right)$	the error function of $q$ and $q^{\prime }$
$HuberLoss\left(E\left(q,q^{\prime }\right)\right)$	the loss function of $E\left(q,q^{\prime }\right)$

.

3. Proposed Method

3.1. Automatic Locating Calibration Board

In this section, we describe how to automatically locate the point clouds of the calibration board ($P{C}_{board}$) and the background wall ($P{C}_{wall}$). First, to reduce the influence of ground points on 3D clustering, we use $patchwork$ [36] to remove the ground. This is a robust algorithm for removing the ground, which can reduce the $ring$ filtering parameters. The effect is illustrated in Figure 3a. We can cluster all objects according to Euclidean distance (Figure 3b), and each object after clustering is called$P{C}_i$. The normal vector $V_{P{C}_i}$ of $P{C}_i$ is obtained by plane fitting, shown in Figure 3c. Performing rotation transformation according to $V_{P{C}_i}$ and $V_X\left(1,0,0\right)$ to transform the $P{C}_i$ to the YOZ plane, which is called $PC{\prime }_i$. We take the geometric centroid of $PC{\prime }_i$ as the center and generate a box with the same size as our calibration board in the YOZ plane, as shown in the Figure 3d,e. We calculate the matching $Score$ between each $PC{\prime }_i$ and the preset box, and redefine $PC{\prime }_i$ as the set of $P{C}_{in}$ and $P{C}_{out}$:

$$PC{\prime }_i=P{C}_{in}{{\displaystyle \cup }}^{ }P{C}_{out}$$

(1)

where $P{C}_{in}$ is the set of red points in Figure 3d, $P{C}_{out}$ is the set of black and blue points, and the $Score\left(PC{\prime }_i\right) \mathrm{is}$:

$$Score\left(P{C^{\prime }}_i\right)=ScoreA\times ScoreB\times ScoreC$$

(2)

where $ScoreA=N(P{C}_{in})/N(P{C}_{out}$); $ScoreB$ is the ratio of the short side ($w$) to the long side ($h$) of the $Box\left(P{C}_{in}\right)$, $ScoreB=w/h$; $ScoreC$ is the ratio of $S\left(Box(P{C}_{in}\right))$ to the actual calibration board square, $ScoreC$=$S\left(Box(P{C}_{in}\right))$/$S\left(board\right)$. Since our calibration board area is just 1.0 m², $ScoreC$ = $S\left(Box(P{C}_{in}\right))$.

Let us take Score as the evaluation function to screen $PC{\prime }_i$. Three parts of $Score$ can filter out the objects with a large or small number of points to ensure that the cluster where the calibration board is located can be filtered out with a high score. Figure 3e is the one with the highest Score among many point clouds. In Section 4.2, we will analyze the impact of ScoreA, ScoreB, and ScoreC on the automatic locating calibration board. It can be found that the connecting rod of the calibration board is also selected, shown in Figure 3e. The results of Figure 3c,d are ideal. In fact, due to the Euclidean clustering, the connecting rod of the calibration board will also be clustered into one category. We need to discard it for it should not participate in the calculation. Finally, we can delete the ring with a small number of point clouds by traversing each ring, as shown in the Figure 3f. The obtained red point cloud is $P{C}_{board}$. We select a 1.0 m × 1.0 m × 1.0 m box in the center of $P{C}_{board}$ to ensure that the $P{C}_{wall}$ (blue point clouds and green point clouds in Figure 3g) is framed together. Therefore, the calibration board is required to be placed within 1.0 m in front of the background wall. So far, we have acquired $P{C}_{board}$ and $P{C}_{wall}$ through automatic locating, and we can further refine the edge points by the method described in next section.

3.2. Approximation Edge Fitting

This section introduces an approximation edge fitting method, which can better fit the problems of inaccurate edge estimation caused by θ. As shown in Figure 4a, due to θ of LiDAR, two adjacent laser beams, one on the target and the other outside the target, cannot really scan the edge of the object. It can be seen from Figure 4b,c that although the laser spots of LiDAR have randomness within a certain range, overall, the blue laser spots with θ = 0.4° are far from the real edge. Figure 4d–f show the sparsity of different θs The distance between the measured value and the real edge is closely related θ and the test distance, and the larger the θ, the worse the edge estimation of the target. Except the θ, vertical angle resolution, related to the number of rings of the LiDAR, may also be another factor. From Figure 4d–f, we can find that under the same number of lines of the LiDAR, the measurement results under different θs have different density levels of point clouds. Therefore, in the following sections, we take the LiDAR with 64 lines as an example (details are listed in Section 4.1) to research the calibration of different θs, while ensuring the same vertical angle resolution of the LiDAR, paying attention to the effect improvement in the case of different θs.

As shown in the Figure 5a, for the partial schematic diagram of the calibration board scanned by LiDAR, we use $rin{g}_i$ as an example, where the green background edge point is $p_i^b\left(x_i^b,y_i^b,z_i^b\right)$, the red foreground edge point on the board is $p_i^f\left(x_i^f,y_i^f,z_i^f\right)$, the direction vector (blue dotted line) between these two points is $V_i\left({\stackrel{\rightharpoonup }{x}}_i,{\stackrel{\rightharpoonup }{y}}_i,{\stackrel{\rightharpoonup }{z}}_i\right)$, where ${\stackrel{\rightharpoonup }{x}}_i=x_i^b-x_i^f , {\stackrel{\rightharpoonup }{y}}_i=y_i^b-y_i^f, {\stackrel{\rightharpoonup }{z}}_i=z_i^b-z_i^f$. The set of background edge points is $P^b,p_i^b\in P^b$; the set of edge points of the foreground calibration object is $P^f,p_i^f\in P^f$.

Due to θ, the red line $l_f$ fitted by $P^f$ is still a certain distance from the actual black edge line, so the quadrilateral fitted with these edge points is often smaller than the actual size, and with the increase of test distance or θ, the difference between the quadrilateral and the real size is larger (the verification is shown in Section 4.2). It is necessary to expand $P^f$ to approach the actual edge. First, we should get the points on the calibration board ($P{C}_{board}$) and background wall ($P{C}_{wall}$), which have been achieved in Section 3.1. Then, the edge points on the background, i.e., yellow points in Figure 5c, should be projected to the calibration board plane according to the direction of the laser beam of LiDAR. We match the foreground point (red point in Figure 5b) and background point (green point in Figure 5b) of each scanning $rin{g}_i$, which can be easily realized by the sorting algorithm. At this time, we can obtain the background edge point $P_i^b$ and calibration board edge point $P_i^f$. We choose $V_i$, the direction of the blue arrow in Figure 5a,c as the outward expansion direction. The outward length of the expansion is the distance between two points multiplied by a $scale$ coefficient, $scale\in \left(0 ~ 0.5\right)$, and the blue middle edge point $p_i^m\left(x_i^m,y_i^m,z_i^m\right)$ in Figure 5c can be obtained:

$$p_i^m=p_i^f+scale\ast V_i$$

(3)

where $p_i^m$ must be located between $p_i^f$ and $p_i^b$. By traversing all the $ring$s, we can obtain the remaining edge point set $P^m$ after approximation. We use $P^m$ as the new edge points to fit the edge line and obtain the blue edge line $l_1$, shown in Figure 5b. Similarly, the edge line $l_2, l_3$, and $l_4$ can also be computed.

As for the selection of $scale$, we use the actual size of the calibration board (Figure 2b) for reference. Taking $scale=$ 0.1 as the initial value. The newly expanded edge points are used to fit a rectangle by random sample consensus (RANSAC). The termination condition of RANSAC fitting is the difference between the calculated rectangle size and the actual physical size. If this difference is less than the th, the outward expansion can be stopped, and the edge of the current outward expansion should be saved. If the difference is larger than th, it indicates that there is still a gap between the measured result and the actual size, and it needs to be further expanded. That is, scale needs to continue to increase with a step size (ss):

$$scale=scale+ss$$

(4)

where scale is not recommended to exceed 0.5, because when scale is greater than 0.5, $p_i^m$ is closer to the green point in Figure 5. If scale is less than 0.5, $p_i^m$ is closer to the red point in Figure 5. We prefer the edge point to be the point on the calibration board rather than outside. The $scale$ is expanded to 0.4 in ss of 0.01 in the later experiment.

In the subsequent experiments of this paper, th = 0.005 m is selected as the experimental data. The smaller the $th$, the higher the estimation of the edge point accuracy, but the longer the detection time. By repeating the above process, a rectangular edge point cloud with an error of no more than 5 mm from the actual physical size can be measured (we will verify this in Section 4.2).

3.3. 3D Virtual Mask Matching

This section introduces a scheme for obtaining 3D feature points of the calibration board by using virtual mask matching. The edge point clouds,$P{C}_{edge}$and $P{C}_{board}$, are combined into a new point clouds $(P{C}_{calib}),$ $P{C}_{calib}=P{C}_{edge}$∪$P{C}_{board}.$ We generate a virtual mask $P{C}_{mask}$ with $N$($P{C}_{mask}$) random points in YOZ plane with the same size as our board, shown in Figure 6, set $N$($P{C}_{mask}$) = 10⁴. We use the iterative closest point (ICP) method to match the $P{C}_{mask}$with the $P{C}_{calib}$, and the matching result is a rigid transformation matrix$H$. The feature points $P\left(p_1^L,p_2^L,\dots p_{20}^L\right)$ we need can be calculated by $P\prime *H$, where $P\prime \left(p_1^{\prime L},p{\prime }_2^L,\dots p{\prime }_{20}^L\right)$ are feature points set on the $P{C}_{mask}$.

3.4. 2D Corner Points Detected in Image

We design an automatic matching method according to the characteristics of our calibration board, which can automatically find corners in images and finely process corners. It is mainly divided into two steps: corner points locating and corner points refinement. First, we use a region growing algorithm (RGA) to locate 2D feature points roughly. We divide the input image into c × c cells, and use the center point of each cell as the seed point of RGA, as shown in Figure 7a. In the next experiment, we choose $c$ = 10, that is, the original image is divided into 10 equal parts horizontally and vertically, and 100 seed points are obtained. Since our calibration board is pure white, we only grow the seed points whose gray value is greater than ε, and set ε to 127, half of the maximum grayscale value, which can quickly filter out invalid seed points. Then, rectangle detection is carried out for the RGA results of each seed point, shown in Figure 7b. It is obvious that multiple rectangles are detected at the same location, shown in Figure 7c, which is not conducive to the selection of corner points. Finally, we perform the k-mean clustering method on all rectangles. The clustering index is the distance between the center points of each rectangular box. The k-value of clustering is five, which exactly corresponds to five rectangles to be detected. The clustering result is five clusters of rectangular boxes, and each cluster has at least one rectangle. The mean value of corner points of each cluster of rectangular boxes is taken as the corner points for roughly locating position, as red points shown in Figure 7d.

After obtaining the corner points of roughly locating, it can be clearly found that the locating position is not accurate (red points in Figure 7d–f). This is because the rectangle detection has poor robustness and accuracy. It is necessary to further refine the corner points. Taking $q_1$ and $q_2$, the green points in Figure 7g as an example, we connect $q_1$ and $q_2$ to obtain a straight line $l_{q_1}^{q_2}$, count all the points close to $l_{q_1}^{q_2}$, which are circled by a yellow ellipse in Figure 7g. We use these points to refit the line and get the refined$l{\prime }_{q_1}^{q_2}$. We solve the connecting lines between the other points in the same way and estimate the refined corner points by the intersection of straight lines. As described in the green dots, it can be found that the green points after the fine processing are more accurate than the red points detected before. In Figure 7g, for convenience of representation, we visualize the result in pixels. However, in fact, the coordinates of these corner points are sub-pixel. Finally, according to the order of feature points in Figure 7h, each refined corner point is sorted to facilitate pairing with 3D feature points.

3.5. Optimization Equation Modelling

In this section, an optimization equation is established to find the best intrinsic parameters of the camera and solve the extrinsic parameters of LiDAR and camera. The optimization process is intended to find the best LiDAR–camera extrinsic parameters and camera intrinsic parameters to minimize the re-projection error between the 3D–2D feature points. Through the methods of Section 3.1, Section 3.2, Section 3.3 and Section 3.4, we can obtain 3D and 2D feature points after refinement processing. We use $P_i\left(p_{i,1}^L,p_{i,2}^L,\dots p_{i,20}^L\right)$ and $Q_i\left(q_{i,1},q_{i,2},\dots q_{i,20}\right)$ to represent 3D/2D feature points set measured at the i th position, $i\in \left(1,n\right)$, n is the number of test positions. The projection process from 3D points to 2D points needs to go through a rigid transformation model, pinhole imaging model, and distortion model. The detailed formulas are as follows:

$$\left[\begin{array}{c}{x}^C\\ \begin{array}{c}{y}^C\\ z^C\end{array}\\ 1\end{array}\right]=\left[R^{*}|T^{*}\right]\left[\begin{array}{c}{x}^L\\ \begin{array}{c}{y}^L\\ z^L\end{array}\\ 1\end{array}\right]$$

(5)

$$\left[\begin{array}{c}{q}^u\\ q^v\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_x& 0& u_x\\ 0& f_y& v_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x{\prime }^C\\ y{\prime }^C\\ 1\end{array}\right]$$

(6)

where $x{\prime }^C=x^C/z^C,y{\prime }^C=y^C/z^C, R^{*}$ represents the best rotation matrix, $T^{*}$ represents the best translation matrix, and $K^{*}=\left\{f_x,f_y,u_x,v_y\right\}$ represents the best camera intrinsic, $q^u$ and $q^v$ are the coordinates of a 2D point in the image. Now we consider the distortion model of the camera, and we use $D=\left\{k_1,k_2\right\}$ to represent it. The final imaging of the camera is as follows:

$$r^2=x{\prime }^C\ast x{\prime }^C+y{\prime }^C\ast y{\prime }^C$$

(7)

$$\begin{array}{c}{x}^{″C}=x{\prime }^C\left(1+k_1\ast r^2+k_2\ast r^4\right)\\ y^{″C}=y{\prime }^C\left(1+k_1\ast r^2+k_2\ast r^4\right)\end{array}$$

(8)

$$\left[\begin{array}{c}{q}^u\\ q^v\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_x& 0& u_x\\ 0& f_y& v_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}^{″C}\\ y^{″C}\\ 1\end{array}\right]$$

(9)

According to Equations (5)–(9), we establish the following optimization equations to solve $R^{\ast },T^{\ast },K^{\ast }$, and the best distortion parameter $D^{\ast }:$

$$\left[R^{\ast },T^{\ast },K^{\ast },D^{\ast }\right]=\underset{R,T,K,D}{\arg \min }{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^{20}HuberLoss\{E_{i,j}[q_{i,j}, KD\left(R{p}_{i,j}^L+T\right)]\}$$

(10)

where the form of $E_{i,j}$(·) is:

$$E_{i,j}\left(q,q^{\prime }\right)=\sqrt{{\left(q^u-q{\prime }^u\right)}^2+{\left(q^v-q{\prime }^v\right)}^2}$$

(11)

and $HuberLoss(·)$ is:

$$HuberLoss\left(\delta \right)=\{\begin{array}{cc}\delta & \delta \le 1\\ 2\sqrt{\delta }-1& \delta >1\end{array}$$

(12)

Therefore, our optimization problem can be transformed into the optimization of f(Q,P). The initial value of the camera’s intrinsic parameters is given in advance, which can come from the coarse calibration result of the camera, or the initial value given by the camera manufacturer. We use ceres optimizer (http://ceres-solver.org/, accessed on 1 December 2022) to optimize the intrinsic and extrinsic parameters. In the optimization, we set the maximum iteration as 500 and stopping tolerance as 10⁻¹⁰.

$$f\left(Q,P\right)={\displaystyle \sum }_{\begin{array}{c}{P}_i\in P\\ Q_i\in Q\\ i=1\end{array}}^n{\displaystyle \sum }_{\begin{array}{c}{p}_{i,j}^L\in P_i\\ q_{i,j}\in Q_i\\ j=1\end{array}}^{20}HuberLoss\{E_{i,j}[q_{i,j}, KD\left(R{p}_{i,j}^L+T\right)]\}$$

(13)

4. Experiments and Discussions

4.1. Proposed Experiment System

We setup the system with a LeiShen LiDAR with 64 lines and DaHeng camera with resolution of 1920 × 1200, shown in Figure 8a, and the field of view (FOV) of these two sensors is shown in Figure 8b. The size of the calibration target has been shown in Figure 2b. The maximum value of LiDAR’s detection range is 100.0 m. The wavelength is 905 nm. The rotation speeds of our LiDAR can be set to 300.0 rpm (i.e., 5.0 Hz), 600.0 rpm (i.e., 10.0 Hz), and 1200.0 rpm (i.e., 20.0 Hz), and the corresponding θ are 0.1°, 0.2°, and 0.4°, respectively. We simulate the LiDAR with different θs by switching the rotation speed of the LiDAR motor. We mainly analyze the influence of θ on calibration, so other variables are fixed. In the controlled experiments of different methods in this section, only the θ is different.

4.2. Experiments

• Experiment 1: 3D feature points detected

This section verifies the accuracy of approximation edge points. We take the measurement results of the size of the calibration board as the evaluation index and analyze the evaluation index by the original edge points and the approximated edge points.

Table 2. Calibration board size measurement results.

Distance(m)	θ(°)	BB		Scale = 0.1		Scale = 0.2		Scale = 0.4
		w (mm)	h (mm)	w (mm)	h (mm)	w (mm)	h (mm)	w (mm)	h (mm)
7.0	0.1	1012.34	996.22	984.11	981.82	997.90	991.14	998.51	997.95
	0.2	994.32	1004.66	979.91	987.12	993.41	997.94	997.20	996.22
	0.4	1003.72	993.731	974.41	988.20	994.90	997.01	997.81	999.17
10.0	0.1	1003.50	967.51	981.64	970.10	992.73	995.74	998.20	999.46
	0.2	984.19	1004.14	970.04	983.17	990.73	997.76	997.21	998.03
	0.4	992.16	1009.49	981.64	970.10	992.73	995.74	994.20	996.76
12.0	0.1	977.75	979.99	960.40	971.44	982.07	990.14	998.00	994.41
	0.2	978.03	979.76	954.70	988.00	991.43	996.01	996.24	996.90
	0.4	975.39	969.98	954.55	965.54	981.14	989.33	997.29	996.71

The w and h, respectively, represent the lengths of two sides of a rectangle fitted with the currently refined edge points; BB represents the 2D bounding box of this point clouds.

displays the measurement results under three θs when the test distance is 7.0, 10.0, and 12.0 m. Without any processing, due to the influence of θ, the length of the quadrilateral was lower than the actual value (1000.0 mm × 1000.0 mm), which also means that the selection of edge points was not accurate. The size of the bounding box (BB) will be much larger than the actual size or much smaller, which indicates the instability of initial edge points. After the approximation processing, when scale = 0.4, the measurement error reached within 5.0 mm. Compared with the results of BB, our method is more accurate for the estimation of edge points. The result can be seen that our approximation processing can make the measurement result reach the true size. In Figure 9, we can find that the edge estimation results at different θs are preferable to the initial values. In the following experiments, we will apply this approximation refining strategy to [18,19,22]. The final calibration results are also greatly improved.

• Experiment 2: Automatic locating results

The experiment in this section is to verify the effect of automatic locating of the calibration board. After clustering in one frame data, nine objects, shown in Figure 10 and

Table 3. Matching score of different objects.

θ(°)	Object	N	ScoreA	ScoreB	ScoreC	Score	IsTrue
0.1	Figure 10a	34,884	0.031361	0.997999	0.996798	0.031198	False
	Figure 10b	4180	0.136842	0.992114	0.987158	0.134019	False
	Figure 10c	1116	0.870968	0.994576	0.982246	0.850865	True
	Figure 10d	312	0.282051	0.134227	0.250382	0.009479	False
	Figure 10e	283	0.745583	0.461738	0.302227	0.104046	False
	Figure 10f	263	0.365019	0.279073	0.349486	0.035601	False
	Figure 10g	243	0.786008	0.400130	0.255705	0.080421	False
	Figure 10h	202	0.876238	0.338324	0.215452	0.063871	False
	Figure 10i	201	0.671642	0.323524	0.214348	0.046576	False
0.2	No. 1	9514	0.04225	0.99512	0.98902	0.04158	False
	No. 2	7917	0.04724	0.99461	0.99401	0.04670	False
	No. 3	2004	0.12974	0.99707	0.98869	0.12789	False
	No. 4	555	0.81441	0.98390	0.94908	0.76051	True
	No. 5	156	0.30769	0.17733	0.18889	0.01030	False
	No. 6	142	0.76056	0.46460	0.30043	0.10616	False
	No. 7	132	0.36363	0.29602	0.31777	0.03420	False
	No. 8	106	0.94340	0.39433	0.36821	0.13698	False
	No. 9	101	0.66336	0.34315	0.21464	0.04886	False
0.4	No. 10	3887	0.02727	0.96699	0.96381	0.02541	False
	No. 11	2276	0.06678	0.48287	0.80706	0.02602	False
	No. 12	1689	0.13025	0.99809	0.99416	0.12924	False
	No. 13	529	0.24952	0.99838	0.99115	0.24692	False
	No. 14	447	0.30648	0.98646	0.97205	0.29389	False
	No. 15	318	0.24842	0.52587	0.65521	0.08559	False
	No. 16	279	0.87813	0.97943	0.95880	0.82464	True
	No. 17	49	0.91837	0.72605	0.19707	0.13140	False
	No. 18	48	0.33333	0.32933	0.30258	0.03321	False

The Nos. 1–18 represent different objects after Euclidean clustering and N indicates the number of points belonging to one point cloud; the bold data indicate the optimal value in the current column; the item “IsTrue” indicates whether it is the calibration board.

, can be obtained. It can be found that objects with many points, shown in Figure 10a,b, will lead to the reduction of overall final Score due to the reduction of ScoreA. For objects with a small number of point clouds, the final Score will reduce due to the decrease of ScoreB and ScoreC, shown in Figure 10d–i. Only our target board, whose ScoreA, ScoreB, and ScoreC all maintain the value close to 1.0, has the highest final Score. It can be said that our locating index is not affected by the number and shape of point clouds. For this reason, in Euclidean clustering, there is no need to specify the parameter of the minimum number of cluster point clouds. In this experiment, we set the minimum number of clusters as 30. Therefore, our method does not require the user to set additional parameters. In addition, the distance parameter of clustering needs to be specified additionally. The distance between LiDAR and the calibration board is often placed between 5.0–10.0 m. Within this distance range, the distance between adjacent points is not sensitive to the change of distance between the calibration board and the LiDAR origin point. We can set it as 30.0 cm, which requires that our calibration board cannot be interfered by other objects within 30.0 cm around during the calibration process.

• Experiment 3: Re-projection error analyses

This section is the evaluation of calibration accuracy. We compare the proposed method with other comparable methods. For each method, we select at least 100 pairs of feature points for a re-projection error test. Figure 11a–r are the error distribution diagrams of three methods. We can show the error distribution of different methods more clearly in the form of error circles. Figure 12a–o show the experimental results of different frame numbers.

Table 4. Calibration results of different methods.

Re-Projection Result
Method	θ (°)	Error (pixel)	Improved (pixel)	Gain (%)	Intrinsic Parameters	Extrinsic Parameters	Stability	Automation Level
Proposed method	0.1	0.9350	-	-	√	√	High	High
	0.2	1.0942	-	-
	0.4	1.1992	-	-
[18]	0.1	2.7422	2.2021	+19.7	√	√	Low	Low
	0.2	3.2414	2.4074	+25.7
	0.4	5.9693	3.0073	+49.6
[22]	0.1	2.1751	1.0416	+52.1	×	√	Middle	Middle
	0.2	2.9751	2.0465	+31.2
	0.4	4.8978	3.3172	+32.3
[19]	0.1	2.2684	1.6406	+27.7	√	√	Middle	Low
	0.2	2.6452	1.9234	+27.3
	0.4	4.4114	3.7310	+15.4
[20]	0.1	1.3317	-	-	×	√	Middle	Low
	0.2	1.3982	-	-
	0.4	3.1492	-	-
[24]	0.1	0.7433	-	-	×	√	Middle	High
	0.2	2.4157	-	-
	0.4	3.5512	-	-

is the qualitative and quantitative analysis comparing the results.

After testing a large number of feature points, our proposed method has a maximum error of no more than 2.0 pixels at each angle resolution, shown in Figure 11a–c. The method of [18] used half of the average distance between adjacent scanning points, as a new position of an edge point was greatly affected by θ. It can also show up from the error circles of Figure 11d–f, which is not suitable for using in LiDAR with large θ. It can be observed in the re-projection error analysis that the projection error is about five times of ours in Table 4. On the one hand, it is related to the poor estimation of edge points, and on the other hand, it is related to the linear equation solution without considering the camera distortion. We apply the approximation edge point fitting method to solve corner points and solve the calibration parameters with the same linear equation. Finally, we can find that the re-projection error is significantly reduced, shown in Table 4. Ref. [22] used a two-stage calibration method, and the deviation of camera intrinsic parameter calibration affected the final calibration effect. In addition, the feature point pairs of [22] are manually selected, which is extremely inconvenient for large-scale applications. Ref. [19] obtained the point clouds at different positions by means of distance statistical filtering, which to a certain extent solved the error caused by LiDAR in ranging to the calibration system. However, the amount of data processed is too heavy (300 frames at least). Compared with the method that we only use one frame of data, our method is obviously better. Similarly, we use our edge refinement method to replace the steps of extracting corner points in [19,22], and conduct experiments on the re-projection error of corner points. It can be found in Table 4 that the accuracy had been improved. Ref. [22] ignored the influence of θ on edge estimation, whilst Ref. [19] only considered the accuracy of edge estimation from one side. Our method considers four sides of the calibration board at the same time, and takes the actual size of the four sides as the reference index for fine processing, which has a higher accuracy. Through the comparison experiment with [18,19,22], we can find that improving the accuracy of corner points can effectively improve the accuracy of calibration. The edge refinement method proposed in this paper is helpful to improve the accuracy of corner points. The work of [20] applied two chessboards and one auxiliary calibration object to calibrate the extrinsic parameters of the camera and LiDAR system, which was not convenient. Ref. [24] had high accuracy, and the 3D feature points obtained were the same as the actual physical size. However, when applied to LiDAR with large θ, it is not as stable as our method. This is because when θ is large, the distance between the edge point of the calibration board and the real edge is often large, which will affect the matching accuracy. Finally, we make a qualitative evaluation of various methods in Table 4. We can clearly find from Figure 11a–c that, using our method, the re-projection error will not change significantly with the increase of θ, and the overall error remains within 2.0 pixels (red circles). In general, our method is superior to other methods in stability and automation level.

Figure 12a–o analyze the re-projection results under different numbers of frames. It is obvious that other methods need at least 3–5 frames of data to converge, and with the increase of θ, the error of final convergence will also increase. The proposed method has a low re-projection error in the first frame, and with the increase of the number of frames, the error change is not significant. By increasing the position of the calibration board, more data can be obtained. We use 2–10 positions to carry out the calibration experiment, shown in

Table 5. Mean and standard deviation of re-projection errors for 2–10 groups of images and point clouds.

θ(°)	NP	2	3	4	5	6	7	8	9	10
0.1	Mean	1.535	1.469	1.426	1.301	1.283	1.274	1.272	1.249	1.21
	Std	0.696	0.703	0.549	0.658	0.565	0.529	0.547	0.492	0.476
0.2	Mean	2.109	1.945	1.637	1.470	1.457	1.408	1.324	1.303	1.256
	Std	1.007	0.682	0.740	0.623	0.697	0.644	0.737	0.624	0.552
0.4	Mean	2.301	1.858	1.750	1.745	1.620	1.453	1.388	1.311	1.28
	Std	1.372	0.720	0.712	0.630	0.720	0.802	0.790	0.755	0.586

NP represents number of test positions; Mean is the average error; and Std means the standard deviation.

. It can be clearly found that with the increase of data volume, the average value and variance of re-projection error are reduced, which also shows the robustness of the proposed method. Our calibration board can add 20 more pairs of feature points without adding a position, which is larger than the data volume of other methods, and the results are more convincing.

• Experiment 4: Parameter consistency check

The calibration of LiDAR and camera is a no ground true problem, so the measurement means of calibration effect is not only the re-projection error, but also the data consistency. We collect data from a total of 100 positions, randomly select data from 10 positions for calibration experiments, randomly select 100 times, and record the intrinsic parameters {$f_x,f_y,u_x,v_y,k_1,k_2$} of camera and extrinsic parameters {$R\left\{r_x,r_y,r_z\right\}, T\left\{t_x,t_y,t_z\right\}$}. The mean and Std are shown in the following

Table 6. Mean and Std of optimization parameters.

θ(°)		fx	fy	ux(pixel)	vy(pixel)	k1	k2	rx(°)	ry(°)	rz(°)	tx(cm)	ty(cm)	tz(cm)
0.1	Mean	2825.75	2817.57	969.026	597.901	−0.220	0.187	89.767	−89.738	−0.321	−6.0651	9.6231	−1.4969
	Std	1.947	2.382	6.471	2.873	0.001	0.010	0.042	0.032	0.103	0.5986	0.3348	0.1995
0.2	Mean	2826.23	2818.55	970.405	595.841	−0.230	0.161	90.001	−89.775	−0.582	−6.0065	10.0176	−1.5216
	Std	2.827	3.463	5.824	3.604	0.009	0.000	0.069	0.035	0.111	0.4551	0.3172	0.4093
0.4	Mean	2826.85	2818.18	966.150	602.771	−0.236	0.201	90.106	−89.690	−0.624	−5.9884	8.4497	−1.7355
	Std	3.146	3.955	9.0291	5.042	0.007	0.001	0.044	0.042	0.171	0.4098	0.4856	0.4395

; the change of intrinsic parameters is less than 1.0%, the change of rotation matrix is about 0.1°, and the change of translation matrix is about 4.0 mm, which indicates that the calibration parameters of our method are consistent.

• Experiment 5: Qualitative analysis of calibration accuracy

The following experiment is a qualitative measurement of calibration effect. In order to verify the detection effect of UAVs under low-altitude remote sensing, we select different reference objects 3.0–100.0 m away from the LiDAR for qualitative evaluation. We project the point clouds onto the image plane and color them according to the distance. Figure 13 is the projection of the indoor environment and Figure 14, Figure 15 and Figure 16 are the projection of the outdoor environment. The farthest projection distance in Figure 14, Figure 15 and Figure 16 is about 100.0 m, which is enough to explain that our method has not offset in this distance. The buildings and vehicles can also show that the proposed method can add accurate depth information for nearby objects. The 3D reconstruction effect of Figure 17 also shows that the proposed calibration method is helpful to the advanced application of the LiDAR–camera system.

4.3. Discussions

Low-altitude UAVs have been used in many remote sensing fields [37] in recent years such as urban management, crop monitoring, and industrial patrol inspection, etc. It is essential to apply multi-sensor fusion technology to low-altitude UAVs. Calibration is the basic requirement for multi-sensor platforms where data need to be represented in a common coordinate system for the purpose of analysis and information fusion [38]. The target-based calibration method has become the first choice for the calibration methods because of its stability and high accuracy. This paper designs a special calibration target and its corresponding calibration algorithm, which can achieve high precision and great calibration efficiency. While considering the calibration accuracy, we also consider the implementation efficiency of this method. The fully automatic data processing flow and joint optimization scheme make this method suitable for large-scale commercial applications.

At present, there are few methods to consider the effect of LiDAR with different θs on calibration. Experiment 3 in this paper shows that LiDAR with different θs has significant differences in calibration. The method of fine edge processing based on background point clouds proposed in this method is to solve the problem of insufficient calibration accuracy caused by large θ. After using the proposed method, the calibration accuracy and robustness have been improved to some extent, as analyzed in Table 4, and the demand for data volume has also been decreased, as shown in Figure 12 and Table 5. In addition, the method in this paper focuses on large-scale application, so every step of the design is to minimize the artificial setting of parameters, so as to improve the implementation efficiency of calibration. The average processing time of this method for one frame of data, including point cloud and image, is 7.5 s. It has higher time complexity than the traditional method, which is one aspect that can be optimized in the future. However, we are more focused on the improvement of calibration accuracy and automation level. Higher accuracy and less human participation are our final goals. In terms of optimizing functions, we need to give the initial extrinsic parameters in advance, which can be easily given by PnP [39,40].

Here, we discuss the advantages of the calibration board we designed. We have designed a special hollow monochromatic plane plate, which has the following advantages: (1) the monochromatic flat plate avoids the interference of different colors when the LiDAR is ranging, which is conducive to plane fitting and plane normal vector extraction. (2) The hollow calibration object is convenient for automatically locating its position in the point clouds by using geometric information. As for these two advantages, we conduct the following comparative experiments. We paste a small piece of checkerboard on the calibration object and find that the depth discontinuity with other positions appears at the place where the LiDAR scans, shown in Figure 18. This phenomenon is related to the intensity information performance of LiDAR itself. When we use the LiDAR to detect the calibration object, the chessboard in black and white areas will make the ranging depth uneven [18,19], while the monochromatic calibration board can reflect a flat plane. Our method does not depend on the intensity information [41] or additional sensors [42], and only uses geometric information to locate the calibration board; thus, it can be used with any types of LiDAR device and has better universality.

The advantages of the method described in this paper are as follows: (1) The refined edge points enable the LiDAR with large horizontal angle resolution (θ = 0.4) to complete calibration with high accuracy; (2) The two-step corner points detection method designed by us can realize fully automatic data processing, greatly reduce manual participation, and does not need parameters to be manually set in the calibration process; (3) Our method does not require a high amount of calibration data. For 5–10 positions, one frame of data per position can complete the calibration with high accuracy. This is also related to the results of our 2D and 3D feature point refinement processing because the feature points with high accuracy can reduce the amount of data involved in optimization. The experimental scheme in this paper also has some shortcomings. For example, our research is based on LiDAR with 64 lines, and has not been verified on the LiDAR with fewer lines (32 lines or less) [43,44], and the vertical angle resolution needs to be further considered, which requires further research and improvements. In future work, we will apply the results of joint calibration of LiDAR and visible light camera to the application of low-altitude UAVs remote sensing, ground object recognition, depth information completion, and 3D terrain reconstruction, etc. Our application scenario is to use UAV with a visible light camera and LiDAR to conduct high-precision 3D mapping in a small range. The joint calibration of these sensors is the pre-work of this application. In this paper, this calibration method is completed in an indoor environment, as shown in Figure 2a. As for the verification of an outdoor environment, in order to show the projection effect at different distances, we chose the scene with very rich features, the nearest distance is 3.0 m, and the farthest distance is 100.0 m, as shown in Figure 14, Figure 15 and Figure 16. The final projection effect shows that the projection of our method at different distances is accurate. As for the limitations of our method on the application level, the efficiency of the algorithm is lower than that of a pure image algorithm, which is also related to the increase of 3D data. In addition, we pay more attention to the ground object detection within 100.0 m of low-altitude, and do not consider the objects with a distance greater than 100.0 m, which is related to the inherent properties of our LiDAR.

5. Conclusions

In this paper, we propose a novel, fully automatic joint calibration method of LiDAR and visible light camera. Compared with the existing methods, we propose an approximate edge fitting technology, which has obvious improvements based on the existing methods and has definite accuracy improvement. Our method does not need to use the intensity information of LiDAR, which has better universality for the method itself and can also be applied when the LiDAR’s horizontal resolution is sparse, or the intensity information is poor. For the automatic process, we propose a strategy of first locating and then refinement, which can fully realize 2D and 3D corner points detection. In the calibration process, the user does not need manual intervention, which greatly increases the efficiency of calibration. In the aspect of effect display, we project the 3D point clouds onto the image plane and show the calibration effect in different test distances. The 2D image with depth information can be used for advanced visual applications. We also color the 3D point cloud according to the calibrated parameters, which is instructive for 3D reconstruction.