Background technology
Rebuild the three-dimensional model of object in the real world, by more and more widely be applied in aspects such as virtual emulation, computer animation, reverse engineering, computer-aided design (CAD) and digital museum.Along with the continuous development of 3-D scanning equipment, the cloud data that collects based on scanning device, i.e. depth image is rebuild the model of real-world object, becomes a kind of popular three-dimensional modeling method gradually.The geometric configuration of real-world object is often complicated, and the visual angle of spatial digitizer is limited, therefore, in order to obtain the whole geological information of three-dimensional object surface, need be in a plurality of different viewpoints scannings, the depth image that will collect each time is stitched together again, reverts to a complete cloud data, and this process is exactly the registration of depth image.Because the depth image that obtains in the scanning of different points of view place is in respectively in the different coordinate systems, so key of registration two amplitude deepness images, find out the kinematic matrix between them exactly, under the effect of this kinematic matrix, all summits in second amplitude deepness image are transformed in the coordinate system of first width of cloth image.
In the research of deepness image registration method, Besl, P.J., McKay, N.D.:A method forregistration of 3-d shapes.IEEE Trans.Pattern Anal.Mach.Intell.14 (2) (1992) 239-256 and Chen, Y., Medioni, people such as G.:Object modelling by registration of multiplerange images.Image Vision Comput.10 (3) (1992) 145-155 have at first proposed ICP (the Iterative Closest Point) algorithm of registration two amplitude deepness images.The ICP algorithm is by an iterative process, constantly reduce two points converge close between predefined distance function, up to reaching certain threshold values, thereby calculate the kinematic matrix between two amplitude deepness images.The ICP algorithm is a kind of not based on the alternative manner of feature, and it needs initial motion estimated value preferably.At present, it mainly is used to finish the meticulous registration between two amplitude deepness images.
In recent years, how to realize the registration of the depth image of robotization, become the new focus of registration area research.Autoregistration for two amplitude deepness images, under the situation that does not have initial position estimation and other information (as the position of each scanner or the angle of rotation etc.), generally need to finish automatically earlier the rough registration of two amplitude deepness images, promptly estimate the initial value that moves between them; And then use ICP or other algorithms to be optimized.In the research in this respect, Chen, C.S., Hung, Y.P., Cheng, J.B.:A fastautomatic method for registration of partially-overlapping range images.In:ICCV ' 98:Proceedings of the Sixth Inter-national Conference on ComputerVision, Washington, DC, USA, IEEE Computer Society (1998) 242 utilizes the space length of point-to-point transmission as constraint condition, find out characteristic of correspondence point in two width of cloth images, calculate the approximate value of moving between this two width of cloth image then.Huber, D.F.:Automatic three-dimensional modelingfrom reality.PhD thesis, Carnegie Mellon University (2002) Chair-MartialHebert. adopts the feature on image rotating (Spin Image) descriptive model surface, then in conjunction with statistical method, rough registration two amplitude deepness images.Gelfand, N., Mitra, N.J., Guibas, L.J., Pottmann, the feature that H.:Robust global registration.In:Symposium on Geometry Processing. (2005) 197-206 then proposes to use the method for integration to come the descriptive model surface with the eigenwert that it is calculated that out, is sought out the characteristic of correspondence point again.Sun, Y., Paik, J., Koschan, A., Page, D., Abidi, M.:Point fingerprint:a new 3-d object representation scheme.Systems, Manand Cybernetics, Part B, IEEE Transactions on 33 (4) (2003) 712-717 have adopted a kind of self-defining fixed point fingerprint to seek corresponding point between two amplitude deepness images.
In above-mentioned existing each method, the method for Chen is not considered the feature of model surface, and only relies on geometrical constraint to do global search, and registration efficient is not high.Although the image rotating that Huber proposes is the feature on descriptive model surface preferably, it belongs to the feature of higher-dimension, and calculated amount is big and expend a large amount of storage spaces.The feature descriptor that people such as Gelfand and Sun proposes more complicated still on calculating, and when doing Feature Points Matching, generally rely on the information of feature self.But in actual applications, often there is similar area in model surface, only relies on local feature to be difficult to accurately distinguish these zones.Also have, existing method often adopts the method for uniform sampling or stochastic sampling when determining feature point set, and the former tends to introduce the number of characteristics point, thereby further increases the calculated amount of eigenwert; And the latter can cause the instability of the unique point quality of choosing, thereby influences final registration results.
Summary of the invention
Technology of the present invention is dealt with problems: overcome the deficiencies in the prior art, a kind of automatic deepness image registration method is provided, this method has improved the speed and the accuracy of registration, and calculated amount is little simple with calculating simultaneously.
Technical solution of the present invention: automatic deepness image registration method, its characteristics are that step is as follows:
(1) depth image is carried out trigonometric ratio and handle, convert cloud data to triangle grid data;
(2), find out the frontier point in the depth image according to described triangle grid data;
(3) eigenwert of non-border vertices in the compute depth image;
(4) according to the eigenwert of non-border vertices, extract the feature point set that has the protuberate feature in the depth image;
(5) for two amplitude deepness images subject to registration, seek the matching relationship between their unique point, i.e. corresponding point are found out the summit of at least three pairs of correspondences;
(6) calculate the estimated value of moving between this two amplitude deepness image according to described corresponding point;
(7) adopt result in the improved ICP algorithm optimization step (6), finish the accurate registration of two amplitude deepness images.
In the described step (1), the method that cloud data is converted to triangle grid data is: every amplitude deepness image is projected on the two dimensional surface, use the Delaunay triangulation of two dimension then, annexation between acquisition is had a few, again this annexation is shone upon back three dimensions, just finished the trigonometric ratio of cloud data is handled.
In the described step (2), according to described triangle grid data, the method of finding out the frontier point in the depth image is: under the situation of known triangle grid data, by traveling through whole triangle grid data, find out all isolated limits, be that those belong to a leg-of-mutton limit, the set on these isolated limits is exactly the border of grid data, and two summits that belong to these isolated limits are frontier point.
In the described step (4), the method of extracting unique point in any depth image is: at first adopt Octree segmentation depth image, then in each leaf node of Octree, the summit of selected characteristic value maximum is as the unique point of this node, the set of the unique point of all leaf nodes constitutes the feature point set of this amplitude deepness image.
Principle of the present invention: core of the present invention is choosing of unique point and mates.In fact, as long as between two amplitude deepness images subject to registration, find three pairs of characteristic of correspondence points, just can utilize their spatial coordinates calculation to go out the estimated value of the kinematic matrix between this two width of cloth image.Step of the present invention (3), (4) and (5) are exactly computation of characteristic values, selected characteristic point and the process of doing Feature Points Matching, after calculating the estimated value of kinematic matrix, just can further optimize registration results with improved ICP algorithm in the step (7).In addition,,, can effectively avoid their adverse effects like this, improve registration accuracy registration results so step (1), (2) have been rejected the frontier point in the depth image because accurately the eigenwert at computation bound point place is very difficult.
The present invention's beneficial effect compared with prior art is:
(1) the present invention rejects frontier point earlier before registration, has avoided the inaccurate eigenwert in frontier point place to the negative effect in the Feature Points Matching process, thereby has improved the efficient and the accuracy of rough registration; Simultaneously, reject speed and the accuracy that frontier point has also improved accurate registration.
(2) the present invention adopt sample point principal curvatures as feature descriptor, compare with the method for Huber, Gelfand and Sun, this feature descriptor is simple, intuitive more, is easy to calculate and relatively, makes that calculated amount is little and it is simple to calculate.
(3) feature point selection method based on Octree of the present invention on the basis that has guaranteed the selected characteristic point mass, has effectively reduced the quantity of unique point, makes that calculating is simpler.
(4) the present invention is in the matching process of unique point, combines contrast characteristic's value and two kinds of methods of how much uniformity tests of the overall situation, and the matching result of effectively having rejected the mistake that is caused by the similar area unique point has improved the accuracy of registration more.
(5) improved classical ICP algorithm, in iterative process, adjusted apart from threshold values dynamically, thereby accelerated the speed of convergence of whole algorithm, improved registration speed according to error matrix.
Embodiment
(1) depth image is carried out trigonometric ratio and handle, convert cloud data to grid data.
Depth image is a cloud data, has only comprised the coordinate information on summit, and does not have the annexation between them, i.e. topology information.Operation such as searching frontier point in the later step and calculating summit eigenwert for convenience needs carry out the processing of gridding to them earlier.
The method of gridding is: to every amplitude deepness image, establish P={P
1..., P
nBe wherein all summits.For Pi=(x
i, y
i, z
i), ignore P
iThe Z coordinate, be about to P
iProject to the P on the two dimensional surface
i'=(x
i, y
i), establish P '={ P
1' ..., P
n' the new point set that obtains for projection, use Delaunay triangulation (the B.Delauna y of two dimension then, Sur la sphere vide, Izvestia Akademii Nauk SSSR, OtdelenieMatematicheskikh i Estestvennykh Nauk, 7:793-800,1934), obtain the annexation between the summit among the P ', if P
i' and P
j' between have annexation, P then
iAnd P
jBetween also should have annexation, so the annexation on the two dimensional surface is mapped directly in the three dimensions, thereby has finished the trigonometric ratio of cloud data.
(2), find out the frontier point in the depth image according to triangle gridding.
Because near the geological information the depth image frontier point is imperfect, so be difficult to calculate accurately their eigenwert.In order to guarantee the accuracy of registration results, need before operations such as the eigenwert of calculating the summit and corresponding vertex coupling, find out and reject these frontier points.
Under the situation of known triangle grid data, by traveling through whole triangle grid data, find out all isolated limits, promptly those belong to a leg-of-mutton limit, the set on these isolated limits is exactly the border of grid data, and two summits that belong to these isolated limits are frontier point.
Concrete step is: at first define a null set E, then to each the tri patch T in the grid data
i(i=1 ..., n; N is the tri patch sum), to T
iEach bar limit e
j(j=1,2,3) judge whether comprise e among the E
j, if comprised e among the E
j, show e
jBelonged to another tri patch, so e
jNot isolated limit, in E, delete e
j, if do not comprise e among the E
j, with e
jAdd E.After so handling all tri patchs, remaining limit is the isolated limit of grid data among the E, and the summit that belongs to these isolated limits is frontier point.
(3) eigenwert of non-border vertices in the compute depth image.
The present invention adopts the derivative characteristic of body surface as feature descriptor.Curvature is the main description means of body surface derivative characteristic, wherein comprises principal curvatures, mean curvature, Gaussian curvature etc. again, owing to containing other curvature in the principal curvatures, the principal curvatures at place, employing summit is as this summit eigenwert.
If p is a point on the curved surface S, consider all curve C of ordering by p
i, can calculate the curvature k of every curve at p point place
iMaking in these curvature maximum is k
1, that minimum is k
2, k then
1, k
2The principal curvatures of having formed the p point place of curved surface S.Curvature in the classical differential geometry theory is defined on the regular surface, and for grid data, the present invention adopts surface fitting method to calculate the curvature at place, summit.At first represent curved surface, with least square method curved surface is carried out match again with polynomial approximation.For any summit p of surface mesh, k closest approach establishing p is { p
1, p
2..., p
k, treat polynomial fitting be f (x, y)=z, with the coordinate substitution f of k point (x y)=z, obtains system of linear equations a: Ax=b, and wherein A is a k * n matrix, n be f (x, y) in the number of coefficient, b be the vector of this k the z coordinate formation of putting.The least square solution of equation be f (x, y)=fitting coefficient of z.
General, use quadric surface f (x, y)=a
0+ a
1X+a
2Y+a
3Xy+a
4x
2+ a
5y
2As treating fitting surface, when calculating the principal curvatures at any summit p point place, at first quadric surface is expressed as parametric form r (u, v)=(u, v, a
0+ a
1U+a
2V+a
3Uv+a
4u
2+ a
5v
2), then calculate the first kind and second class basic parameter: the E=r
uR
u, F=r
uR
v, G=r
vR
v, L=r
UuN, M=r
UvN, N=r
VvN calculates the Weingarten matrix then
Calculate its eigenwert at last and be principal curvatures (κ to be calculated
1, κ
2).
(4), choose feature point set with protuberate feature according to the eigenwert on summit.
Use Octree to segment each amplitude deepness image, and in each leaf node of Octree, the summit of selected characteristic value maximum is as the unique point of this node, and the set of the unique point of all Octree leaf nodes promptly constitutes the feature point set of this amplitude deepness image.The number of the degree of depth of Octree decision unique point and whole registration process consuming time, promptly the degree of depth is big more, and the unique point of extraction is many more, and registration results is accurate more, but the registration time is long more.Usually, the degree of depth of Octree gets 3 or 4, can obtain registration results preferably in the rational time.Fig. 1 and Fig. 2 are respectively that to adopt the degree of depth be 3 and 4 the Octree segmentation depth image and the result of extract minutiae.
(5) for each unique point of first width of cloth image, the point that the search principal curvatures is the most close in all summits of second width of cloth image is as candidate's corresponding point, and overall geometrical constraint conditions such as employing Euclidean distance and angle, rejects wrong corresponding relation; Find out at least three pairs of characteristic of correspondence points at last.
Concrete steps are: at first, a point the principal curvatures on each summit in second amplitude deepness image being regarded as two-dimensional space, make up a kd-tree, then to all unique points in first width of cloth image, in kd-tree the search second width of cloth image in the immediate unique point of their eigenwerts, composition characteristic point is right, establishes (p
i, q
i) i=1 ..., n} is that n between this two amplitude deepness image is to the feature corresponding point; Then this n is carried out following two how much uniformity tests to the feature corresponding point:
|||p
i-p
j||-||q
i-q
j|||<ε
1
|p
i·p
j-q
i·q
j|<ε
2
The point that satisfies condition is to being character pair point, ε wherein
1And ε
2Be preassigned threshold values.
(6) calculate the estimated value of moving between this two amplitude deepness image according to these corresponding point.
The present invention adopts Horn, method among B.:Closed-form solution of absolute orientation usingunit quaternions.Journal of the Optical Society of America 4 (1987) 629-642 calculates kinematic matrix between them according at least three pairs of corresponding point between two amplitude deepness images.When exist not unique corresponding point to the time, need calculate the kinematic matrix under every kind of corresponding relation, and calculate under this kinematic matrix conversion the distance between two amplitude deepness images, find out distance kinematic matrix hour then, and with it as the result of registration roughly.
(7) use result in a kind of improved ICP algorithm optimization step (6), finish the accurate registration of two amplitude deepness images.This method is at Besl, P.J., McKay, done following improvement on the ICP classic algorithm among N.D.:A method for registration of3-d shapes.IEEE Trans.Pattern Anal.Mach.Intell.14 (2) (1992) 239-256: before doing iteration, extract the frontier point of depth image earlier, in the iterative process of back, no longer consider these points, thereby improved the accuracy of registration results; In iterative process each time, adjust apart from threshold values dynamically according to error matrix, thereby accelerated convergence of algorithm speed.
Concrete steps are as follows:
It is unit matrix that the first step, initialization, the present invention are established the transformation matrix T that accurate registration asks, and initial distance threshold ε is specified by the user.
Second step, choose reference point: in first width of cloth image, evenly choose non-frontier point as the reference point, because comprise the intensive sampling point in the depth image of scanning, in order to improve counting yield, only select some of them point to participate in accurate registration, according to the experimental result of this paper, about 10% conducts that can choose whole points are with reference to point.The 3rd step, setpoint distance threshold value: the error behind the calculating last iteration
If new distance threshold ε=3e.
The 4th step, coupling: for each reference point is mated corresponding point in second width of cloth image.Because error reduces after each iteration, when next iteration, has reduced the hunting zone like this, improves search speed.Using kd-tree closest approach in the detection range threshold range in second width of cloth image is frontier point as corresponding point if do not have point or closest approach in the threshold value, and then this reference point is a Null Spot, otherwise is the corresponding point of effective reference point.
The 5th step, the minimum error function.Calculate rigid transformation matrix T to first width of cloth image
k, make it to first width of cloth image transformation after error
Minimum.Wherein translation transformation is point set p
iAnd q
iThe translation transformation that center of gravity overlaps behind the translation transformation, uses hypercomplex number to calculate rotational transform.If p
Ixp
IyRepresent p respectively
iX, y coordinate, q
Ixq
IyRepresent q respectively
iX, y coordinate,
In order to calculate rotational transform, compute matrix
(y z), has guaranteed it is optimal transformation to the eigenvalue of maximum characteristic of correspondence vector of matrix N for w, x, and the matrix representation of establishing the represented rotational transform of q is R as the hypercomplex number q=of rotational transform
3 * 3, t
3 * 1Represent translation vector, obtain the rigid transformation matrix of current iteration by rotational transform and translation transformation
Total transformation matrix of being asked becomes T=T
kT.If error, changeed for the 2nd step greater than assign thresholds in advance, otherwise, withdraw from.
As shown in Figure 3, the registration process with the Buddha model describes.The Buddha model comprises 15 original depth image data altogether, establishes them by scanning sequency and is followed successively by R
1..., R
15, there is the coincidence zone in promptly adjacent two amplitude deepness images.If T
i(i=1 ..., 14) be R
iAnd R
I+1Between kinematic matrix, promptly pass through R
kT
K-1(k=2 ..., 14) can be with R
kTransform to R
K-1Coordinate system in.
Adopt this paper method, successively to R
iAnd R
I+1(i=1...14) do registration in twos, obtain the kinematic matrix T between them
i, pass through R then
I+1'=R
I+1T
iT
I-1... T
1, with R
I+1Be transformed into R
1Coordinate system in, R
I+1' be R
I+1At R
1The image of correspondence in the coordinate system.Merge R at last
1, R
2' ..., R
15', promptly can obtain the complete cloud data of Buddha model.
Wherein Fig. 3 a is R
1And R
2State before the registration; Fig. 3 b is R
1And R
2State behind the registration; Fig. 3 c is R
2And R
3State before the registration; Fig. 3 d is final registration results.