CN103729846B - LiDAR point cloud data edge detection method based on triangular irregular network - Google Patents

Abstract

The invention discloses a LiDAR point cloud data edge detection method based on a triangular irregular network. The method comprises the steps of reading original LiDAR point cloud data of an object region, eliminating gross error noise points, conducting Delaunay triangulation on the LiDAR point cloud data with the gross error noise points eliminated, calculating and saving the deformation quantity of each spatial triangle in the triangular irregular network, calculating deformation quantity threshold values of the triangles, marking the spatial triangles with the corresponding deformation quantity larger than the corresponding deformation quantity threshold value as marginal triangles, calculating long-narrow triangles, comparing the long-narrow degrees Li with long-narrow degree threshold values of the triangles, and determining marginal points and obtaining an edge image of the LiDAR point cloud data. According to the LiDAR point cloud data edge detection method based on the triangular irregular network, elevation-jump edge point detection is conducted by means of triangle deformation during triangulation, blank area edge point detection is conducted by means of the long-narrow triangles, and edge detection in original data is achieved to reserve more information of the LiDAR point cloud data.

Description

LiDAR point cloud data edge detection method based on TIN

Technical field

The invention belongs to remote sensing application research field, be specifically related to a kind of LiDAR point cloud number based on TIN According to edge detection method.

Background technology

Airborne LiDAR system can directly obtain Three Dimensional Ground data, has high accuracy, high density, high efficiency and low The advantage of cost, the available required various image products of mass data utilizing it quickly to obtain, therefore at Modern Surveying & Mapping In play an increasingly important role.Image border is one of key property of image, is to compare valuable letter in image Breath.The marginal information detected is provided precondition as characteristic information, registration and fusion for image.Therefore, Utilize airborne LiDAR system quick obtaining high-precision LiDAR point cloud data, and detect that preferable edge has weight The meaning wanted.LiDAR point cloud data refer to that airborne LiDAR system obtains earth's surface by launching and receiving laser pulse 3-dimensional point coordinates in high precision.

The edge detection method of traditional LiDAR point cloud data, is usually and discrete LiDAR point cloud data is advised Then grid resampling, obtains digital depth (elevation) matrix, more different depth values carries out grey level quantization and stretching change Changing, obtain the distance gray scale image that same gray level image is the same, this process is referred to as " Range Imaging ".Lai Xudong, Wan Youchuan Carry out medium filtering Deng to depth image, utilize Robert operator, Sobel operator and Prewitt operator to carry out edge and carry Take；Wang great Ying, Cheng Xinwen, Pan Huibo etc. use optimal threshold to carry out after binaryzation depth image to use Canny and Log operator carries out rim detection；Xu Jingzhong, ten thousand children rivers etc. utilize adaptive medium filtering and adaptive to depth image The Canny operator answering threshold value carries out rim detection；It is swollen that depth image is utilized in mathematical morphology by Wu Hangbin, Liu Chun etc. Swollen and caustic solution carries out sequential computing.And the image obtained is carried out edge extracting and edge vectorization, obtain eachly Edge corresponding to thing and data point.

By the edge detection method of above-mentioned various LiDAR point cloud data is analyzed, find that they there are as follows Defect: carry out during Range Imaging to bring error during resampling so that some rivers or the limit of some cloud white space Edge cannot be detected；Can bring round-off error when grey level quantization and stretching conversion, the point that elevation difference is less is entering Same gray level may be quantified as during row grey level quantization or gray difference is the least, it is also possible to by elevation less for difference It is quantified as different gray levels, therefore all can bring certain error when rim detection so that some endpoint detections are not To or obtain mistake marginal point.

To sum up, although existing carrying out the edge detection method of LiDAR point cloud data based on depth image and can reach certain Effect, but but lost some data messages during depth image, the marginal information obtained is sufficiently complete, Therefore the marginal information obtained can not obtain other useful informations such as grade of LiDAR point cloud data, is unfavorable for follow-up place Reason.

Summary of the invention

The deficiency of rim detection is carried out, purpose herein for the existing depth image that LiDAR point cloud data is converted to It is, it is provided that a kind of LiDAR point cloud data edge detection method based on TIN, the method is directly former In beginning data, difference condition based on TIN intermediate cam shape deformation quantity carries out LiDAR point cloud data edges detection. The method makes full use of the endpoint detections that the change of shape of triangle during triangulation carries out having elevation sudden change, profit The endpoint detections of data white space is carried out, it is achieved that in initial data, carry out rim detection reservation with long-narrow triangular mesh The more information of LiDAR point cloud data.Present invention could apply to airborne LiADR cloud data rim detection.

In order to realize above-mentioned task, the present invention adopts the following technical scheme that and is solved:

A kind of LiDAR point cloud data edge detection method based on TIN, specifically includes following steps:

Step 1: read target area original LiDAR point cloud data；

Step 2: excluding gross error noise spot；

Step 3: the LiDAR point cloud data after excluding gross error noise step 2 obtained carry out Delaunay triangle and cut open Point, generate the TIN in xoy plane；Plus forming xyz after the elevation information of each point on TIN Space TIN TIN；

Step 4: calculate the deformation quantity of each spatial triangle in TIN and preserve；

Step 5: calculate the deformation quantity threshold value of triangle；

Step 6: deformation quantity step 4 obtained is labeled as edge triangles more than the spatial triangle of deformation quantity threshold value； If deformation quantity is more than deformation quantity threshold value, enters step 8, otherwise enter step 7；

Step 7: calculate long-narrow triangular mesh, the relatively long and narrow degree L of each triangle_iWith narrow length threshold Thresd, if L_iMore than Thresd, then this triangle is long-narrow triangular mesh, is edge triangles, wherein, Thresd by this triangular marker Not less than 100；

Step 8: determine marginal point, obtains the edge image of LiDAR point cloud data, specifically comprises the following steps that

801: the determination of thick edge point:

The each edge triangles obtained for step 6 utilizes formula 8 to calculate elevation threshold value H (G), and respectively to each limit Summit in edge triangle carries out choosing of thick edge point:

$H\left(G\right)=\frac{H\left(A\right)+H\left(B\right)+H\left(C\right)}{3}---\left(8\right)$

Wherein, H (A), H (B), H (C) are the elevation information of A, B, C point respectively, unit: m；H (G) is Elevation threshold value, unit: m；

If the three of edge triangles summit A, B, C have any two points A, the elevation of B be all higher than H (G) and ≤ σ, σ take 0.5-2.5 to | H (A)-H (B) |, then by A, B 2 is all taken as thick edge point, otherwise, by A, in B, C 3 The point of elevation maximum is taken as thick edge point；

Three summits of each edge triangles obtaining step 7 are defined as thick edge point；

: 802: delete the isolated point in all thick edges point, using remaining thick edge point as final marginal point output, I.e. obtain the edge image of LiDAR point cloud data.

Further, the rough error noise spot of described step 2 utilizes formula 1 to be calculated:

N=∑ (| P_i-P_j| ＜ σ) (1)

Wherein, P_i, P_jRepresenting thick edge point midpoint i and the coordinate of some j, n represents in original LiDAR point cloud data With coordinate P_iThe distance number of point less than σ, n takes 0-10, σ and takes 10-20m.

Further, the deformation quantity of each spatial triangle in TIN that calculates of described step 4 refers to: to often Individual triangle utilizes formula 2 to calculate:

$s=C\×\underset{i=1}{\overset{3}{\Σ}}|\sin A_i^{\′}-\sin A_i|---\left(2\right)$

Wherein, s is the deformation quantity of spatial triangle, and C is constant, takes the arbitrary value more than 1, A_iRepresent that xoy puts down The interior angle of face triangle, A '_iWith A after expression interpolation elevation information formation spatial triangle_iCorresponding interior angle, i=1,2,3.

Further, described step 5 calculates the specifically comprising the following steps that of deformation quantity threshold value of triangle

501: read the deformation quantity set S of the deformation quantity s composition of all spatial triangles；Set deformation quantity threshold value T^k's Initial value T⁰:

$T^0=\frac{S_{\min }+S_{max}}{2}---\left(3\right)$

Wherein, S_minAnd S_maxRepresent the minima in deformation quantity set S and maximum respectively；

502: triangle deformation quantity s all of in deformation quantity set S is divided into deformation quantity more than T^kSet A and deformation Amount is less than or equal to T^kSet B two parts；

503: utilize formula 4, formula 5 to calculate A, B two-part averaged deformation amount g_AAnd g_B:

$g_A=\frac{1}{n}\underset{S\left(i\right)>T^k}{\Σ}S\left(i\right)---\left(4\right)$

$g_B=\frac{1}{m}\underset{S\left(i\right)\≤T^k}{\Σ}S\left(i\right)---\left(5\right)$

Wherein, S (i) is the i-th triangle deformation quantity in deformation quantity set S, n and m is A and B intermediate cam respectively The number of shape deformation quantity s；

504: utilize formula 6 to calculate new threshold value T^k+1:

$T^{k+1}=\frac{g_A+g_B}{2}---\left(6\right)$

505: Rule of judgment | T^k+1-T^k|≤δ, k=0,1,2 ... whether set up, δ is allowable error, herein error model Enclosing C*0.01-C*0.05, wherein, C is the constant in formula 2；Condition is set up, then deformation quantity threshold value T^kIt is updated to step 504 T obtained^k+1, then perform step 6；Otherwise, threshold value is still T^k, this seasonal k=k+1, and return step 503。

Further, the calculating long-narrow triangular mesh of described step 7 refers to:

Utilize the long and narrow degree of all trianglees in the space triangular net that formula 7 calculation procedure 3 obtains；

L_i=| A_iG_i|²+|B_iG_i|²+|C_iG_i|² (7)

Wherein, L_iRepresent the long and narrow degree of i-th triangle in the triangulation network；A_i, B_i, C_iIt is respectively i-th triangle Three summits, G_iFor i-th barycenter oftriangle.

Further, isolated point described in the 802 of described step 8 utilizes the formula 9 of step 2 to calculate:

N=∑ (| P_i-P_j| ＜ σ) (9)

Wherein, P_iRepresent thick edge point midpoint i XY coordinate, n represent in thick edge point with coordinate P_iDistance little In the number of the point of σ, n takes 0-10, σ and takes 10-20m.

Compared with prior art, the invention have the advantages that

First, the present invention is to when having the atural object edge of elevation sudden change, and application Delaunay Triangulation is set up irregular The triangulation network, utilizes the deformation quantity size judgment criterion as edge triangles of triangle, by deformation quantity more than the three of threshold value Corner mark position edge triangles, processes edge triangles and obtains marginal point.The method is without by LiDAR point cloud number According to generating depth image, thus the MARG when being converted into depth image is avoided in existing edge detection method to damage Lose, improve precision edge being detected, be also beneficial to follow-up LiDAR data feature extraction and and remote sensing images Fusion with SAR image.

Second, the present invention, when the LiDAR point cloud data such as data white space such as having river are carried out rim detection, utilizes White space edge, as estimating, is detected by long-narrow triangular mesh, eliminate during Data grid loss of data with And the marginal information loss that interpolation causes, maintain the integrity of marginal information.

Accompanying drawing explanation

Fig. 1 is the general flow chart of the LiDAR point cloud data edge detection method based on TIN of the present invention.

Fig. 2 is the design sketch of the l-G simulation test 1 that the method for the present invention detects for LiDAR point cloud data edges.Wherein, Fig. 2 (a) is original LiDAR point cloud data image；Fig. 2 (b) is the image after original image carries out triangulation；Figure 2 (c) is the endpoint detections figure that the method for the present invention obtains；Fig. 2 (d) is by initial data and the design sketch of marginal point superposition.

Fig. 3 is the method l-G simulation test 2 for the data edges detection containing data white spaces such as rivers of the present invention Design sketch.Wherein, Fig. 3 (a) is original LiDAR point cloud data image；Fig. 3 (b) is that original image is carried out triangle Image after subdivision；Fig. 3 (c) is the endpoint detections figure that the method for the present invention obtains；Fig. 3 (d) is by initial data and limit The design sketch of edge point superposition.

Below in conjunction with the drawings and specific embodiments, the present invention is further explained.

Detailed description of the invention

With reference to Fig. 1, the LiDAR point cloud data edge detection method based on TIN of the present invention, specifically include Following steps:

Step 1: read target area original LiDAR point cloud data；

Step 2: excluding gross error noise；

Original point cloud data can produce some rough error noises, in order to obtain cloud data accurately, first has to rough error noise Point is rejected, and rough error noise spot utilizes formula 1 to determine:

N=∑ (| P_i-P_j| ＜ σ) (1)

Wherein, P_i, P_jRepresenting thick edge point midpoint i and the coordinate of some j, n represents in original LiDAR point cloud data With coordinate P_iThe distance number of point less than σ, n takes 0-10, σ and takes 10-20m；

Step 3: the LiDAR point cloud data after excluding gross error noise step 2 obtained carry out Delaunay triangle and cut open Point, the discrete data point that will be randomly distributed generates the TIN in xoy plane；On TIN Plus forming xyz space TIN (TIN) after the elevation information of each point；

Satisfied following 3 criterions of Delaunay Triangulation:

(1) the minimum angle of all trianglees in triangulation T is maximized, to avoid the occurrence of the most flat triangle；

(2) circumscribed circle of arbitrary triangle of triangulation T is empty circle, does not i.e. comprise any other point in this circumscribed circle, Reach local optimum to process；

(3) the angle sum on the subdivision diagonal both sides that tetragon is corresponding is not more than 180 °.

Step 4: utilize following formula calculate the deformation quantity of each spatial triangle in TIN and preserve；

$s=C\×\underset{i=1}{\overset{3}{\Σ}}|\sin A_i^{\′}-\sin A_i|---\left(2\right)$

Wherein, s is the deformation quantity of spatial triangle, and C is that (C is the span in order to increase s to constant, takes and is more than The arbitrary value of 1), A_iRepresent the interior angle of xoy plane triangle, A '_iRepresent that adding elevation information forms spatial triangle Afterwards with A_iCorresponding interior angle, i=1,2,3.

Step 5: calculate the deformation quantity threshold value of triangle, specifically comprise the following steps that

501: read the deformation quantity set S of the deformation quantity s composition of all spatial triangles；Set deformation quantity threshold value T^k's Initial value T⁰:

$T^0=\frac{S_{\min }+S_{max}}{2}---\left(3\right)$

Wherein, S_minAnd S_maxRepresent the minima in deformation quantity set S and maximum respectively；

502: triangle deformation quantity s all of in deformation quantity set S is divided into deformation quantity more than T^kSet A and deformation Amount is less than or equal to T^kSet B two parts；

503: utilize formula 4, formula 5 to calculate A, B two-part averaged deformation amount g_AAnd g_B:

$g_A=\frac{1}{n}\underset{S\left(i\right)>T^k}{\Σ}S\left(i\right)---\left(4\right)$

$g_B=\frac{1}{m}\underset{S\left(i\right)\≤T^k}{\Σ}S\left(i\right)---\left(5\right)$

Wherein, S (i) is the i-th triangle deformation quantity in deformation quantity set S, n and m is in A and B three respectively The number of dihedral deformation quantity s；

504: utilize formula 6 to calculate new threshold value T^k+1:

$T^{k+1}=\frac{g_A+g_B}{2}---\left(6\right)$

505: Rule of judgment | T^k+1-T^k|≤δ, k=0,1,2 ... whether set up, δ is allowable error, herein error Scope C*0.01-C*0.05, wherein, C is the constant in formula 2；Condition is set up, then deformation quantity threshold value T^kIt is updated to The T that step 504 obtains^k+1, then perform step 6；Otherwise, threshold value is still T^k, this seasonal k=k+1, and return Return step 503.

Step 6: deformation quantity step 4 obtained is more than deformation quantity threshold value T^kSpatial triangle be labeled as edge triangle Shape；If deformation quantity is more than deformation quantity threshold value, enters step 8, otherwise enter step 7；

Step 7: calculating long-narrow triangular mesh:

Utilize the long and narrow degree of all trianglees in the space triangular net that formula 7 calculation procedure 3 obtains；Relatively each triangle is narrow Length L_iWith narrow length threshold Thresd, if L_iMore than Thresd, Thresd not less than 100；Then this triangle is Long-narrow triangular mesh, is edge triangles by this triangular marker；

L_i=| A_iG_i|²+|B_iG_i|²+|C_iG_i|² (7)

Wherein, L_iRepresent the long and narrow degree of i-th triangle in the triangulation network；A_i, B_i, C_iIt is respectively i-th triangle Three summits, G_iFor i-th barycenter oftriangle.

Step 8: determine marginal point, obtains the edge image of LiDAR point cloud data, specifically comprises the following steps that

801: the determination of thick edge point:

The each edge triangles obtained for step 6 utilizes formula 8 to calculate elevation threshold value H (G), and respectively to each limit Summit in edge triangle carries out choosing of thick edge point:

$H\left(G\right)=\frac{H\left(A\right)+H\left(B\right)+H\left(C\right)}{3}---\left(8\right)$

Wherein, H (A), H (B), H (C) are the elevation information of A, B, C point respectively, unit: m；H (G) is Elevation threshold value, unit: m；

If the three of edge triangles summit A, B, C have any two points A, the elevation of B be all higher than H (G) and ≤ σ, σ take 0.5-2.5 to | H (A)-H (B) |, then by A, B 2 is all taken as thick edge point, otherwise, by A, in B, C 3 The point of elevation maximum is taken as thick edge point.

Three summits of each edge triangles obtaining step 7 are defined as thick edge point.

: 802: deleting the isolated point in all thick edge points obtained, isolated point utilizes the formula 9 of step 2 to calculate； Using remaining thick edge point as final marginal point output, i.e. obtain the edge image of LiDAR point cloud data.

N=∑ (| P_i-P_j| ＜ σ) (9)

Wherein, P_i, P_jRepresent thick edge point midpoint i and some j coordinate, n represent in thick edge point with coordinate P_i's The number of the distance point less than σ, n takes 0-10, σ and takes 10-20m.

In order to the effect of the present invention is described, inventor has carried out following l-G simulation test.

L-G simulation test 1:

The LiDAR point cloud data blank to no data in the present invention carry out rim detection simulated effect.Simulated conditions: MATLAB7.0 software.Follow the technical scheme of the invention described above, see Fig. 2, normal data Sample23 is carried out Rim detection emulates.In this test, in step 2, the value of n is 6, and the value of σ is 15m；Step 4 and step 5 In the value of constant C be 1000；In step 5, the value of δ is 10；In step 6, the threshold value of long-narrow triangular mesh is 100； In step 7, the value of σ is 1.5m.

From Fig. 2 (d) it can be seen that the method for the present invention remains the full detail of LiDAR point cloud, follow-up Feature extraction, registration fusion etc. processes.And utilize gray-scale map after LiDAR point cloud data are carried out depth image Some data can be lost as carrying out the method for rim detection, some marginal informations can be reduced during gray processing, and do not have Keep the effective information of LiDAR point cloud data, be unfavorable for follow-up process.

L-G simulation test 2:

The LiDAR data containing the data white spaces such as river in the present invention is carried out the emulation of rim detection.Simulated conditions: MATLAB7.0 software.Follow the technical scheme of the invention described above, with reference to Fig. 3, normal data Sample61 is carried out edge Detection simulation.In this test, in step 2, the value of n is 6, and the value of σ is 15m；Constant C in step 4 and step 5 Value be 1000；In step 5, the value of δ is 10；In step 6, the threshold value of long-narrow triangular mesh is 100；σ in step 7 Value is 1.5m.

From Fig. 3 (d) it can be seen that the method for the present invention can be good at detecting the marginal point of data white space, profit River marginal point is detected as estimating with long-narrow triangular mesh.And it is laggard that LiDAR point cloud data are carried out depth imageization White space can be eliminated by row interpolation, and the edge detection method of recycling gray level image there will be when detecting and can't detect Edge, river or the edge of some mistakes detected.

Claims (6)

1. a LiDAR point cloud data edge detection method based on TIN, it is characterised in that specifically include following steps:

Step 1: read target area original LiDAR point cloud data；

Step 2: excluding gross error noise spot；

Step 3: the LiDAR point cloud data after excluding gross error noise step 2 obtained carry out Delaunay Triangulation, generates the TIN in xoy plane；Plus forming xyz space TIN TIN after the elevation information of each point on TIN；

Step 4: calculate the deformation quantity of each spatial triangle in TIN and preserve；

Step 5: calculate the deformation quantity threshold value of triangle；

Step 6: deformation quantity step 4 obtained is labeled as edge triangles more than the spatial triangle of deformation quantity threshold value；If deformation quantity is more than deformation quantity threshold value, enters step 8, otherwise enter step 7；

Step 7: calculate long-narrow triangular mesh, the relatively long and narrow degree L of each triangle_iWith narrow length threshold Thresd, if L_iMore than Thresd, then this triangle is long-narrow triangular mesh, is edge triangles by this triangular marker, and wherein, Thresd is not less than 100；

Step 8: determine marginal point, obtains the edge image of LiDAR point cloud data, specifically comprises the following steps that

801: the determination of thick edge point:

The each edge triangles obtained for step 6 utilizes formula 1 to calculate elevation threshold value H (G), and the summit in each edge triangles carries out thick edge point choose respectively:

Wherein, H (A), H (B), H (C) are the elevation information of A, B, C point respectively, unit: m；H (G) is elevation threshold value, unit: m；

If the three of edge triangles summit A, B, having any two points A in C, the elevation of B is all higher than H (G) and | H (A)-H (B) |, and≤σ, σ take 0.5-2.5m, then by A, B 2 is all taken as thick edge point, otherwise, by A, the point of 3 middle elevations maximums of B, C is taken as thick edge point；

Three summits of each edge triangles obtaining step 7 are defined as thick edge point；

802: delete the isolated point in all thick edges point, using remaining thick edge point as final marginal point output, i.e. obtain the edge image of LiDAR point cloud data.

2. LiDAR point cloud data edge detection method based on TIN as claimed in claim 1, it is characterised in that the rough error noise spot of described step 2 utilizes formula 2 to be calculated:

N=∑ (| P_i-P_j| ＜ σ) (2)

Wherein, P_iRepresent the XY coordinate of original LiDAR point cloud data midpoint i, n represent in original LiDAR point cloud data with coordinate P_iThe distance number of point less than σ, n takes 0-10, σ and takes 10-20m.

3. LiDAR point cloud data edge detection method based on TIN as claimed in claim 1, it is characterized in that, the deformation quantity of each spatial triangle in TIN that calculates of described step 4 refers to: utilize formula 3 to calculate each triangle:

Wherein, s is the deformation quantity of spatial triangle, and C is constant, takes the arbitrary value more than 1, A_iRepresent the interior angle of xoy plane triangle, A '_iWith A after expression interpolation elevation information formation spatial triangle_iCorresponding interior angle, i=1,2,3.

4. LiDAR point cloud data edge detection method based on TIN as claimed in claim 1, it is characterised in that described step 5 calculates specifically comprising the following steps that of the deformation quantity threshold value of triangle

501: read the deformation quantity set S of the deformation quantity s composition of all spatial triangles；Set deformation quantity threshold value T^kInitial value T⁰:

Wherein, S_minAnd S_maxRepresent the minima in deformation quantity set S and maximum respectively；

502: triangle deformation quantity s all of in deformation quantity set S is divided into deformation quantity more than T^kSet A and deformation quantity less than or equal to T^kSet B two parts；

503: utilize formula 5, formula 6 to calculate A, B two-part averaged deformation amount g_AAnd g_B:

Wherein, S (i) is the i-th triangle deformation quantity in deformation quantity set S, n and m is the number of A and B intermediate cam shape deformation quantity s respectively；

504: utilize formula 7 to calculate new threshold value T^k+1:

505: Rule of judgment | T^k+1-T^k|≤δ, k=0,1,2 ... whether setting up, δ, herein range of error C*0.01-C*0.05, wherein, C is the constant in formula 2 if being allowable error；Condition is set up, then deformation quantity threshold value T^kIt is updated to the T that step 504 obtains^k+1, then perform step 6；Otherwise, threshold value is still T^k, this seasonal k=k+1, and return step 503.

5. LiDAR point cloud data edge detection method based on TIN as claimed in claim 1, it is characterised in that the calculating long-narrow triangular mesh of described step 7 refers to:

Utilize the long and narrow degree of all trianglees in the space triangular net that formula 8 calculation procedure 3 obtains；

L_i=| A_iG_i|²+|B_iG_i|²+|C_iG_i|² (8)

Wherein, L_iRepresent the long and narrow degree of i-th triangle in the triangulation network；A_i, B_i, C_iIt is respectively three summits of i-th triangle, G_iFor i-th barycenter oftriangle.

6. LiDAR point cloud data edge detection method based on TIN as claimed in claim 1, it is characterised in that isolated point described in the 802 of described step 8 utilizes the formula 9 of step 2 to calculate:

N=∑ (| P_i-P_j| ＜ σ) (9)

Wherein, P_i, P_jRepresent thick edge point midpoint i and some j coordinate, n represent in thick edge point with coordinate P_iThe distance number of point less than σ, n takes 0-10, σ and takes 10-20m.