CN106373088A - Quick mosaic method for aviation images with high tilt rate and low overlapping rate - Google Patents

Abstract

The invention proposes a quick mosaic method for aviation images with high tilt rate and low overlapping rate. By utilizing the method, the robustness of registration of the images with high tilt rate and low overlapping rate can be improved, the mosaic error can be reduced, and the synthesis efficiency of a wide image can be enhanced. The method is implemented by the following technical scheme of firstly, calculating a position distribution relationship of the images by utilizing imaging assistance data of the aviation images, establishing an image position relationship matrix, and calculating a reference image; secondly, calculating overlapped regions between adjacent images on adjacent paths through an adjacent image path between a shortest path search image and the reference image, and extracting significant feature points of the overlapped regions; thirdly, normalizing the feature points or the regions into a vertex, constructing a multi-scale image by utilizing different adjacent distance point sets, and performing matching by applying image matching to calculate a spatial transformation matrix; and finally, calculating corresponding coordinates of synthesized image pixels in a source image by using the spatial transformation matrix, and calculating the pixel values of the synthesized image by using bilinear interpolation, thereby obtaining a final mosaic image.

Description

Tilt greatly the quick joining method of low Duplication aerial image

Technical field

The present invention relates to a kind of image processing field is with regard to photograph cartography, computer vision, image procossing and computer The image split-joint method of graphics is and in particular to a kind of big inclination low Duplication aerial chart based on significance key point diagram coupling As quick joining method.

Background technology

Image mosaic (image mosaic) technology is exactly by the several images having lap it may be that different time, The image that different visual angles obtain, is combined into the technology of a large-scale seamless high-definition picture, is by one group of mutual overlapping portion The image sequence dividing carries out spatial match be aligned, forms, after resampling synthesis, the wide viewing angle that a width comprises each image sequence information The technology of scene, complete, high-resolution new images.The image mosaic research of early stage is used for cartography of taking a picture, mainly always Be to take photo by plane in a large number or satellite image integration.In order to solve the limitation of lens shooting angle, with computer technology and The exhibition of digital image processing techniques, automatic Image Stitching technology be increasingly becoming photogrammetry, computer vision, image procossing and The study hotspot of computer graphicss, and be widely used in survey of deep space, remote sensing image processing (multispectral classification, environmental monitoring, Change detection, image mosaic, weather forecast, especially big image in different resolution synthesis, Earth Information emerging system (gis)), medical image Analysis (ct and nmr Data Integration etc.), cartography map rejuvenation, computer vision target positioning, automated quality control, video pressure Reduce the staff the fields such as code, virtual reality technology, super-resolution reconstruction.The method of image mosaic is a lot, and different algorithm steps have Different, but process substantially is identical.In general, image mosaic mainly comprises the steps that

A) Image semantic classification.Including the basic operation of Digital Image Processing, such as denoising, edge extracting, histogram treatment etc., set up The matching template of image and image is carried out with certain conversion, such as Fourier transformation, wavelet transformation etc. operation.

B) image registration.It is exactly using certain matching strategy, find out the template in image to be spliced or characteristic point and joining Examine corresponding position in image, and then determine the transformation relation between two width images.

C) set up transformation model.According to the corresponding relation between template or characteristics of image, calculate in mathematical model Each parameter value, thus set up the mathematical transformation model of two width images.

D) uniform coordinate conversion.According to the mathematics transformation model set up, image to be spliced is transformed into the seat of reference picture In mark system, complete uniform coordinate conversion.

E) merge reconstruct.Overlapping region with stitching image is carried out merging the smooth and seamless panorama sketch obtaining splicing reconstruct Picture.

The registration of adjacent image and splicing be panorama picture formation technology key, about image registration techniques research so far Existing very long history, its main method has following two: the minimum method of the luminance difference based on two width images and based on spy The method levied.Joining method using feature based template matching characteristic point.It is certain that the method allows image to be spliced to have Tilt and deform, overcome the problem that when obtaining image, axle center must be consistent, allow there is certain aberration between adjacent image simultaneously. The splicing of panorama sketch mainly includes following 4 steps: image pre-splicing, that is, determine the more accurate position that two width adjacent images overlap Put, the search being characterized a little lays the foundation.The extraction of characteristic point, that is, after essentially coinciding position determination, finds spy to be matched Levy a little.Image matrix transformation and splicing, set up the transformation matrix of image according to match point and realize the splicing of image.It is finally The smoothing processing of image.

In remote sensing technology field, can be to from the same area using the image registration techniques in image mosaic technology Two width or multiple image are compared it is also possible to have distortion ground image using image mosaic technology by what remote sensing satellite photographed It is spliced into more accurately complete image, as the foundation studied further.The research of image mosaic technology is most of at this stage Be based on satellite remote sensing images, aerial mapping image, digital camera pan-shot image, its image-forming condition generally have special will Ask, the essentially positive firing angle degree imaging of such as satellite remote sensing images, imaging platform is relatively stable；Aerial mapping image is essentially near Like positive firing angle degree imaging, image-forming range is nearer, and the image Duplication that mapping requires is also higher；The image-forming range of pan-shot image Relatively near it is desirable to image Duplication also higher, above-mentioned image-forming condition figure below image distortion is less, the visual angle between image overlapping region Difference is little.But in image is scouted, due to the restriction of many shooting conditions, such as farther out, imaging angle is larger for image-forming range, imaging Platform mobility strong, image Duplication low it is impossible to obtain high-quality image.Therefore, traditional joining method solve remote, Many difficult points are yet suffered from, the concrete problem existing has during image mosaic problem under big inclination angle, low Duplication image-forming condition:

(1) the image mosaic problem that remote, big inclination angle image-forming condition brings

Contemporary aircraft due to remote execution task needs it is necessary under big inclination angle at a distance shooting image, at a distance Big inclination angle shoots and inevitably causes pattern distortion.Reconnaissance flight device needs to carry out fast reserve due to task, not only schemes As Duplication relatively low, there is also visual angle difference between the sequence image of continuous imaging, and shooting distance farther out, inclination angle larger In the case of, the pattern distortion difference that small visual angle difference causes is also bigger.Therefore, compare and just penetrating shooting, incline greatly at a distance The two width images shooting under oblique angle, the difference of the image overlapping region being caused due to visual angle difference is become much larger.Traditional images are spelled The research connecing is mostly based on satellite application, joining method be applied to just penetrate or approximate orthograph picture splicing, for inclining greatly The distortion that oblique angle causes does not consider well, it is dfficult to apply to the splicing of aerial reconnaissance image.

The key technology of splicing is image registration.The image registration techniques of early stage mainly adopt point match method, this kind of method Speed is slow, precision is low, and usually needs manually to choose initial matching point, cannot adapt to large data volume.Image registration is figure As the basis merged, and the amount of calculation of image registration algorithm is typically very big.Traditional can resisting is affine, yardstick, brightness become The method for registering changed is matching process, the characteristic point of extraction or the feature of characteristic area based on local feature region or characteristic area Vector is the statistical nature of distinguished point based surrounding neighbors pixel, coupling when by calculate characteristic vector between similarity away from From, think when similarity distance meets threshold requirement Feature Points Matching success.Therefore, its match is successful whether depends on feature The similarity of point regional area, but in the splicing of big tilted image, tilt to cause changing greatly of image local area, only rely on The similarity mode of regional area is difficult to successfully.

(2) the image mosaic problem that low Duplication imaging brings

Duplication between the image being obtained due to the needs of aircraft fast reserve, imaging is relatively low, and low Duplication means figure As overlapping region is less, the effective coverage putting forward characteristic point is little, and the characteristic point quantity of extraction is necessarily less, and traditional method characteristic point Matching rate less, lead to obtain feature point pairs few, there is the risk that cannot resolve space conversion matrices.

(3) wide cut aerial image splicing speed issue

Image snoop requests produce Image Intelligence in time, higher to the splicing rate request of image.But the picture of aerial image leads to Often larger, amount of images is more, consumes more time in feature extraction, coupling search, wide format images synthesis, and with The increasing of stitching image quantity, the time of consumption dramatically increases.Therefore, how to reduce the time of feature extraction, accelerate feature Join the speed of search, improve the efficiency of wide format images synthesis, be the difficult point solving the quick splicing of aerial image.

Content of the invention

It is an object of the invention to the image split-joint method solving at this stage does not consider remote, big inclination angle, low overlap The aerial image distortion that rate image-forming condition causes is serious, and the local feature of overlapping region differs greatly, and characteristic point to be matched is relatively Few problem, and the problem that aerial image picture is larger, the quantity more image mosaic time causing is longer, provide one kind to draw Enter to be imaged assistance data and calculate effective overlapping region, by the salient region Detection and Extraction key point of overlapping region, using pass Key point builds multiple dimensioned figure and the method based on figure coupling, improves the probability of success and the splicing of special image-forming condition hypograph splicing Speed, realizes tilting greatly the quick joining method of low Duplication aerial image.

In order to obtain above-mentioned technique effect, the technical solution used in the present invention is: a kind of low Duplication aerial chart of inclination greatly The quick joining method of picture is it is characterised in that comprise the steps:

Read the sequence aerial image of sensor front end collected by camera and the imaging assistance data of each image first, according to one-tenth As assistance data calculates the position distribution relation of each image, set up image positional relationship matrix；According to image positional relationship matrix Calculate the corresponding image in center of all image overlay areas, and the reference map as composite diagram coordinate system；Then pass through the shortest Routing algorithm searches for the adjacent map path between image and reference map, calculates the position between adjacent image on neighborhood paths and closes System, calculates the overlapping region between adjacent image according to position relationship and Duplication, extracts overlapping region in two width images respectively Significant characteristics point and region, set up significant characteristics set of regions, the characteristic point that concentrate characteristic area or region normalization For summit, set up the non-directed graph of different scale with the point set of different adjacencies；Carry out feature in two width images with figure feature again The figure coupling of point, obtains matching double points collection, one by one registering adjacent two width images, calculates all neighbor map two-by-two on neighborhood paths Space conversion matrices between picture, transmit, using neighborhood paths, the space conversion matrices calculating between image and reference map；With Space conversion matrices calculate respective coordinates in source images for the composograph pixel and calculate composograph using bilinear interpolation Pixel value, obtain final stitching image.

The present invention has the advantages that compared to prior art.

The present invention utilizes the imaging assistance data of aerial image to calculate the position distribution relation between image, decreases image The number of times of registration two-by-two；By calculating suitable image as the reference map of the frame of reference, reduce image and reference coordinate Average viewing angle difference between system, thus reducing the averaged deformation error of composite diagram, overcomes the figure that low Duplication imaging brings As Bonding Problem.

The present invention adopts Shortest Path Searching to calculate the adjacent map path between image and reference map, by extracting image Significance key point sets up the multiple dimensioned figure of key point, decreases the time of feature extraction, accelerates the speed of characteristic matching search Degree, improves the efficiency of wide format images synthesis, reduces wide cut aerial image splicing speed difficulty.

The present invention is directed to the problem that prior art exists, and the imaging assistance data introducing aerial image shortens big image mosaic Time, using the Stable distritation characteristic when multi-angle is imaged for the significance key point, reduce the spelling tilting greatly low Duplication image Connect error；Introduce the robustness that figure matching technique improves image registration.By salient region Detection and Extraction key point, solve figure The few problem of the characteristic point that causes as low Duplication；Build large-scale structure feature by using key point and based on figure coupling Method, solves the problems, such as that the local feature that remote, big inclination angle imaging the causes coupling that differs greatly is more difficult；It is imaged by introducing Assistance data calculates image overlapping region to reduce the feature extraction and matching time, solves the quick Bonding Problem of aerial image.Make Calculate the space conversion matrices between adjacent image on adjacent map path one by one with figure matching process, transmit meter using neighborhood paths Space conversion matrices between nomogram picture and reference map；Finally calculate composograph pixel in source images using space conversion matrices In respective coordinates and using bilinear interpolation calculate composograph pixel value, obtain final stitching image, solve long distance The image mosaic problem brought from, big inclination angle image-forming condition.

Brief description

In order to be more clearly understood that the present invention, now pass through embodiment of the present invention, referring concurrently to accompanying drawing, to describe this Bright, wherein:

Fig. 1 is the image mosaic flow chart of the present invention low Duplication aerial image of big inclination.

Fig. 2 is that the present invention calculates image positional relationship matrix schematic diagram.

Fig. 3 is the space conversion matrices flow chart that the present invention calculates between image and reference map.

Fig. 4 is the schematic diagram that the present invention describes graph structure feature.

Specific embodiment

Refering to Fig. 1.According to the present invention, read the sequence aerial image of sensor front end collected by camera and every width figure first The imaging assistance data of picture, calculates the position distribution relation of each image, sets up image positional relationship square according to imaging assistance data Battle array；According to the corresponding image in center of image positional relationship matrix calculus all images overlay area, and as composite diagram coordinate The reference map of system；Then the adjacent map path between image and reference map is searched for by shortest path first, calculate neighborhood paths Position relationship between upper adjacent image, calculates the overlapping region between adjacent image according to position relationship and Duplication, respectively Extract significant characteristics point and the region of overlapping region in two width images, set up significant characteristics set of regions, by characteristic area collection In characteristic point or region be normalized to summit, set up the non-directed graph of different scale with the point set of different adjacencies；Again with figure Feature carries out the figure coupling of characteristic point in two width images, obtains matching double points collection, registering adjacent two width images, calculate neighbour one by one Connect the space conversion matrices between all adjacent images two-by-two on path, transmitted using neighborhood paths calculate image and reference map it Between space conversion matrices；Calculate respective coordinates in source images for the composograph pixel with space conversion matrices and using double Linear interpolation calculates the pixel value of composograph, obtains final stitching image.

Refering to Fig. 2.With camera primary optical axis and ground intersection point p as picture centre, set up c-xyz and represent with aircraft center C is expressed as the northeast sky coordinate system of initial point, according to the spatial relationship of image i imaging, calculates figure using assistance data during imaging As the geographical coordinate at the center of i, then calculate the position relationship matrix of image according to assistance data, by each picture centre coordinate meter The position distribution relation of nomogram picture, sets up position relationship matrix, and wherein, o-xyz is earth axes, θ,Middle expression image i becomes As the angle of pitch and the azimuth of phase owner optical axis cp, when (x, y, h) represents image i imaging, aircraft center is in ground o-xyz Position in coordinate system.

Refering to Fig. 3.In calculating the space conversion matrix between each image and reference map, calculate adjacent on neighborhood paths first Position relationship between image；Then, the overlapping region between adjacent image is calculated according to position relationship and Duplication；Then divide Indescribably take significant characteristics point and the region of overlapping region in two width images, set up significant characteristics set of regions；By characteristic area The characteristic point concentrated or region are normalized to summit, set up the non-directed graph of different scale, structure using the point set of different adjacencies Build the multiple dimensioned figure of significant characteristics；Then the architectural feature of figure is described；Reuse figure feature and carry out characteristic point in two width images Figure coupling, obtain matching double points collection；Calculate the space conversion matrices between adjacent two width images followed by matching double points；? Calculate the space conversion matrices between all adjacent images two-by-two on neighborhood paths afterwards, calculated by the transmission of space reflection Space conversion matrices between image and reference map.

Refering to Fig. 4.In the architectural feature of description figure, with characteristic point p point in image and d closest adjacent feature point k₁,k₂,…k_d-1,k_dBetween line constitute the non-directed graph of diagram p, the up time between the side forming using the Neighbor Points of p and p Pin angle theta₁,θ₂,…θ_d-1,θ_dThe architectural feature of description figure.

In the quick splicing of the low Duplication aerial image of big inclination,

(1) s1: read sequence image and assistance data

Read n width sequence image { i_n, n=1,2 ..., n and image imaging assistance data Assume image center and aircraft center superposition.

(2) s2: calculate image positional relationship matrix r

S21: according to image i_nSpatial relationship during imaging, using assistance data me_nCalculate image i_nGeographical coordinate

S22: calculate the position relationship matrix of n width image, wherein, θ_n,Represent in the northeast with aircraft center c as initial point In its c-xyz coordinate system, the n-th width image is imaged the angle of pitch and the azimuth of phase owner optical axis cp, x_n,y_n,h_nIt is illustrated respectively in Aircraft center position coordinates during the n-th width image imaging in the o-xyz coordinate system of ground.

S221: initially set up the position relationship matrix r of a n × n, and be initialized as 0,

S222: the longitudinal center's distance calculating between image is d_h=d × h/2, wherein image picture elements resolution are d, and image is big Little wherein w is picture traverse for w × h, and h is picture altitude.

S223: according to the y of image_nN width image is sorted by coordinate from small to large, then according to distance interval d_hIt is divided into l Individual set, i.e. c_l={ i_n, l=1,2 ... l.

S224: the index of n width image is mapped to position relationship matrix r

It is first depending in same a line that image index is mapped to position relationship matrix by image latitude y-coordinate: to image i_n, n= 1,2 ..., n, if i_n∈c_l, then indexed the l row that n is mapped to position relationship matrix r, when multiple image belongs to same Set c_lWhen, multiple image is both mapped in l row, specific algorithm is as follows.

For the multiple image in l row, it is ranked up according to image longitude x coordinate, and is mapped to each row of l row In: corresponding image collection c in the l row to position relational matrix r_l, according to the x of image_nCoordinate is sorted from small to large, And from left to right map in the l row of matrix.

(3) s3: calculate suitable figure as the reference map i of the frame of reference_c

Position relationship r according to image_ijWith image i_nLatitude coordinates (x_n,y_n), calculate suitable reference map i_c.Consider Position relationship matrix r represents the spatial relationship of image, when the barycenter of chosen position relational matrix is as reference map, other images Average viewing angle difference and reference map between is less, and the anamorphose of splicing synthesis is less, circular

${\mathrm{base}}_i=\frac{1}{n}\σget\_row\_index(r_{ij}\&notequal;0)$

${\mathrm{base}}_j=\frac{1}{n}\σget\_row\_index(r_{ij}\&notequal;0)$

get_row_index(r_ij≠ 0) represent and meet r_ijCorresponding element line index in matrix r, get_colum_ when ≠ 0 index(r_ij≠ 0) represent and meet r_ijColumn index in matrix r for corresponding element when ≠ 0.Then reference map isCorresponding Image, that is,

(4) s4: calculate image i_nWith reference map i_cBetween adjacent map path

Calculate image i_nWith reference map i_cBetween adjacent shortest path, the image that it passes through is followed successively byWill The index of these images is as image i_nTo reference map i_cBetween transition matrix bang path l_n, that is,

l_n=(n₀,n₁,n₂,…,n_m,c),n_m∈{1,2,…,n}

Wherein m represents path l_nLength, index n₀Correspondence image i_n, index c correspondence benchmark image i_c, index n_mCorresponding figure Picture and index n_m+1Corresponding image is adjacent.

(5) s5: calculate image i_nWith reference map i_cBetween space conversion matrix h_n.

1) s51: calculate path l_nOn adjacent imageWithBetween position relationship

Represent y direction Shang You overlapping region between the corresponding image of adjacent lines in position relationship matrix r, remember the position between them Relation is right_left or left_right, represents x direction Shang You overlapping region, remember him between the corresponding image of adjacent column Between position relationship be top_bottom or bottom_top, computation rule is as follows: assumes imageWithIn r In corresponding element be respectivelyWith

If i₂=i₁+ 1, represent imageWithBetween position relationship be

If i₁=i₂+ 1, represent imageWithBetween position relationship be

If j₂=j₁+ 1, represent imageWithBetween position relationship be

If j₁=j₂+ 1, represent imageWithBetween position relationship be

2) s52: calculate l_nUpper adjacent imageWithBetween overlay chart pictureWith

FoundationWith adjacent imageBetween position relationshipWith Duplication ρ, calculate the weight in image Folded region, computational methods are as follows:

IfThen

IfThen

IfThen

IfThen

3) s53: extract the significant characteristics set of regions of overlay chart picture

Respectively to overlapping imageWithThe significant characteristics set of regions of extraction as followsWith

S531: using wave filter f to overlapping area imageCarry out convolution and extract feature f_i,

f_i=o*f, i ∈ o

S532: to each point i in the o of overlapping region, calculate the significance measure s of i⁽ⁱ⁾

$s^{\left(i\right)}=\underset{i\&notequal;j,j\&element;n\left(i\right)}{\σ}c(f,e^{(i,j)})d^{(i,j)}$

Wherein n (i) is 8 Neighbourhood set of point i, d_ijFeature difference for point i and point j

d^(i,j)=| f_i-f_j|²

Wherein e^(i,j)Represent the irrelevance of human eye retina, the distance dependent with block i and j, wherein c (f, e) are oxyopia Degree.

C (f, e)=1/t (f, e)

In formula, t (f, e) is contrast threshold

$t(f,e)=t_0\exp \left\{\mathrm{\α}f\frac{e+e_2}{e_2}\right\}$

F is the spatial frequency of image, and α is spatial frequency constant, and e is Eccentricity, e₂It is the half-resolution degree of eccentricity, t₀It is Minimum contrast threshold value, empirical valueα=0.106, e₂=2.3.

S533: calculate the significant characteristics s of each point i⁽ⁱ⁾Afterwards, do following threshold process:

Work as s⁽ⁱ⁾< threshold, then s⁽ⁱ⁾=0.

Work as s⁽ⁱ⁾> threshold, then s⁽ⁱ⁾=s⁽ⁱ⁾.

Significant characteristics point s by connection⁽ⁱ⁾Merge, unconnected significance feature is rejected, and obtains salient region collection a ={ s^(p), p=1,2 ..., p.

Respectively to overlapping regionWithThe step carrying out s441 to s443, extracts salient region collection

4) s54: in the multiple dimensioned in figure building significant characteristics,

Respectively to overlapping imageWithThe salient region collection of middle extraction Do following process, build the multiple dimensioned figure of significant characteristics

S541: take the salient region s in marked feature set of regions a of overlapping region o of image i^(p)∈ a, p=1,2 ..., p Center as figure summit c^(p).

$c^{\left(p\right)}=\left(\begin{array}{cc}\frac{1}{m_p}\underset{i\&element;s^{\left(p\right)}}{\&sigma;}{x}_i^{\left(p\right)}& \frac{1}{m_p}\underset{i\&element;s^{\left(p\right)}}{\&sigma;}{y}_i^{\left(p\right)}\end{array}\right)$

m_pRepresent characteristic area s^(p)The number at midpoint,Represent coordinate in image i for the pixel i.Obtain overlap Region o vertex set c={ c^(p), p=1,2 ..., p }.

S542: build multiple dimensioned neighbour using vertex set c and scheme g.

S5421: calculate summit c using k nearest neighbor algorithm^(p)The kd of ∈ c^thNeighbor vertices set(k-1) d^thNeighbor vertices SetThen summit c^(p)D_kNeighbour's point setFor

$c_{d_k}^{\left(p\right)}=c_k^{\left(p\right)}-c_{k-1}^{\left(p\right)}$

Neighbour's point setIn have d fixing point.

S5422: utilizeBuild summit c^(p)Neighbour's yardstick be d_kRadiation diagram

S5423: define yardstick set d={ d_k, k=1,2,3 ..., k }, wherein d_k=[(k-1) d, kd].Take different chis respectively Degree d_k∈ d, repeat step s4521 to s4522, build summit c^(p)Multiple dimensioned neighbour figure

S5424: all summit c in opposite vertexes collection c^(p)∈ c, p=1,2,3 ..., p are s4521 to s4523 and process, and build summit The neighbour figure g={ g of collection c^(p), p=1,2,3 ..., p }.

Respectively to overlapping imageWithIt is step s451 to s452 to process, obtain the multiple dimensioned figure of significant characteristics

5) s55: the architectural feature of description figure

Respectively to overlapping imageWithSignificant characteristics build figureWithDo following method description, obtain Architectural feature to figureWith

S551: utilize point c^(p)D_kNeighbour's figure of yardstickCalculate by point c^(p)With summit k_d-1The side connecting and point c^(p)With summit k_dThe angle on side connecting isWhereinSo that feature has rotational invariance, d is pressed from both sides AnglePoint c is formed after being sorted from small to large^(p)Characteristic vector.

$f_k^{\left(p\right)}=({\mathrm{\&theta;}}_{k_1}^{\left(p\right)},{\mathrm{\&theta;}}_{k_2}^{\left(p\right)},...,{\mathrm{\&theta;}}_{k_d}^{\left(p\right)})$

S552: to point c^(p)Multiple dimensioned neighbour's figureDo s461 process, obtain point c^(p)Multiple dimensioned Graph structure feature

S553: to all summit c^(p)∈ c, p=1,2 ..., the neighbour figure g of p^(p)∈ g is s461 to s462 and processes, and obtains notable Architectural feature f={ the f of the figure g of property feature^(p), p=1,2 ..., p.

Respectively to overlapping imageWithSignificant characteristics build figureWithIt is s461 to s463 to process, Obtain the graph structure feature of significant characteristicsWith

6) s56: characteristic matching

S561: utilize feature point setSet up kd- Tree search tree；

S5611: calculate d dimensional feature collectionIn there is maximum variance Intrinsic dimensionality k^*, as the dimension of left and right subtree division, WhereinRepresentKth_iThe variance of dimensional feature component.

S5612: by the kth of all of point in set^*Dimensional feature componentBe ranked up, take medianAs k^*Divide value in dimension

S5613: utilize k^*WithSet up root node feature set is divided, by kth^*Dimensional feature componentCharacteristic point as right subtree, by kth^*Dimensional feature component Characteristic point as left subtree.

S5614: respectively s4711, s4712, s4713 recursive operation is carried out to the point set in left subtree, right subtree, until left and right Tree is all regardless of subdivided.

S562: find characteristic point in kd-tree search tree using bbf methodNearest neighbor pointWith secondary Neighbor Points

S5621: for pointRoot node from kd-tree search tree starts search and compares, with node in search treeCompare, ifThen select sub- rightchild tree search to the right, Otherwise select to search for left subtree leftchild, and by not selected subtree subtree={ rightchild or then Leftchild } and its root node withBetween characteristic distanceIt is saved in minimum In Priority Queues q, queue q is according to the order arrangement from small to large apart from d.

The execution s5621 operation of s5622: recurrence, until searching for leaf node, then calculates pointPreserve with leafy node Characteristic pointThe distance between

$\begin{array}{c}dis(c^{\left(p_{n_m}\right)},c^{\left(p_{n_{m+1}}\right)})=\left|\right|f_k^{(n_m,p)}-f_k^{(n_{m+1},p)}|{|}_2\\ =\underset{i=1}{\overset{d}{\&sigma;}}{({\mathrm{\&theta;}}_{k_i}^{\left(p_{n_m}\right)}-{\mathrm{\&theta;}}_{k_i}^{\left(p_{n_m+1}\right)})}^2\end{array}$

The possible more than one of characteristic point that wherein leaf node preserves, preserves wherein two minimum characteristic points of distance As nearest neighbor point and time Neighbor Points, and each self-corresponding two distancesWith

S5623: select queue q in corresponding subtree subtree=q of minimum range (0), execution s5621 and s5622 operate into Row retrospective search, when the new nearest neighbor point calculatingDistance less than original nearest neighbor distance dis ' when, will be original Nearest neighbor point replace with time Neighbor Points and beNearest neighbor distance originally replaces with time nearest neighbor distance Dis "=dis ', new nearest neighbor point is replaced original nearest neighbor point isNew nearest neighbor distance is replaced Nearest neighbor distance originally is dis '=dis " '；When the new secondary Neighbor Points calculatingDistance than original secondary neighbour away from From dis, " hour, new secondary Neighbor Points are replaced original secondary Neighbor Points isNew secondary nearest neighbor distance is replaced Secondary nearest neighbor distance originally is dis "=ds " ".Subtree q (0) recalled in queue is deleted after the completion of lookup.

S5624: repeated execution of steps s5623, until it is sky that backtracking step is more than threshold value or queue q, finally gives a little Nearest neighbor pointWith secondary Neighbor Points

S563: judge when satisfaction

$\frac{{\mathrm{dis}}^{\&prime;}}{{\mathrm{dis}}^{\&prime;\&prime;}}>threshold$

When it is believed that pointFor pointMatch point, note matching double points be

S564: rightIn all characteristic pointsRepeat step s472 and s473, obtain To two width adjacent imagesWithThe matching double points set of image

7) s57: using matching double points setCalculate adjacent image using stochastic sampling concordance (rancac) methodWithBetween space conversion matrices

S571: from matching double points collectionIn randomly draw four points to collection Calculate four points between transformation matrix h, and calculate the interior point set of transformation matrices h.

$\left[\begin{array}{cccccccc}{x}_1& y_1& 1& 0& 0& 0& -x_1*x_1^{\&prime;}& -y_1*x_1^{\&prime;}\\ x_2& y_2& 1& 0& 0& 0& -x_2*x_2^{\&prime;}& -y_2*x_2^{\&prime;}\\ x_3& y_3& 1& 0& 0& 0& -x_3*x_3^{\&prime;}& -y_3*x_3^{\&prime;}\\ x_4& y_4& 1& 0& 0& 0& -x_4*x_4^{\&prime;}& -y_4*x_4^{\&prime;}\\ 0& 0& 0& x_1& y_1& 1& -x_1*y_1^{\&prime;}& -y_1*y_1^{\&prime;}\\ 0& 0& 0& x_2& y_2& 1& -x_2*y_2^{\&prime;}& -y_2*y_2^{\&prime;}\\ 0& 0& 0& x_3& y_3& 1& -x_3*y_3^{\&prime;}& -y_3*y_3^{\&prime;}\\ 0& 0& 0& x_4& y_4& 1& -x_4*y_4^{\&prime;}& -y_4*y_4^{\&prime;}\end{array}\right]*\left[\begin{array}{c}a\\ b\\ c\\ d\\ e\\ f\\ g\\ h\end{array}\right]=\left[\begin{array}{c}{x}_1^{\&prime;}\\ x_2^{\&prime;}\\ x_3^{\&prime;}\\ x_4^{\&prime;}\\ y_1^{\&prime;}\\ y_2^{\&prime;}\\ y_3^{\&prime;}\\ y_4^{\&prime;}\end{array}\right]$

Solve overdetermined equation using method of least square, obtain transformation matrix h

$h=\left[\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& 1\end{array}\right]$

S572: execution step s571n time, in record, the interior point set of the point maximum single sample of number is as concordance matching double points Collection

S573: using concordance matching double points collectionCalculate the space conversion matrices between matching double points

8) s58: calculate image i_nWith reference map i_cBetween transition matrix

To path l_n=(n₀,n₁,n₂,…,n_m, c) on adjacent two-by-two imageWithExecution s51 to s57 is processed, respectively Calculate the transformation matrix between all adjacent imagesThen image i_nWith reference map i_cBetween Transition matrix h_nFor

$h_n={\&pi;}_{m=0}^{m+1}{h}_{n_m}$

WhereinFor unit matrix.

To each image i_nDo the process of step s4, s5 respectively, calculate image i_nWith reference map i_cBetween space conversion matrices h_n, n= 1,2,…,n.

(6) s6: image synthesizes

S61: calculate each image i_nFour angles of mapping area in the composite image coordinate map_n={ (tlx_n,tly_n), (trx_n,try_n),(brx_n,bry_n),(blx_n,bly_n), wherein tlx_n=c, tly_n=f

${\mathrm{trx}}_n=\frac{a*w_n+c}{g*w_n+1},{\mathrm{try}}_n=\frac{d*w_n+f}{g*w_n+1}$

${\mathrm{brx}}_n=\frac{a*w_n+b*h_n+c}{g*w_n+h*h_n+1},{\mathrm{bry}}_n=\frac{d*w_n+e*h_n+f}{g*w_n+h*h_n+1}$

${\mathrm{blx}}_n=\frac{b*h_n+c}{h*h_n+1},{\mathrm{bly}}_n=\frac{e*h_n+f}{h*h_n+1}$

Wherein w_nFor the width of ground n width image, h_nHeight for the n-th figure.

S62: according to all of mapping area { map_n, calculate the size of composograph

Lx=min { tlx_n,blx_n}

Ty=min { tly_n,try_n}

Rx=max { trx_n,brx_n}

By=max { bly_n,bry_n}

Then the size of composograph is

Width=rx-lx

Height=by-ty

S63: calculate h_nInverse matrix hinv_n

S64: calculate mapping area map in composograph using bilinear interpolation_nIn pixel value

For mapping area map_nIn pixel (x, y) ∈ map_n, corresponding artwork i_nMidpoint is

$x^{\&prime;}=\frac{a*x+b*y+c}{g*x+h*y+1},$

$y^{\&prime;}=\frac{d*x+e*y+f}{g*x+h*y+1}$

Xint=floor (x '), x "=x '-xint, "=y '-yint, wherein floor () represent for yint=floor (y '), y Round downwards, then in composograph, the corresponding pixel value of pixel (x, y) is

F (x, y)=x " * y " * i (xint, yint)+(1-x ") * y " * i (xint+1, yint)+x " * (1-y ") * i (xint, yint +1)+(1-x″)*(1-y″)*i(xint+1,yint+1).

Claims (9)

1. a kind of quick joining method of the low Duplication aerial image of big inclination is it is characterised in that comprise the steps:

Read the sequence aerial image of sensor front end collected by camera and the imaging assistance data of each image first, according to one-tenth As assistance data calculates the position distribution relation of each image, set up image positional relationship matrix；According to image positional relationship matrix Calculate the corresponding image in center of all image overlay areas, and the reference map as composite diagram coordinate system；Then pass through the shortest Routing algorithm searches for the adjacent map path between image and reference map, calculates the position between adjacent image on neighborhood paths and closes System, calculates the overlapping region between adjacent image according to position relationship and Duplication, extracts overlapping region in two width images respectively Significant characteristics point and region, set up significant characteristics set of regions, the characteristic point that concentrate characteristic area or region normalization For summit, set up the non-directed graph of different scale with the point set of different adjacencies；Carry out feature in two width images with figure feature again The figure coupling of point, obtains matching double points collection, one by one registering adjacent two width images, calculates all neighbor map two-by-two on neighborhood paths Space conversion matrices between picture, transmit, using neighborhood paths, the space conversion matrices calculating between image and reference map；With Space conversion matrices calculate respective coordinates in source images for the composograph pixel and calculate composograph using bilinear interpolation Pixel value, obtain final stitching image.

2. as claimed in claim 1 big tilt low Duplication aerial image quick joining method it is characterised in that: with camera Primary optical axis is the center p of image with the focus on ground, sets up c-xyz and represents in the northeast being expressed as initial point with aircraft center c Its coordinate system, according to the spatial relationship of image i imaging, calculates the geographical seat at the center of image i using assistance data during imaging Mark, is then calculated the position distribution relation of image, sets up position relationship matrix by each picture centre coordinate, and wherein, o-xyz is ground Areal coordinate system, θ,Represent that image i is imaged the angle of pitch and the azimuth of phase owner optical axis cp, (x, y, h) represents image i imaging When position in the o-xyz coordinate system of ground for the aircraft center.

3. as claimed in claim 2 big tilt low Duplication aerial image quick joining method it is characterised in that: with image Middle characteristic point p point and d closest adjacent feature point k₁, k₂... k_d-1, k_dBetween line composition p point non-directed graph, make Angle theta clockwise between the side being formed with Neighbor Points and the p of p point₁, θ₂... θ_d-1, θ_dThe architectural feature of description figure.

4. as claimed in claim 1 big tilt low Duplication aerial image quick joining method it is characterised in that: inclining greatly Tiltedly in the quick splicing of low Duplication aerial image, read n width sequence image { i_n, n=1,2 ..., n } and the imaging of image auxiliary Help dataAccording to image i_nSpatial relationship during imaging, using assistance data me_nCalculate figure As i_nGeographical coordinate

Calculate the position relationship matrix r of n width image, wherein, θ_n、Represent in the northeast sky with aircraft center c as initial point In c-xyz coordinate system, the n-th width image is imaged the angle of pitch and the azimuth of phase owner optical axis cp, x_n, y_n, h_nIt is illustrated respectively in ground Aircraft center position coordinates during the n-th width image imaging in the o-xyz coordinate system of face, n represents sequence image quantity to be spliced.

5. as claimed in claim 4 big tilt low Duplication aerial image quick joining method it is characterised in that: set up one The position relationship matrix r of individual n × n, and it is initialized as 0, the longitudinal center's distance calculating between image is d_h=d × h/2, its In, image picture elements resolution is d, and image size is w × h, and wherein, w is image pixel width, h image pixel height.According to figure The latitude y of picture_nN width image is sorted by coordinate from small to large, then according to distance interval d_hIt is divided into l set c_l= i_n, l=1,2 ... l, to image i_n, n=1,2 ..., n, if i_n∈c_l, then indexed n and be mapped to relational matrix r's L row, when multiple image belongs to same set c_lWhen, multiple image is both mapped in l row；For many in l row Width image, it is ranked up according to image longitude x coordinate, and is mapped in each row of l row: the l row to relational matrix r In corresponding image collection c_l, according to the x of image_nCoordinate is sorted from small to large, and from left to right maps to the l of matrix In row.

6. as claimed in claim 5 big tilt low Duplication aerial image quick joining method it is characterised in that: calculate figure As i_nWith reference map i_cBetween adjacent shortest path, the image of process is followed successively byRope by these images Draw as image i_nTo reference map i_cBetween transition matrix bang path l_n, l_n=(n₀, n₁, n₂..., n_m, c), n_m∈ 1, 2 ..., n }, wherein, m represents path l_nLength, index n_oCorrespondence image i_n, index c correspondence benchmark image i_c, index n_mCorresponding Image with index n_m+1Corresponding image is adjacent.

7. as claimed in claim 1 big tilt low Duplication aerial image quick joining method it is characterised in that: extract During the significant characteristics set of regions of overlay chart picture, respectively to overlapping imageWithThe significance of extraction as follows special Levy set of regionsWithUsing wave filter f to overlapping area imageCarry out convolution and extract feature f_i, f_i=o*f, i ∈ o, then to each point i in the o of overlapping region, calculates the significance measure s of i⁽ⁱ⁾, s⁽ⁱ⁾= ∑_{I ≠ j, j ∈ n (i)}C (f, e^{(i, j)})d^{(i, j)}；Wherein, n (i) is 8 Neighbourhood set of point i, d_ijFeature difference d for point i and point j^{(i, j)} =| f_i-f_j|², e^{(i, j)}Represent the irrelevance of human eye retina, the distance dependent with block i and j, c (f, e) is visual acuity, e It is Eccentricity.

8. as claimed in claim 1 big tilt low Duplication aerial image quick joining method it is characterised in that: build The multiple dimensioned in figure of significant characteristics, respectively to overlapping imageWithThe salient region collection of middle extractionDo following process, build notable The multiple dimensioned figure of property featureTake the salient region s in marked feature set of regions a of overlapping region o of image i^(p)The center of ∈ a, p=1,2 ..., p is as the summit c of figure^(p),? To overlapping region o vertex set c={ c^(p), p=1,2 ..., p }, build multiple dimensioned neighbour using vertex set c and scheme g, m_pRepresent special Levy region s^(p)The number at midpoint,Represent coordinate in image i for the pixel i.

9. as claimed in claim 1 big tilt low Duplication aerial image quick joining method it is characterised in that: figure In the feature description joined, respectively to overlapping imageWithSignificant characteristics build figureWithDo as follows Process, obtain the architectural feature of figureWithTake point c^(p)With summit k_d-1The side connecting and point c^(p)With summit k_dConnect The angle on side isAnd after being sorted from small to large, form point c^(p)Characteristic vector? Arrive point c^(p)Multiple dimensioned graph structure featureTo all summit c^(p)∈ c, p=1, 2 ..., p neighbour figure g^(p)∈ g process, obtains the architectural feature f={ f of the figure g of significant characteristics^(p), p=1,2 ..., p.