CN104915988A - Photogrammetric dense point cloud generation method - Google Patents

Abstract

The invention discloses a photogrammetric dense point cloud generation method comprising the steps that exterior orientation elements of left and right images are obtained by utilizing a photogrammetric backward intersection method; global matching is performed on the left and right images and the corresponding relation between the left and right image points is established so that all the image points are enabled to possess matching points, parallax graphs used for describing the left and right images are acquired, and the coordinate values of the homologous image points of the left and right images in an image space coordinate system are calculated on the basis of the parallax graphs of the left and right images; and three-dimensional position information of all the corresponding ground points are recovered by a forward intersection method according to the exterior orientation elements of left and right images and the coordinates of the homologous image points of the left and right images in the image space coordinate system so that dense three-dimensional point clouds are generated. With application of the disclosed method, a matching problem of image no-texture area image points is solved, LiDAR point cloud data holes can be compensated and thus technical basis is established for construction of a digital city and a smart city.

Description

A kind of photogrammetric point of density cloud generation method

Technical field

The present invention relates to technical field of mapping, particularly relate to a kind of photogrammetric point of density cloud generation method.

Background technology

Digital city and smart city are built, and propose requirements at the higher level for the technology how rapidly and efficiently obtaining City surface three-dimensional spatial information.Airborne LIDAR (laser radar) is a kind of initiatively earth observation systems, is integrated with the technology such as laser ranging, inertia measurement (IMU), GNSS Differential positioning.The laser pulse energy of LIDAR sensor emission partly penetrates the woods and blocks, and directly obtains the cloud data of high-precision three-dimensional earth's surface landform.Have that automaticity is high, little by weather effect, data are with short production cycle, precision high.Airborne LIDAR data, after Correlation method for data processing, can generate high-precision digital terrain model, contour map.Airborne LiDAR point cloud data also for the target classification such as City Building, road and extraction, but utilize on-board LiDAR data to carry out the automatic identification and extraction of City Building, larger with extraction comparison difficulty such as the buildingss based on remote sensing image.Most of LiDAR point cloud Data Management Analysis, mainly concentrates on the geometric characteristic analysis (as curved surface and principal curvatures analysis, TIN analysis etc.) of cloud data.Modern LiDAR sensor research is focusing more on multiecho and the analysis of the multi-sensor information fusion such as texture information, strength information.Such as based on the classification tree algorithm of multiecho data, be applied to target classification and extractions such as City Building, road and trees.There is data void holes in the cloud data measured due to LiDAR, needs to utilize other technologies to carry out compensation data.

Digital photogrammetry comprises that satellite three line scanner is photogrammetric, aviation three line scanner is photogrammetric, unmanned plane is photogrammetric, the mainly multidisciplinary technical method such as appliance computer, digitized video process, Image Matching, pattern-recognition.Traditional photography is measured Image Matching and is mostly feature-based matching method, is to utilize automatic corresponding image points coupling and aerotriangulation to resolve and obtain object dimensional locus.When image existing weak texture area, traditional algorithm effectively can not extract image feature point and carry out reliable Feature Matching, thus cannot resolve object dimensional locus; That is the method for feature based coupling, cannot solve cloud data cavity problem.

Summary of the invention

The object of this invention is to provide a kind of photogrammetric point of density cloud generation method, solve the matching problem of image without the picture point of texture region, LiDAR point cloud data void holes can be realized and compensate, for digital city and smart city construction lay the foundation.

The object of the invention is to be achieved through the following technical solutions:

A kind of photogrammetric point of density cloud generation method, comprising:

Photogrammetric resect is utilized to ask for the elements of exterior orientation of left and right image;

Corresponding relation between the picture point of left and right is set up by carrying out global registration to left and right image, all picture points are made to have its match point, thus the disparity map obtained for describing left and right image, and calculate the coordinate of left and right image corresponding image points in image space coordinate system based on the disparity map of left and right image;

According to elements of exterior orientation and the coordinate of left and right image corresponding image points in image space coordinate system of described left and right image, and calculate the three-dimensional coordinate of all ground points in earth coordinates by forward intersection, thus generate intensive three-dimensional point cloud.

Further, the described elements of exterior orientation utilizing photogrammetric resect to ask for left and right image comprises:

Based on collinearity equation, set up error equation:

$\left.\begin{array}{c}{V}_x=a_{11}{\mathrm{\ΔX}}_S+a_{12}{\mathrm{\ΔY}}_S+a_{13}{\mathrm{\ΔZ}}_S+a_{14}\mathrm{\Δ}\mathrm{\φ}+a_{15}\mathrm{\Δ}\mathrm{\ω}+a_{16}\mathrm{\Δ}\mathrm{\κ}-l_x\\ V_y=a_{21}{\mathrm{\ΔX}}_S+a_{22}{\mathrm{\ΔY}}_S+a_{23}{\mathrm{\ΔZ}}_S+a_{24}\mathrm{\Δ}\mathrm{\φ}+a_{25}\mathrm{\Δ}\mathrm{\ω}+a_{26}\mathrm{\Δ}\mathrm{\κ}-l_y\end{array}\right\};$

Wherein, Δ X _s, Δ Y _s, Δ Z _s, Δ φ, Δ ω, Δ κ are unknown number to be asked, and it represents image elements of exterior orientation (X _s, Y _s, Z _s, φ, ω, κ) correction, a ₁₁, a ₁₂..., a ₁₆, a ₂₁, a ₂₂..., a ₂₆for error equation coefficient, l _x, l _yfor constant term, V _x, V _yfor error equation correction;

The matrix form of this error equation is:

V＝AX-l；

Wherein, V=[V _xv _y] ^t, l=[l _xl _x] ^t, $A=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& a_{13}& a_{14}& a_{15}& a_{16}\\ a_{21}& a_{22}& a_{23}& a_{24}& a_{25}& a_{26}\end{array}\right],$ X＝[ΔX _SΔY _SΔZ _SΔφ Δω Δκ] ^T；

The normal equation of above-mentioned matrix is:

A ^TPAX＝A ^TPL；

Wherein, P is the weight matrix of observed reading, then the expression formula of unknown number to be asked is:

X＝(A ^TA) ^-1A ^TL；

Again iterative computation is carried out to the expression formula of unknown number above-mentioned to be asked, if the last difference calculating unknown number to be asked for twice is less than predetermined value, then calculate stopping.

Further, describedly set up corresponding relation between the picture point of left and right by carrying out global registration to left and right image, all picture points are made to have its match point, thus obtain disparity map for describing left and right image, and calculate the left and right coordinate figure of image corresponding image points in image space coordinate system based on the disparity map of left and right image and comprise:

First, set up the energy function of left and right image parallax, it is expressed as:

$E\left(d\right)=\underset{(p,q)\∈N}{\Σ}V(d_p,d_q)+\underset{p\∈P}{\Σ}{D}_p\left(d_p\right);$

Wherein, d represents that the parallax of entire image distributes; N represents four neighborhood point sets of all picture points in image; d _prepresent the parallax value that pixel p distributes; P represents the set of picture point in image;

Non-set of metadata of similar data item D _p(d _p): represent that picture point p parallax is d _ptime non-similarity tolerance, calculated by following formula: D _p(d _p)=λ min (| I _l(x, y)-I _r((x+d _p), y) |, T), in formula, λ is regulatory factor, and T represents truncation error, I _l(x, y), I _r(x+dx, y+dy) is respectively image intensity corresponding to left and right image picture point;

Level and smooth item V (d _p, d _q): represent that two adjacent image points p and q distribute parallax d _pand d _qtime the discontinuous punishment amount of parallax, determined by following formula:

In above formula, a ₁, a ₂with a ₃for preset value, meet a ₃> a ₂> a ₁> 0; I _pwith I _qrepresent the intensity of picture point p and q respectively;

Secondly, utilize belief propagation optimization to calculate and this energy function is minimized, distribute the final parallax obtaining left and right image according to the parallax in minimized energy function;

Finally, the coordinate figure of left and right image corresponding image points in image space coordinate system is calculated according to the disparity map of left and right image.

Further, describedly utilize belief propagation algorithm that this energy function is minimized to comprise:

Information updating, its expression formula is:

$m_{p\→q}^t\left(d_q\right)=\underset{d_p\∈\mathrm{\Ω}}{min}\[V(d_p,d_q)+D_p\left(d_p\right)+\underset{r\∈N\left(p\right)\backslash q}{\Σ}{m}_{r\→p}^{t-1}\left(d_p\right)\];$

In formula, Ω represents disparity search scope, when being the t time iteration, the message that picture point p transmits to picture point q, N (p) q represent that picture point p removes three neighborhood nodes outside picture point q, comprise x ₁, x ₂, x ₃; when representing the t-1 time iteration, the message that r point transmits to picture point q;

$\underset{r\∈N\left(p\right)\backslash q}{\Σ}{m}_{r\→p}^{t-1}\left(d_p\right)=m_{x_1\→p}^{t-1}+m_{x_2\→p}^{t-1}+m_{x_3\→p}^{t-1};$

Confidence calculations, its expression formula is:

$b_p\left(d_p\right)=D_p\left(d_p\right)+\underset{r\∈N\left(p\right)}{\Σ}{m}_{r\→p}^T\left(d_p\right);$

B _p(d _p) represent that the parallax value that picture point p distributes is d _ptime degree of confidence, N (p) represents the four neighborhood nodes of picture point p, comprises q, x ₁, x ₂, x ₃; Wherein:

$\underset{r\∈N\left(p\right)}{\Σ}{m}_{r\→p}^T\left(d_p\right)=m_{x_1\→p}^T+m_{x_2\→p}^T+m_{x_3\→p}^T+m_{q\→p}^T;$

Minimizing degree of confidence makes this energy function minimize, and the expression formula minimizing degree of confidence is:

$d_p^{*}=\underset{d_p\∈\mathrm{\Ω}}{\arg \min }{b}_p\left(d_p\right).$

Further, the described elements of exterior orientation according to described left and right image and the coordinate of left and right image corresponding image points in image space coordinate system, and calculate the three-dimensional coordinate of all ground points in earth coordinates by forward intersection, thus generate intensive three-dimensional point cloud and comprise:

From collinearity equation:

$\left\{\begin{array}{c}\frac{X_A-X_{S1}}{X_1}=\frac{Y_A-Y_{S1}}{Y}=\frac{Z_A-Z_{S1}}{Z_1}=N_1\\ \frac{X_A-X_{S2}}{X_2}=\frac{Y_A-Y_{S2}}{Y_2}=\frac{Z_A-Z_{S2}}{Z_2}=N_2\end{array}\right.;$

In formula, N ₁, N ₂be respectively left and right image picture point projection coefficient; (X _a, Y _a, Z _a) be the three-dimensional coordinate of ground point in earth coordinates; The elements of exterior orientation of left and right image is designated as (X respectively _s1, Y _s1, Z _s1, φ ₁, ω ₁, κ ₁), (X _s2, Y _s2, Z _s2, φ ₂, ω ₂, κ ₂) _,(Xs ₁, Ys ₁, Zs ₁) and (Xs ₂, Ys ₂, Zs ₂) be respectively image photographic center, left and right, (φ ₁, ω ₁, κ ₁) and (φ ₂, ω ₂, κ ₂) be respectively the attitude angle of left and right image; (X ₁, Y ₁, Z ₁) and (X ₂, Y ₂, Z ₂) be respectively the coordinate of left and right image corresponding image points in the auxiliary coordinates of image space, be expressed as:

$\begin{array}{cc}\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]=R_1\·\left[\begin{array}{c}{x_1}^{\′}\\ y_1\\ -f\end{array}\right];& \left[\begin{array}{c}{X}_2\\ Y_2\\ Z_2\end{array}\right]=R_2\·\left[\begin{array}{c}{x_2}^{\′}\\ y_2\\ -f\end{array}\right]\end{array};$

Wherein, the coordinate of left and right image corresponding image points in image space coordinate system, is designated as (x respectively ₁', y ₁,-f), (x ₂', y ₂,-f); R ₁and R ₂be respectively the attitude angle (φ with left image ₁, ω ₁, κ ₁), the attitude angle (φ of right image ₂, ω ₂, κ ₂) relevant rotation matrix;

Calculate the three-dimensional coordinate (X of ground point in earth coordinates _a, Y _a, Z _a), its formula is as follows:

$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\end{array}\right]=\left[\begin{array}{c}{X}_{S1}\\ Y_{S1}\\ Z_{S1}\end{array}\right]+\left[\begin{array}{c}{N}_1{X}_1\\ N_1{Y}_1\\ N_1{Z}_1\end{array}\right]=\left[\begin{array}{c}{X}_{S2}\\ Y_{S2}\\ Z_{S2}\end{array}\right]+\left[\begin{array}{c}{N}_2{X}_2\\ N_2{Y}_2\\ N_2{Z}_2\end{array}\right];$

In formula:

$\begin{array}{cc}\left\{\begin{array}{c}{N}_1=\frac{B_X{Z}_2-B_Z{X}_2}{X_1{Z}_2-X_2{Z}_1}\\ N_2=\frac{B_X{Z}_1-B_Z{X}_1}{X_1{Z}_2-X_2{Z}_1}\end{array}\right.,& \left\{\begin{array}{c}{B}_X=X_{S2}-X_{S1}\\ B_Y=Y_{S2}-Y_{S1}\\ B_Z=Z_{S2}-Z_{S1}\end{array}\right.\end{array}$

The three-dimensional coordinate of all ground points in earth coordinates can be obtained by calculating above, thus generate intensive three-dimensional point cloud according to its three-dimensional coordinate.

As seen from the above technical solution provided by the invention, by asking for left and right image elements of exterior orientation according to photogrammetric resection, then corresponding relation is set up by the coupling by pixel, acquisition can describe left and right image disparity map, finally carry out forward intersection and recover three dimensional local information, generate intensive three-dimensional point cloud; Not only solve the Pixel matching problem of image without texture region, also can be used for LiDAR point cloud data void holes and compensate, for digital city and smart city construction lay the foundation.

Accompanying drawing explanation

In order to be illustrated more clearly in the technical scheme of the embodiment of the present invention, below the accompanying drawing used required in describing embodiment is briefly described, apparently, accompanying drawing in the following describes is only some embodiments of the present invention, for those of ordinary skill in the art, under the prerequisite not paying creative work, other accompanying drawings can also be obtained according to these accompanying drawings.

The process flow diagram of a kind of photogrammetric point of density cloud generation method that Fig. 1 provides for the embodiment of the present invention;

The schematic diagram of the pass-along message that Fig. 2 provides for the embodiment of the present invention;

The schematic diagram of the degree of confidence transmission that Fig. 3 a provides for the embodiment of the present invention;

The schematic diagram of the confidence calculations that Fig. 3 b provides for the embodiment of the present invention.

Embodiment

Below in conjunction with the accompanying drawing in the embodiment of the present invention, be clearly and completely described the technical scheme in the embodiment of the present invention, obviously, described embodiment is only the present invention's part embodiment, instead of whole embodiments.Based on embodiments of the invention, those of ordinary skill in the art, not making the every other embodiment obtained under creative work prerequisite, belong to protection scope of the present invention.

Embodiment

The process flow diagram of a kind of photogrammetric point of density cloud generation method that Fig. 1 provides for the embodiment of the present invention.As shown in Figure 1, the method mainly comprises the steps:

Step 11, photogrammetric resect is utilized to ask for the elements of exterior orientation of left and right image.

In this step, collinearity equation is expressed as:

$\left.\begin{array}{c}x=-f\frac{a_1(X_A-X_S)+b_1(Y_A-Y_S)+c_1(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)}\\ y=-f\frac{a_2(X_A-X_S)+b_2(Y_A-Y_S)+c_2(Z_A-Z_S)}{a_3(X_A-X_S)+b_3(Y_A-Y_S)+c_3(Z_A-Z_S)}\end{array}\right\}$

In above-mentioned collinearity equation: (x, y) is image picpointed coordinate, and f is camera focal length, (X _a, Y _a, Z _a), (X _s, Y _s, Z _s) be respectively ground point, the photo centre coordinate in earth coordinates, (φ, ω, κ) is the attitude angle of image.In formula:

$R(\mathrm{\φ},\mathrm{\ω},\mathrm{\κ})=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$

Be the rotation matrix relevant with image attitude angle (φ, ω, κ), also claim (XS, YS, ZS, φ, ω, κ) to be image elements of exterior orientation.

Taylor series expansion is carried out to collinearity equation and is taken to single order item, can error equation be obtained:

$\left.\begin{array}{c}{V}_x=a_{11}{\mathrm{\ΔX}}_S+a_{12}{\mathrm{\ΔY}}_S+a_{13}{\mathrm{\ΔZ}}_S+a_{14}\mathrm{\Δ}\mathrm{\φ}+a_{15}\mathrm{\Δ}\mathrm{\ω}+a_{16}\mathrm{\Δ}\mathrm{\κ}-l_x\\ V_y=a_{21}{\mathrm{\ΔX}}_S+a_{22}{\mathrm{\ΔY}}_S+a_{23}{\mathrm{\ΔZ}}_S+a_{24}\mathrm{\Δ}\mathrm{\φ}+a_{25}\mathrm{\Δ}\mathrm{\ω}+a_{26}\mathrm{\Δ}\mathrm{\κ}-l_y\end{array}\right\}$

In formula, Δ X _s, Δ Y _s, Δ Z _s, Δ φ, Δ ω, Δ κ are unknown number to be asked, and it represents image elements of exterior orientation (X _s, Y _s, Z _s, φ, ω, κ) correction, a ₁₁, a ₁₂..., a ₁₆, a ₂₁, a ₂₂..., a ₂₆for error equation coefficient, l _x, l _yfor constant term, V _x, V _yfor error equation correction.

Error equation is written as matrix form: V=AX-l;

Wherein: V=[V _xv _y] ^t, l=[l _xl _x] ^t, $A=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& a_{13}& a_{14}& a_{15}& a_{16}\\ a_{21}& a_{22}& a_{23}& a_{24}& a_{25}& a_{26}\end{array}\right],$ X＝[ΔX _SΔY _SΔZ _SΔφ Δω Δκ] ^T；

The normal equation of above-mentioned matrix is: A ^tpAX=A ^tpL;

Wherein, P is the weight matrix of observed reading, the measurement accuracy of reaction observed reading, then the expression formula of unknown number to be asked is:

X＝(A ^TA) ^-1A ^TL；

Again iterative computation is carried out to the expression formula of unknown number above-mentioned to be asked, if the last difference calculating unknown number to be asked for twice is less than predetermined value, then calculate stopping;

Calculate Δ X _s, Δ Y _s, Δ Z _s, Δ φ, Δ ω, after Δ κ, then can obtain corresponding image elements of exterior orientation (X _s, Y _s, Z _s, φ, ω, κ).

Step 12, set up corresponding relation between the picture point of left and right by carrying out global registration to left and right image, all picture points are made to have its match point, thus the disparity map obtained for describing left and right image, and calculate the coordinate figure of left and right image corresponding image points in image space coordinate system based on the disparity map of left and right image.

In the embodiment of the present invention, the disparity map of left and right image is dense disparity map, and its concrete steps are as follows:

First, set up the energy function of left and right image parallax, it is expressed as:

$E\left(d\right)=\underset{(p,q)\∈N}{\Σ}V(d_p,d_q)+\underset{p\∈P}{\Σ}{D}_p\left(d_p\right);$

Wherein, d represents that the parallax of entire image distributes, and N represents four neighborhood point sets of all picture points in image, d _prepresent the parallax value that picture point p distributes, P represents the set of picture point in image.

Non-set of metadata of similar data item D _p(d _p) represent that picture point p parallax is d _ptime non-similarity tolerance, calculated by following formula:

D _p(d _p)＝λmin(|I _l(x,y)-I _r(x+dx,y+dy)|,T)，

In formula, λ is regulatory factor, and T represents truncation error.I _l(x, y), I _r(x+dx, y+dy) is respectively image intensity corresponding to left and right image picture point.

Level and smooth item V (d _p, d _q) represent that two adjacent image points p and q distribute parallax d _pand d _qtime the discontinuous punishment amount of parallax, determined by following formula:

That is, if the parallax of picture point p and q is equal, so the discontinuous punishment amount of parallax is 0; If the disparity difference of picture point p and q is within a pixel, to a less punishment amount a ₁; And if difference is outside a pixel, may be now cause parallax discontinuous owing to blocking, also should give less punishment amount, but the inaccurate phenomenon in edge that Iamge Segmentation draws may be there is here, thus need again to judge two picture points intensity difference (| I _p-I _q|), if intensity difference is greater than threshold value δ, the actual boundary caused is blocked in explanation, and parallax is discontinuous, should give the punishment amount a that is less ₂and if intensity difference is less than threshold value δ, then illustrate parallax discontinuous be because the disparity plane template of distributing is not optimum template, a larger punishment amount a should be given ₃.Wherein, a ₁, a ₂with a ₃for preset value, meet a ₃> a ₂> a ₁> 0; I _pwith I _qrepresent the intensity of picture point p and q respectively.

Now, initial parallax figure can be obtained based on core line: calculate basis matrix F; Ask and the limit e1 of epipolar lines intersect, e2; Select transformation matrix U2, limit e2 is mapped to infinite distance limit e; According to formula solve and obtain U2; Correction is completed according to U1 and U2 resampling; The x coordinate difference correcting secondary homonym picture point is initial parallax.

Then, utilize belief propagation optimized algorithm that energy function is minimized, distribute the final parallax obtaining left and right image according to the parallax in minimization of energy function.

1) information updating (see Fig. 2), its expression formula is:

$m_{p\→q}^t\left(d_q\right)=\underset{d_p\∈\mathrm{\Ω}}{min}\[V(d_p,d_q)+D_p\left(d_p\right)+\underset{r\∈N\left(p\right)\backslash q}{\Σ}{m}_{r\→p}^{t-1}\left(d_p\right)\];$

In formula, Ω represents disparity search scope, when being the t time iteration, the message that picture point p transmits to picture point q, N (p) q represent that picture point p removes three neighborhood nodes outside picture point q, comprise the x in Fig. 2 ₁, x ₂, x ₃;

Equation right-hand member Section 2 is:

$\underset{r\∈N\left(p\right)\backslash q}{\Σ}{m}_{r\→p}^{t-1}\left(d_p\right)=m_{x_1\→p}^{t-1}+m_{x_2\→p}^{t-1}+m_{x_3\→p}^{t-1};$

Wherein, when representing the t-1 time iteration, the message that r point transmits to picture point q;

Message alternating renewal process is as follows:

A, all picture points are divided into two groups: being added by each point ranks number, if odd number, is then A group; Otherwise be B group.Whole pixel set G equal the union of A and B.

Table 1 picture point divides

B, only depending on the value of picture point in its first order neighbors owing to calculating the message of picture point, so when carrying out message iteration, when iterations t is odd number, upgrading the message of B group picture point to A group picture point, the result that the message maintenance of B group picture point is t-1 time, without the need to renewal; When t is even number, upgrade the message of A group picture point to B group picture point, the message of A group picture point remains unchanged.

${\stackrel{\‾}{m}}_{pq}^t=\left\{\begin{array}{cc}{m}_{pq}^t& \left.\begin{array}{cc}\mathrm{if}& p\∈A(P\end{array}\right)\\ m_{pq}^{t-1}& other\end{array}\right.$

C, travel through the information updating of whole pixel set G.

D, return a step, carry out iterative computation.The to the last difference of twice information updating calculating, meet the demands end.

2) confidence calculations.As shown in Fig. 3 a ~ Fig. 3 b, represent the transmission of degree of confidence and the calculating of degree of confidence respectively.

The formula calculating degree of confidence is as follows:

$b_p\left(d_p\right)=D_p\left(d_p\right)+\underset{r\∈N\left(p\right)}{\Σ}{m}_{r\→p}^T\left(d_p\right);$

Wherein, b _p(d _p) represent that the parallax value that picture point p distributes is d _ptime degree of confidence, N (p) represents the four neighborhood nodes of the picture point p shown in Fig. 3 a, comprises q, x ₁, x ₂, x ₃.Wherein:

$\underset{r\∈N\left(p\right)}{\Σ}{m}_{r\→p}^T\left(d_p\right)=m_{x_1\→p}^T+m_{x_2\→p}^T+m_{x_3\→p}^T+m_{q\→p}^T;$

3) minimizing degree of confidence makes this energy function minimize, and the expression formula minimizing degree of confidence is:

$d_p^{*}=\underset{d_p\∈\mathrm{\Ω}}{\arg \min }{b}_p\left(d_p\right).$

Finally, the coordinate figure of left and right image corresponding image points in image space coordinate system is calculated according to the disparity map of left and right image.

The difference of left and right image corresponding image points coordinate figure is contained, i.e. right picpointed coordinate-left picpointed coordinate in the disparity map of left and right image)=parallax (figure); And all calculate based on left image when calculating the parallax of left and right image, namely the picpointed coordinate in left image is known, thus the picpointed coordinate that can calculate according to the disparity map of left and right image in right image, thus obtain the coordinate figure of left and right image corresponding image points in image space coordinate system, be designated as (x respectively ₁', y ₁,-f), (x ₂', y ₂,-f).

Step 13, according to the elements of exterior orientation of described left and right image and the coordinate of left and right image corresponding image points in image space coordinate system, and calculate the three-dimensional coordinate of all ground points in earth coordinates by forward intersection, thus generate intensive three-dimensional point cloud.

In the embodiment of the present invention, after carrying out overall Stereo matching by step 12, all picture points have its match point, and the elements of exterior orientation also having obtained left and right image in a step 11 (is designated as (X respectively _s1, Y _s1, Z _s1, φ ₁, ω ₁, κ ₁), (X _s2, Y _s2, Z _s2, φ ₂, ω ₂, κ ₂)), just can obtain all topocentric three-dimensional coordinate (X by forward intersection _a, Y _a, Z _a), from collinearity equation:

$\left\{\begin{array}{c}\frac{X_A-X_{S1}}{X_1}=\frac{Y_A-Y_{S1}}{Y}=\frac{Z_A-Z_{S1}}{Z_1}=N_1\\ \frac{X_A-X_{S2}}{X_2}=\frac{Y_A-Y_{S2}}{Y_2}=\frac{Z_A-Z_{S2}}{Z_2}=N_2\end{array}\right.;$

In formula, N ₁, N ₂be respectively left and right image picture point projection coefficient; (Xs ₁, Ys ₁, Zs ₁) and (Xs ₂, Ys ₂, Zs ₂) be respectively image photographic center, left and right, (X _a, Y _a, Z _a) be the three-dimensional coordinate of ground point in earth coordinates; (X ₁, Y ₁, Z ₁) and (X ₂, Y ₂, Z ₂) be respectively the coordinate of left and right image corresponding image points in the auxiliary coordinates of image space, that is:

$\begin{array}{cc}\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]=R_1\·\left[\begin{array}{c}{x}_1\′\\ y_1\\ -f\end{array}\right];& \left[\begin{array}{c}{X}_2\\ Y_2\\ Z_2\end{array}\right]=R_2\·\left[\begin{array}{c}{x}_2\′\\ y_2\\ -f\end{array}\right]\end{array};$

Wherein, R ₁and R ₂be respectively the attitude angle (φ with left image ₁, ω ₁, κ ₁), the attitude angle (φ of right image ₂, ω ₂, κ ₂) relevant rotation matrix;

Calculate the three-dimensional coordinate (X of ground point in earth coordinates _a, Y _a, Z _a), its formula is as follows:

$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\end{array}\right]=\left[\begin{array}{c}{X}_{S1}\\ Y_{S1}\\ Z_{S1}\end{array}\right]+\left[\begin{array}{c}{N}_1{X}_1\\ N_1{Y}_1\\ N_1{Z}_1\end{array}\right]=\left[\begin{array}{c}{X}_{S2}\\ Y_{S2}\\ Z_{S2}\end{array}\right]+\left[\begin{array}{c}{N}_2{X}_2\\ N_2{Y}_2\\ N_2{Z}_2\end{array}\right];$

In formula:

$\begin{array}{cc}\left\{\begin{array}{c}{N}_1=\frac{B_X{Z}_2-B_Z{X}_2}{X_1{Z}_2-X_2{Z}_1}\\ N_2=\frac{B_X{Z}_1-B_Z{X}_1}{X_1{Z}_2-X_2{Z}_1}\end{array}\right.,& \left\{\begin{array}{c}{B}_X=X_{S2}-X_{S1}\\ B_Y=Y_{S2}-Y_{S1}\\ B_Z=Z_{S2}-Z_{S1}\end{array}\right.\end{array}$

Obtain the three-dimensional coordinate of ground point in earth coordinates by calculating above again, and then generate ground point three-dimensional coordinate point of density cloud.

The such scheme of the embodiment of the present invention is by asking for left and right image elements of exterior orientation according to photogrammetric resection, then corresponding relation is set up by the coupling by pixel, acquisition can describe left and right image disparity map, finally carry out forward intersection and recover three dimensional local information, generate intensive three-dimensional point cloud; Not only solve the matching problem of image without texture region picture point, LiDAR point cloud data void holes can also be compensated, for digital city and smart city construction lay the foundation.

Through the above description of the embodiments, those skilled in the art can be well understood to above-described embodiment can by software simulating, and the mode that also can add necessary general hardware platform by software realizes.Based on such understanding, the technical scheme of above-described embodiment can embody with the form of software product, it (can be CD-ROM that this software product can be stored in a non-volatile memory medium, USB flash disk, portable hard drive etc.) in, comprise some instructions and perform method described in each embodiment of the present invention in order to make a computer equipment (can be personal computer, server, or the network equipment etc.).

The above; be only the present invention's preferably embodiment, but protection scope of the present invention is not limited thereto, is anyly familiar with those skilled in the art in the technical scope that the present invention discloses; the change that can expect easily or replacement, all should be encompassed within protection scope of the present invention.Therefore, protection scope of the present invention should be as the criterion with the protection domain of claims.

Claims (5)

1. a photogrammetric point of density cloud generation method, is characterized in that, comprising:

Photogrammetric resect is utilized to ask for the elements of exterior orientation of left and right image;

Corresponding relation between the picture point of left and right is set up by carrying out global registration to left and right image, all picture points are made to have its match point, thus the disparity map obtained for describing left and right image, and calculate the coordinate of left and right image corresponding image points in image space coordinate system based on the disparity map of left and right image;

According to elements of exterior orientation and the coordinate of left and right image corresponding image points in image space coordinate system of described left and right image, and calculate the three-dimensional coordinate of all ground points in earth coordinates by forward intersection, thus generate intensive three-dimensional point cloud.

2. method according to claim 1, is characterized in that, the described elements of exterior orientation utilizing photogrammetric resect to ask for left and right image comprises:

Based on collinearity equation, set up error equation:

$\left.\begin{array}{c}{V}_x=a_{11}{\mathrm{\ΔX}}_S+a_{12}{\mathrm{\ΔY}}_S+a_{13}{\mathrm{\ΔZ}}_S+a_{14}\mathrm{\Δ}\mathrm{\φ}+a_{15}\mathrm{\Δ}\mathrm{\ω}+a_{16}\mathrm{\Δ}\mathrm{\κ}-l_x\\ V_y=a_{12}{\mathrm{\ΔX}}_S+a_{22}{\mathrm{\ΔY}}_S+a_{23}{\mathrm{\ΔZ}}_S+a_{24}\mathrm{\Δ}\mathrm{\φ}+a_{25}\mathrm{\Δ}\mathrm{\ω}+a_{26}\mathrm{\Δ}\mathrm{\κ}-l_y\end{array}\right\};$

Wherein, Δ X _s, Δ Y _s, Δ Z _s, Δ φ, Δ ω, Δ κ are unknown number to be asked, and it represents image elements of exterior orientation (X _s, Y _s, Z _s, φ, ω, κ) correction, a ₁₁, a ₁₂..., a ₁₆, a ₂₁, a ₂₂..., a ₂₆for error equation coefficient, l _x, l _yfor constant term, V _x, V _yfor error equation correction;

The matrix form of this error equation is:

V＝AX-l；

Wherein, V=[V _xv _y] ^t, l=[l _xl _x] ^t, $A=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& a_{13}& a_{14}& a_{15}& a_{16}\\ a_{21}& a_{22}& a_{23}& a_{24}& a_{25}& a_{26}\end{array}\right],$ X＝[ΔX _SΔY _SΔZ _SΔφ Δω Δκ] ^T；

The normal equation of above-mentioned matrix is:

A ^TPAX＝A ^TPL；

Wherein, P is the weight matrix of observed reading, then the expression formula of unknown number to be asked is:

X＝(A ^TA) ^-1A ^TL；

Again iterative computation is carried out to the expression formula of unknown number above-mentioned to be asked, if the last difference calculating unknown number to be asked for twice is less than predetermined value, then calculate stopping.

3. method according to claim 1, it is characterized in that, describedly set up corresponding relation between the picture point of left and right by carrying out global registration to left and right image, all picture points are made to have its match point, thus obtain disparity map for describing left and right image, and calculate the left and right coordinate figure of image corresponding image points in image space coordinate system based on the disparity map of left and right image and comprise:

First, set up the energy function of left and right image parallax, it is expressed as:

$E\left(d\right)=\underset{(p,q)\∈N}{\Σ}V(d_p,d_q)+\underset{p\∈P}{\Σ}{D}_p\left(d_p\right);$

Wherein, d represents that the parallax of entire image distributes; N represents four neighborhood point sets of all picture points in image; d _prepresent the parallax value that pixel p distributes; P represents the set of picture point in image;

Non-set of metadata of similar data item D _p(d _p): represent that picture point p parallax is d _ptime non-similarity tolerance, calculated by following formula: D _p(d _p)=λ min (| I _l(x, y)-I _r((x+d _p), y) |, T), in formula, λ is regulatory factor, and T represents truncation error, I _l(x, y), I _r(x+dx, y+dy) is respectively image intensity corresponding to left and right image picture point;

Level and smooth item V (d _p, d _q): represent that two adjacent image points p and q distribute parallax d _pand d _qtime the discontinuous punishment amount of parallax, determined by following formula:

In above formula, a ₁, a ₂with a ₃for preset value, meet a ₃> a ₂> a ₁> 0; I _pwith I _qrepresent the intensity of picture point p and q respectively;

Secondly, utilize belief propagation optimization to calculate and this energy function is minimized, distribute the final parallax obtaining left and right image according to the parallax in minimized energy function;

Finally, the coordinate figure of left and right image corresponding image points in image space coordinate system is calculated according to the disparity map of left and right image.

4. method according to claim 3, is characterized in that, describedly utilizes belief propagation algorithm that this energy function is minimized to comprise:

Information updating, its expression formula is:

$m_{p\→q}^t\left(d_q\right)=\underset{d_p\∈\mathrm{\Ω}}{\min }\[V(d_p,d_q)+D_p\left(d_p\right)+\underset{r\∈N\left(p\right)\backslash q}{\Σ}{m}_{r\→p}^{t-1}\left(d_p\right)\];$

In formula, Ω represents disparity search scope, when being the t time iteration, the message that picture point p transmits to picture point q, N (p) q represent that picture point p removes three neighborhood nodes outside picture point q, comprise x ₁, x ₂, x ₃; when representing the t-1 time iteration, the message that r point transmits to picture point q;

$\underset{r\∈N\left(p\right)\backslash q}{\Σ}{m}_{r\→p}^{t-1}\left(d_p\right)=m_{x_1\→p}^{t-1}+m_{x_2\→p}^{t-1}+m_{x_3\→p}^{t-1};$

Confidence calculations, its expression formula is:

$b_p\left(d_p\right)=D_p\left(d_p\right)+\underset{r\∈N\left(p\right)}{\Σ}{m}_{r\→p}^T\left(d_p\right);$

B _p(d _p) represent that the parallax value that picture point p distributes is d _ptime degree of confidence, N (p) represents the four neighborhood nodes of picture point p, comprises q, x ₁, x ₂, x ₃; Wherein:

$\underset{r\∈N\left(p\right)}{\Σ}{m}_{r\→p}^T\left(d_p\right)=m_{x_1\→p}^T+m_{x_2\→p}^T+m_{x_3\→p}^T+m_{q\→p}^T;$

Minimizing degree of confidence makes this energy function minimize, and the expression formula minimizing degree of confidence is:

$d_p^{*}=\underset{d_p\∈\mathrm{\Ω}}{\arg min}{b}_p\left(d_p\right).$

5. the method according to any one of claim 1-4, it is characterized in that, the described elements of exterior orientation according to described left and right image and the coordinate of left and right image corresponding image points in image space coordinate system, and calculate the three-dimensional coordinate of all ground points in earth coordinates by forward intersection, thus generate intensive three-dimensional point cloud and comprise:

From collinearity equation:

$\left\{\begin{array}{c}\frac{X_A-X_{S1}}{X_1}=\frac{Y_A-Y_{S1}}{Y_1}=\frac{Z_A-Z_{S1}}{Z_1}=N_1\\ \frac{X_A-X_{S2}}{X_2}=\frac{Y_A-Y_{S2}}{Y_2}=\frac{Z_A-Z_{S2}}{Z_2}=N_2\end{array}\right.;$

In formula, N ₁, N ₂be respectively left and right image picture point projection coefficient; (X _a, Y _a, Z _a) be the three-dimensional coordinate of ground point in earth coordinates; The elements of exterior orientation of left and right image is designated as (X respectively _s1, Y _s1, Z _s1, φ ₁, ω ₁, κ ₁), (X _s2, Y _s2, Z _s2, φ ₂, ω ₂, κ ₂) _,(Xs ₁, Ys ₁, Zs ₁) and (Xs ₂, Ys ₂, Zs ₂) be respectively image photographic center, left and right, (φ ₁, ω ₁, κ ₁) and (φ ₂, ω ₂, κ ₂) be respectively the attitude angle of left and right image; (X ₁, Y ₁, Z ₁) and (X ₂, Y ₂, Z ₂) be respectively the coordinate of left and right image corresponding image points in the auxiliary coordinates of image space, be expressed as:

$\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]=R_1\·\left[\begin{array}{c}{x_1}^{\′}\\ y_1\\ -f\end{array}\right];$ $\left[\begin{array}{c}{X}_2\\ Y_2\\ Z_2\end{array}\right]=R_2\·\left[\begin{array}{c}{x_2}^{\′}\\ y_2\\ -f\end{array}\right];$

Wherein, the coordinate of left and right image corresponding image points in image space coordinate system, is designated as (x respectively ₁', y ₁,-f), (x ₂', y ₂,-f); R ₁and R ₂be respectively the attitude angle (φ with left image ₁, ω ₁, κ ₁), the attitude angle (φ of right image ₂, ω ₂, κ ₂) relevant rotation matrix;

Calculate the three-dimensional coordinate (X of ground point in earth coordinates _a, Y _a, Z _a), its formula is as follows:

$\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\end{array}\right]=\left[\begin{array}{c}{X}_{S1}\\ Y_{S1}\\ Z_{S1}\end{array}\right]+\left[\begin{array}{c}{N}_1{X}_1\\ N_1{Y}_1\\ N_1{Z}_1\end{array}\right]=\left[\begin{array}{c}{X}_{S2}\\ Y_{S2}\\ Z_{S2}\end{array}\right]+\left[\begin{array}{c}{N}_2{X}_2\\ N_2{Y}_2\\ N_2{Z}_2\end{array}\right];$

In formula:

$\left\{\begin{array}{c}{N}_1=\frac{B_X{Z}_2-B_Z{X}_2}{X_1{Z}_2-X_2{Z}_1}\\ N_2=\frac{B_X{Z}_1-B_Z{X}_1}{X_1{Z}_2-X_2{Z}_1}\end{array}\right.,\left\{\begin{array}{c}{B}_X=X_{S2}-X_{S1}\\ B_Y=Y_{S2}-Y_{S1}\\ B_Z=Z_{S2}-Z_{S1}\end{array}\right.$

The three-dimensional coordinate of all ground points in earth coordinates can be obtained by calculating above, thus generate intensive three-dimensional point cloud according to its three-dimensional coordinate.