CN106485737A - Cloud data based on line feature and the autoregistration fusion method of optical image - Google Patents

Abstract

The invention discloses the autoregistration fusion method of a kind of cloud data based on line feature and optical image, including：Cloud data is filtered processing using Mathematical Morphology Method；Being respectively adopted self adaptation supports weight dense stereo algorithm and Delaunay Triangulation algorithm to determine the depth map of optical image and cloud data；Extract the line feature of depth map by Hough algorithm, slightly mated as similarity measure with line length ratio using line feature two wire clamp angle；Reject error matching points pair by two step RANSAC algorithms, obtain camera position parameter estimation；Carry out the color and vein mapping of cloud data and optical image, the 3-D view after being merged.The present invention no relies on to priori such as GPS/INS, and robustness is good, and autoregistration degree is high.

Description

Cloud data based on line feature and the autoregistration fusion method of optical image

Technical field

The invention belongs to a kind of improved heterologous image registration algorithm, particularly a kind of cloud data based on line feature and light Learn the autoregistration fusion method of image.

Background technology

The LiDAR point cloud data registration that according to the registering primitive that used can be divided into based on gray scale registering with optical image Method and the method for registering of feature based.Wherein feature registration is divided into a feature, line feature and region feature registration again.Point Feature be feature based Image registration method in most commonly used, but point feature registration has limitation it is impossible to adapt to universal Commonly used.As soon as line feature is the more higher leveled characters of ground object of analogy point, for the extraction of feature, extract meaningful Line aspect ratio extract significant feature easily, and, all there are rich in natural environment or artificial environment Rich linear ground object, such as road, the edge in river, contour line of target etc..

But traditional line feature registration algorithm has obvious defect.Extract straightway from profile point, using straight line The distribution of the slope of section and edge contour is come to determine the method for transformation parameter be a weight in terms of based on linear feature registration Big breakthrough, is to be laid a good foundation based on the method for registering of linear feature.But concentration, shape are compared for building distribution The more similar city three-dimensional reconstruction of shape, when the similar straightway of a plurality of position of appearance, slope, can cause matching error. A kind of Line feature method, three straight lines being constituted with edges such as building, houses and the correlation of other 3CS Angle asks for similarity measure, finally using two step RANSAC methods, registration result is corrected.This method is completeer Whole comprehensive used line feature, but complexity is higher, operation time is long, and for extracting on remote sensing image of taking photo by plane For being the building at rule of similarity edge, matching characteristic reduces, and error hiding rate can uprise.Algorithm above generally needs Rely on the initial position priori such as GPS/INS, for the registration between three-dimensional, algorithm complex is very high, and depends on three The precision of dimension reconstruct.

Content of the invention

It is an object of the invention to provide a kind of improved heterologous image registration algorithm, particularly a kind of point based on line feature Cloud data and the autoregistration fusion method of optical image.

The technical solution realizing the object of the invention is：

A kind of cloud data based on line feature and the autoregistration fusion method of optical image, comprise the following steps：

Step 1, is filtered to cloud data processing；

Step 2, is respectively adopted self adaptation and supports weight dense stereo algorithm and Delaunay Triangulation algorithm to determine optics Image and the depth map of cloud data；

Step 3, extracts the line feature of depth map by Hough algorithm, using line feature two wire clamp angle and line length ratio as Similarity measure is slightly mated；

Step 4, rejects error matching points pair by two step RANSAC algorithms, obtains camera position parameter estimation；

Step 5, carries out the color and vein mapping of cloud data and optical image, the 3-D view after being merged.

The present invention compared with prior art, its remarkable advantage：

(1) cloud data based on a feature of the present invention and the autoregistration fusion method of optical image are without GPS/INS Etc. priori, high degree of automation, robustness is good；

(2) present invention is carried out on the basis of visible ray depth map, improves the registration accuracy of image；

(3) present invention adopts improved line feature registration, reduces amount of calculation on the premise of ensureing registration accuracy, tool There is higher registering accuracy.

Brief description

Fig. 1 is the autoregistration fusion method flow chart with optical image for the cloud data based on line feature of the present invention.

Fig. 2 (a) is the cloud data of the embodiment of the present invention, and Fig. 2 (b) is visible images, and Fig. 2 (c) is to merge 3-D view afterwards.

Specific embodiment

In conjunction with Fig. 1, the cloud data based on a feature of the present invention and the autoregistration fusion method of optical image, including Following steps：

Step 1, using mathematical morphology filter method cloud data is filtered process：

Morphologic filtering removes cancelling noise as a kind of wave filter being mutated thought based on the local discrepancy in elevation from bottom to top, purpose Point, redundancy and elevation rough error；Corrosion and dilation operation are the bases that morphological images are processed, for " minimizing " The size of character shape in (corrosion) or " increase " (expansion) image；Corrosion based on gray level image and dilation operation It is the minimum after selecting image pixel value phase separation in the neighborhood of structural element definition or max pixel value；

If LiDAR point cloud data observation value sequence is p (x, y, z), then at x and y, the expansive working to height z is defined as：

In formula, (x_p, y_p, z_p) representing point in p neighborhood of a point window w, w is also called size of structure element；Neighborhood window Can be one-dimensional straight line it is also possible to make rectangle or the other shapes of two dimension；The result of dilation operation is in neighborhood window Highest elevation value；

Corrosion, as the dual operator expanding, is defined as：

Expansion and corrosion are combined, obtain opening operation and the closed operation that can be directly used for LiDAR filtering, during opening operation Data wire is corroded, then is expanded, and closed operation contrast.

Using linear structure element, a certain vertical section of LiDAR point cloud data is carried out after opening operation filtering, if trees target It is smaller in size than construction unit, then can be removed after erosion operation, building then obtains weight in dilation operation Build.

Step 2, it is respectively adopted self adaptation and supports weight dense stereo algorithm and Delaunay Triangulation algorithm to determine optics Image and the depth map of cloud data, specially：

For cloud data, the three-dimensional point cloud discrete data P ' after morphologic filtering, first pass through Delaunay triangle Change the gridding of algorithm construction LiDAR point cloud data, then mapped net according to the color that height value z carries out gray level The LiDAR data formatted is converted into depth map；

For optical image I, in order to obtain its depth map, we support weight from based on related in color and self adaptation Dense stereo coupling ask for optical image I and be adjacent the optimum Stereo matching side of image I ' (there is overlapping region) Method.The method combines color similarity, and in Euclidean distance similarity and color, associated similarity is determining in match window The weight size of pixel；Simultaneously in order to eliminate the different impact to images match result of illumination, match point is first carried out Carry out Matching power flow relation calculating again after rank conversion.Finally, three step optimizations are carried out to the initial parallax figure calculating, pick Remove and the different parallax mistakes causing such as blocked, repeat by image, thus obtaining final parallax result.It is specially：

Step 2-1, for image I, its RGB triple channel color vector be I_R(x, y), I_G(x, y), I_B(x, y), wherein X, y represent the coordinate figure on image slices vegetarian refreshments ranks direction, and related rin vector representation in color is：

Rin=[rce₁, rce₂, rce₃]

In formula, rce is defined as：

rce₁=I_R(x, y)-I_G(x, y)

rce₂=I_G(x, y)-I_B(x, y)

rce₃=I_R(x, y)-I_R(x, y)

Support weight between pixel p and q to be matched is expressed as：

W (p, q)=f (Δ c_pq, Δ d_pq, Δ r_pq)

=f (Δ c_pq)·f(Δd_pq)·f(Δr_pq)

Wherein, Δ c_pq, Δ d_pqWith Δ r_pqRepresent that color distortion between pixel p and q to be matched, Euclidean distance are poor respectively Relevant difference in different and color, f () represents difference intensity, and it is expressed as：

In formula, τ_c, τ_d, τ_rIt is respectively fixed constant, in present embodiment, take τ_c=5, τ_d=17.5, τ_r=5.

So, the support weight between pixel p and q to be matched is expressed as：

Produce corresponding matching error in order to eliminate two width images because illumination is different, the point of coupling is first carried out rank change Change, then carry out Matching power flow calculating again, specific rank conversion expression formula is as follows：

Wherein I_pFor the gray value of current matching pixel, I_qFor the gray value of pixel in window, τ₁, τ₂For rank conversion Class condition, in practice for fixed constant, takes τ in present embodiment₁=9, τ₂=2.

Optical image left view I and optical image right view I ' match point p andCorresponding point in match window For q andIts Matching power flow value isIt is expressed as：

In formula, N represents match window size, d_maxFor the maximum disparity value of matched pixel point,Represent left and right figure As in q andThe matching value at place：

Step 2-2, according to the victor is a king (WTA) algorithm, when left images scan for coupling, using polar curve about Restraint, the parallax of final p point is：

Step 2-3, left and right consistency check and medium filtering are carried out to initial parallax, estimated with reference to Epipolar geometry by disparity map Meter obtains the depth map D of optical image_I.

Step 3, extract the line feature of depth map by Hough algorithm, using line feature two wire clamp angle and line length ratio as Similarity measure is slightly mated, and comprises the following steps that：

Step 3-1, rim detection is carried out using Canny operator, the edge image generating is carried using Hough transform Cut-off line：

(1) obtain linear feature, adjust Hough transform parameter, obtain linear feature；

(2) adjacent segments AB, CD head and the tail end-point distances are less than first threshold Th₁Line segment to be connected to become same straight Line AD, by straight line extended line in Second Threshold Th₂Intersecting line segment l in pixel₁、l₂Connect intersecting；

(3) line segment extracting is carried out line length to compare, if line length d is less than the 3rd threshold value Th₃, then this line segment be ignored；

Step 3-2, in the line feature that cloud data depth map and visible images depth map extract, with two wire clamp angles With line length ratio as similarity measure；It is specially：

First line length ratio r is compared, if through comparing left view line to be matched feature line length ratio r₁Treat with right view Matched line feature line length ratio r₂Difference d_rIn the 4th threshold value T_rIn the range of, then prestore and carry out the angle changing rate of next step, If through comparing left view line to be matched feature angle α₁With right view line to be matched feature angle α₂Difference d_αIn the 5th threshold Value T_αIn the range of, then for the line feature of a pair of coupling；

If with line feature F₁The feature of coupling has multiple, then to line feature F₁The line feature of neighborhood and matched line are special Line feature near levying is carried out two wire clamp angles and is compared with line length ratio, draws the line feature of Optimum Matching, finally gives many In the best match line pair of 3, the end points choosing each pair line pair is asked for similarity transformation parameter, is realized the rough registration of image.

Step 4, by two step RANSAC algorithms reject error matching points pair, obtain camera position parameter estimation, specifically Step is as follows：

Step 4-1, according to spatial neighbor degree it is assumed that matching characteristic to being divided into different groups, that is, piece image is divided It is slit into different windows, window size is set to 4s*4s pixel, s is positive integer；Window is from the upper left corner of image from a left side Start from top to bottom to slide to the right side, slip unit is s, thus carrying out next step segmentation, so that the coupling in each group is special Levy and intersection be might have to the matching double points in quantity each group in a suitable threshold range；

Step 4-2, the matching characteristic that each is organized calculate to carrying out a step RANSAC algorithm, reject point not in the know：Repeatedly During generation, two matching characteristics calculate corresponding homography matrix to being sampled simultaneously, if other matching characteristic is to symbol Closing this matrix is considered as then intra-office point, and otherwise for point not in the know, the matching characteristic electing intra-office point as mates spy to for basic Levy；

Step 4-3, one again the basic matching characteristic in different groups is carried out towards an overall step RANSAC algorithm Calculate；The basic matching characteristic of each two group is sampled the corresponding homography matrix of calculating, if other basic matching characteristic Meeting this matrix is considered as then intra-office point, and otherwise for point not in the know, the fraction equation below of each iteration calculates. Assume to have point in j group, interior number of kth group is n_k, iterate to calculate every time must be divided into：

Represent the weight of kth group, when interior point concentrates on several groups, weight is larger；When interior point is distributed in institute When having group, weight will decline；I=1～j, expression group number；

It is referred to as the matching characteristic pair of intra-office point after step 4-4, all calculating by step 4-3 of reservation, be referred to as point not in the know Matching characteristic to then again being calculated again using fraction highest homography matrix, to retain more features pair；Traversal After end, all of intra-office point is the matching characteristic pair of our needs.

Because, in Feature Points Matching, mapping model is the feature in the characteristic point a plane to an other plane The projective rejection of point, reacts for projection matrix；In the match point eliminating mistake through two steps RANSAC, i.e. point not in the know Afterwards, intra-office point is used to calculate mapping matrix.

Step 5, be based on above two step RANSAC algorithms, we will reject mistake matching double points, obtain high accuracy Phase position parameter estimation, carry out the re-projection between cloud data and optical image, complete cloud data and optics shadow The color and vein mapping of picture, the 3-D view after being merged.

With reference to specific embodiment, the invention will be further described.

Embodiment

The present embodiment adopts airborne LiDAR point cloud data and optical image, mainly covers portion centers area, comprises to lead Want building, short jungle etc., have chosen part airborne LiDAR point cloud data, with optical image, there is overlapping region Data, is processed using improved line feature registration algorithm.

Fig. 2 (a) is cloud data, and Fig. 2 (b) is visible images, and Fig. 2 (c) is the 3-D view after merging, Image be can be seen that according to the image after merging and achieve the three-dimensional reconstruction to city, make point cloud chart picture have visible ray The color information of image.

The cloud data based on line feature for the present invention and the autoregistration fusion method of optical image, without GPS/INS etc. In the case of priori, algorithm high degree of automation, robustness is good, and amount of calculation is little, just has higher registration simultaneously Really rate.

Claims (6)

1. a kind of cloud data based on line feature and the autoregistration fusion method of optical image are it is characterised in that wrap Include following steps：

Step 1, is filtered to cloud data processing；

Step 2, is respectively adopted self adaptation and supports weight dense stereo algorithm and Delaunay Triangulation algorithm to determine optics Image and the depth map of cloud data；

Step 3, extracts the line feature of depth map by Hough algorithm, using line feature two wire clamp angle and line length ratio as Similarity measure is slightly mated；

Step 4, rejects error matching points pair by two step RANSAC algorithms, obtains camera position parameter estimation；

Step 5, carries out the color and vein mapping of cloud data and optical image, the 3-D view after being merged.

2. the autoregistration fusion method of the cloud data based on line feature according to claim 1 and optical image, It is characterized in that, in step 1, cloud data is filtered processing using Mathematical Morphology Method.

3. the autoregistration fusion method of the cloud data based on line feature according to claim 1 and optical image, It is characterized in that, the depth map of cloud data, detailed process is determined in step 2 using Delaunay Triangulation algorithm For：

First pass through the gridding that Delaunay Triangulation Algorithm constructs cloud data, then gray scale is carried out according to height value z The color mapping of level, the cloud data after gridding is converted into depth map D_P.

4. the autoregistration fusion method of the cloud data based on line feature according to claim 1 and optical image, It is characterized in that, determine the depth map of optical image in step 2, support power using based on related in color and self adaptation The weight dense stereo matching algorithm Stereo matching that to determine optical image I and the adjacent image I ' with overlapping region between optimum； Detailed process is：

Step 2-1, for optical imagery I, its RGB triple channel color vector is I_R(x, y), I_G(x, y), I_B(x, y), Wherein x, y represent the coordinate figure on image slices vegetarian refreshments ranks direction；Related rin vector representation in color is：

Rin=[rce₁, rce₂, rce₃]

In formula：

rce₁=I_R(x, y)-I_G(x, y)

rce₂=I_G(x, y)-I_B(x, y)

rce₃=I_B(x, y)-I_R(x, y)

Support weight between pixel p and q to be matched is expressed as：

W (p, q)=f (Δ c_pq, Δ d_pq, Δ r_pq)

=f (Δ c_pq)·f(Δd_pq)·f(Δr_pq）

Wherein, Δ c_pq, Δ d_pqWith Δ r_pqRepresent that color distortion between pixel p and q to be matched, Euclidean distance are poor respectively Relevant difference in different and color, f () represents difference intensity, and it is expressed as：

$f\left({\mathrm{\Δc}}_{pq}\right)=\exp (-\frac{{\mathrm{\Δc}}_{pq}}{{\mathrm{\τ}}_c})$

$f\left({\mathrm{\Δd}}_{pq}\right)=\exp (-\frac{{\mathrm{\Δd}}_{pq}}{{\mathrm{\τ}}_d})$

$f\left({\mathrm{\Δr}}_{pq}\right)=\exp (-\frac{{\mathrm{\Δr}}_{pq}}{{\mathrm{\τ}}_r})$

In formula, τ_c, τ_d, τ_rIt is respectively fixed constant；

So, the support weight between pixel p and q to be matched is expressed as：

$\begin{array}{c}w(p,q)=f({\mathrm{\Δc}}_{pq},{\mathrm{\Δd}}_{pq},{\mathrm{\Δr}}_{pq})\\ =f\left({\mathrm{\Δc}}_{pq}\right)\·f\left({\mathrm{\Δd}}_{pq}\right)\·f\left({\mathrm{\Δr}}_{pq}\right)\\ =\exp (-(\frac{{\mathrm{\Δc}}_{pq}}{{\mathrm{\τ}}_c}+\frac{{\mathrm{\Δd}}_{pq}}{{\mathrm{\τ}}_d}+\frac{{\mathrm{\Δr}}_{pq}}{{\mathrm{\τ}}_r}\left)\right)\end{array}$

The point of coupling is carried out rank conversion, rank conversion expression formula is as follows：

$R\left(q\right)=\left\{\begin{array}{cc}-2& I_p-I_q<-{\mathrm{\&tau;}}_1\\ -1& -{\mathrm{\&tau;}}_1\&le;I_p-I_q\&le;-{\mathrm{\&tau;}}_2\\ 0& -{\mathrm{\&tau;}}_2<I_p-I_q\&le;{\mathrm{\&tau;}}_2\\ 1& {\mathrm{\&tau;}}_2<I_p-I_q\&le;{\mathrm{\&tau;}}_1\\ 2& I_p-I_q>{\mathrm{\&tau;}}_1\end{array}\right.$

Wherein I_pGray value for p, I_qGray value for q, τ₁, τ₂Convert class condition for rank；

Optical image left view I and optical image right view I ' match point be p andCorresponding point in match window be q andIts Matching power flow value isIt is expressed as：

$E(p,{\stackrel{\&OverBar;}{p}}_d)=\frac{{\mathrm{\&Sigma;}}_{q,{\stackrel{\&OverBar;}{q}}_d\&Subset;N}w(p,q)w({\stackrel{\&OverBar;}{p}}_d,{\stackrel{\&OverBar;}{q}}_d)e(q,{\stackrel{\&OverBar;}{q}}_d)}{{\mathrm{\&Sigma;}}_{q,{\stackrel{\&OverBar;}{q}}_d\&Subset;N}w(p,q)w({\stackrel{\&OverBar;}{p}}_d,{\stackrel{\&OverBar;}{q}}_d)}$

In formula, N represents match window size, d_maxFor the maximum disparity value of matched pixel point,Represent left and right figure As in q andThe matching value at place；

Step 2-2, according to the victor is a king algorithm, when left images scan for coupling, using epipolar-line constraint, final p Point parallax be：

$d_p=\arg \max E(p,{\stackrel{\&OverBar;}{p}}_d)$

Step 2-3, carries out left and right consistency check and medium filtering to initial parallax, is estimated with reference to Epipolar geometry by disparity map Meter obtains the depth map D of optical image_I.

5. the autoregistration fusion method of the cloud data based on line feature according to claim 1 and optical image, It is characterized in that, step 3 is specially：

Step 3-1, rim detection is carried out using Canny operator, the edge image generating is carried using Hough transform Cut-off line：

(1) adjust Hough transform parameter, obtain linear feature；

(2) adjacent segments AB, CD head and the tail end-point distances are less than first threshold Th₁Line segment to be connected to become same straight Line AD, by straight line extended line in Second Threshold Th₂Intersecting line segment l in pixel₁、l₂Connect intersecting；

(3) line segment extracting is carried out line length to compare, if line length d is less than the 3rd threshold value Th₃, then this line segment be ignored；

Step 3-2, in the line feature that cloud data depth map and visible images depth map extract, with two wire clamp angles With line length ratio as similarity measure, specially：

First line length ratio r is compared, if through comparing left view line to be matched feature line length ratio r₁Treat with right view Matched line feature line length ratio r₂Difference d_rIn the 4th threshold value T_rIn the range of, then prestore and carry out the angle changing rate of next step, If through comparing left view line to be matched feature angle α₁With right view line to be matched feature angle α₂Difference d_αIn the 5th threshold Value T_αIn the range of, then for the line feature of a pair of coupling；

If with line feature F₁The feature of coupling has multiple, then to line feature F₁The line feature of neighborhood and matched line are special Line feature near levying is carried out two wire clamp angles and is compared with line length ratio, draws the line feature of Optimum Matching, finally gives many In the best match line pair of 3, the end points choosing each pair line pair is asked for similarity transformation parameter, is realized the rough registration of image.

6. the autoregistration fusion method of the cloud data based on line feature according to claim 1 and optical image, It is characterized in that, step 4 is specially：

Step 4-1, according to spatial neighbor degree it is assumed that matching characteristic to being divided into different groups, that is, piece image is divided It is slit into different windows, window size is set to 4s*4s pixel, s is positive integer；Window is from the upper left corner of image from a left side Start from top to bottom to slide to the right side, slip unit is s, thus carrying out next step segmentation；

Step 4-2, the matching characteristic that each is organized calculate to carrying out a step RANSAC algorithm, reject point not in the know：Repeatedly During generation, two matching characteristics calculate corresponding homography matrix to being sampled simultaneously, if other matching characteristic is to symbol Closing this matrix is considered as then intra-office point, otherwise for point not in the know；The matching characteristic electing intra-office point as mates spy to for basic Levy；

Step 4-3, again the basic matching characteristic in different groups is calculated using a step RANSAC algorithm；Every two Individual group of basic matching characteristic is sampled the corresponding homography matrix of calculating, if other basic matching characteristic meets this square Battle array is considered as then intra-office point, otherwise for point not in the know；The fraction equation below of each iteration obtains：

Assume to have point in j group, interior number of kth group is n_k, iterate to calculate every time must be divided into：

$G={\mathrm{\&Sigma;}}_{i=1}^j{n}_k\&CenterDot;\sqrt{n_k}$

For the weight of kth group, i=1～j, expression group number；

It is considered the matching characteristic pair of intra-office point in step 4-4, reservation step 4-3, be referred to as the matching characteristic of point not in the know To then being calculated again using fraction highest homography matrix；After traversal terminates, all of intra-office point is our needs Matching characteristic pair；

In Feature Points Matching, mapping model is the characteristic point in the characteristic point a plane to another one plane Projective rejection, mapping matrix is determined by intra-office point, obtains phase position parameter estimation.