Path-following based Point Matching using Similarity Transformation
PPath-following based Point Matching using SimilarityTransformation
Wei LianDept. of Computer Science, Changzhi UniversityChangzhi, Shanxi, P.R. China, 046031E-mail: [email protected]
Abstract
To address the problem of 3D point matching where the poses of two point sets are unknown, weadapt a recently proposed path following based method to use similarity transformation insteadof the original affine transformation. The reduced number of transformation parameters leadsto more constrained and desirable matching results. Experimental results demonstrate betterrobustness of the proposed method over state-of-the-art methods.
Point matching is a challenging problem with applications in computer vision and patternrecognition. To solve this problem, the robust point matching (RPM) algorithm [1] usesdeterministic annealing for optimization. But it needs regularization to avoid undesirablematching results and has the tendency of aligning the mass centers of two point sets. Toaddress this issue, Lian and Zhang reduced the objective function of RPM to a concave functionof the point correspondence variable and used the branch-and-bound (BnB) algorithm foroptimization [2, 3]. These methods are more robust, but their worse case time complexityis exponential due to use of BnB. To address this issue, Lian used the path following (PF)strategy [4] to optimize the objective function of [3] by adding a convex quadratic term to theobjective function and dynamically changing the weights of the terms [5].But in the case of 3D matching, the method of [5] is experimentally shown to only performwell when the transformation is regularized, while there are problems where the poses of twopoint sets are unknown which call for matching methods where the transformations are not a r X i v : . [ c s . C V ] J a n egularized. The reason that the method of [5] performs poorly is that it uses affine transfor-mation which has large number of parameters, thus resulting in high degree of transformationfreedom and unconstrained matching results. To address this issue, we modify the method touse similarity transformation whose number of parameters is considerably smaller. Suppose we are to match two point sets X = { x i , ≤ i ≤ m } and Y = { y j , ≤ j ≤ n } in R d .For this problem, RPM uses the following mixed linear assignment − least square model:min Φ( P, s, R, t )= (cid:88) i,j p ij (cid:107) y j − sRx i − t (cid:107) − µ (cid:62) m P n = 1 (cid:62) m P (cid:101) y + s (cid:101) x (cid:62) P n − s · tr( RX (cid:62) P Y ) + 1 (cid:62) m P n (cid:107) t (cid:107) − t (cid:62) ( Y (cid:62) P (cid:62) m − sRX (cid:62) P n ) − µ (cid:62) m P n (1) s.t. P n ≤ m , (cid:62) m P ≤ n , P ≥ , (2) s ≤ s ≤ s (3)Here we use similarity transformation with R , s and t being rotation matrix, scale change andtranslation vector. The constants s ≥ s ≥ s . Thematching matrix P = { p ij } has p ij = 1 if two points i , j are matched and 0 otherwise. Thelast term in Φ is used to regularize the number of correct matches with µ being the balancingweight. (cid:107)·(cid:107) is the l norm of a vector and tr() denotes the trace of a square matrix. 1 n representsthe n -dimensional vector of all ones. The matrices X (cid:44) (cid:104) x , . . . , x m (cid:105) (cid:62) , Y (cid:44) (cid:104) y , . . . , y n (cid:105) (cid:62) andvectors (cid:101) x (cid:44) (cid:104) (cid:107) x (cid:107) , . . . , (cid:107) x m (cid:107) (cid:105) (cid:62) , (cid:101) y (cid:44) (cid:104) (cid:107) y (cid:107) , . . . , (cid:107) y n (cid:107) (cid:105) (cid:62) .It’s easily seen that given the values of P , s and R , Φ is a convex quadratic function of t .Hence, the optimal t minimizing Φ can be obtained via ∂ Φ ∂t = 0 to be (cid:98) t = (cid:62) m P n ( Y (cid:62) P (cid:62) m − sRX (cid:62) P n ) .Substituting (cid:98) t into Φ to eliminate t yields an energy function with reduced number of variables:Φ( P, s, R ) = 1 (cid:62) m P (cid:101) y − µ (cid:62) m P n − (cid:62) m P n (cid:107) Y (cid:62) P (cid:62) m (cid:107) + s ( (cid:101) x (cid:62) P n − (cid:62) m P n (cid:107) X (cid:62) P n (cid:107) ) − s · tr( R ( X (cid:62) P Y − (cid:62) m P n X (cid:62) P n (cid:62) m P Y )) (4)2
Optimal s, R minimizing Φ( P, s, R ) Let matrix A (cid:44) X (cid:62) P Y − (cid:62) m P n X (cid:62) P n (cid:62) m P Y and let
U SV (cid:62) be the singular value decomposition of A (cid:62) , where S is a diagonal matrix andthe columns of U and V are orthogonal unity vectors. Then given s >
0, the optimal rotationmatrix R minimizing Φ in (4) is (cid:98) R = U diag( (cid:104) , . . . , , det( U V (cid:62) ) (cid:105) ) V (cid:62) [6], where diag( · ) denotesconverting a vector into a diagonal matrix and det( · ) is the determinant of a square matrix.Substituting (cid:98) R into (4) to eliminate R yields a (possibly concave) quadratic program in singlevariable s . Given the range of s as s ≤ s ≤ s , one can easily solve this quadratic program bycomparing the function values at the boundary points s , s and the extreme point. P We aim to obtain an objective function only in one variable P , which can be achieved byminimizing Φ with respect to s and R , i.e.:Φ( P ) (cid:44) min s,R Φ( P, s, R ) (5)For Φ( P ), the following results can be established: Proposition 1 Φ( P ) is concave under constraints (3) . Proof:
Based on the aforementioned derivation, we have Φ(
P, s, R ) = min t Φ( P, s, R, t ). Con-sequently, we have Φ( P ) = min s,R Φ( P, s, R ) = min s,R,t Φ( P, s, R, t )For each s , R and t , Φ( P, s, R, t ) is apparently a linear function of P . We see that Φ( P ) is thepoint-wise minimum of a family of linear functions, and thus is concave, as illustrated in Fig.1. The fact that Φ( P ) is concave makes it easier for the PF algorithm to be applied to theminimization of our objective function as it requires two terms, a concave and a convex term,to be provided. Proposition 2
There exists an integer solution for any local minima (including the globalminimum) of function Φ( P ) under constraints (2) and (3) . Φ( P, s , R , t )Φ( P, s , R , t )Φ( P, s , R , t ) Figure 1: Point-wise minimization of a family of linear functions Φ(
P, s, R, t ) (dashed straightlines) with respect to parameters s , R and t results in a concave function (solid piecewisestraight line). Proof:
The polytope formed by constraint (2) satisfies the total unimodularity property [7],which means that the coordinates of the vertices of this polytope are integer valued. Wealready proved that Φ( P ) is concave under constraints (3). It is well known that any localminima (including the global minimum) of a concave function over a polytope can be obtainedat one of its vertices. Thus, the proposition follows.This result implies that minimization of Φ( P ) by simplex-like algorithms results in integervalued solution. This is important as it avoids the need of discretizing solutions which cancause error and poor performance [8].To facilitate optimization of Φ, we needs to convert P into a vector. We define the vector-ization of a matrix as the concatenation of its rows , denoted by vec(). Let p (cid:44) vec( P ). Toget the form of Φ in terms of vector p , new denotations are needed. Letvec( X (cid:62) P Y ) (cid:44) Bp, X (cid:62) P n (cid:44) Cp, Y (cid:62) P (cid:62) m (cid:44) Dp, (cid:101) x (cid:62) P n (cid:44) a (cid:62) p, (cid:62) m P (cid:101) Y (cid:44) b (cid:62) p Based on the fact vec( M M M ) = ( M ⊗ M (cid:62) )vec( M ) for any matrices M , M and M , wehave matrices B = X (cid:62) ⊗ Y (cid:62) , C = X (cid:62) ⊗ (cid:62) n , D = 1 (cid:62) m ⊗ Y (cid:62) and vectors a = (cid:101) x ⊗ n , b = 1 m ⊗ (cid:101) y .Here ⊗ denotes the Kronecker product. With the above preparation, Φ( P ) can be written interms of vector p asΦ( p ) =( b − µ mn ) (cid:62) p − (cid:62) mn p (cid:107) Dp (cid:107) + min s,R { s ( a (cid:62) p − (cid:62) mn p (cid:107) Cp (cid:107) ) − s · tr( R [mat( Bp ) − (cid:62) mn p Cpp (cid:62) D (cid:62) ]) } (6)where mat() denotes converting a vector into a matrix, which can be seen as inverse of theoperator vec(). This is different from the conventional definition.
4o facilitate optimization of Φ, we need to get the formula of the gradient of Φ. AsΦ involves minimization operations, it’s difficult to directly derive the formula of ∂ Φ ∂p . Toaddress this issue, we appeal to the result of Danskin’s theorem [9] (page 245 therein), whichin our case states that if Φ( p, s, R ) is concave in p for each s and R (this can be provedanalogously as the proof of Proposition 1) and the feasible regions of s and R are compact,then Φ( p ) = min s,R Φ( p, s, R ) has gradient: ∂ Φ( p ) ∂p = ∂ Φ( p, (cid:98) s, (cid:98) r ) ∂p = b − µ mn − (cid:62) mn p D (cid:62) Dp + (cid:107) Dp (cid:107) (1 (cid:62) mn p ) mn + (cid:98) s ( a − (cid:62) mn p C (cid:62) Cp + 1(1 (cid:62) mn p ) (cid:107) Cp (cid:107) mn ) − (cid:98) s { B (cid:62) (cid:98) r − (cid:62) mn p [ D (cid:62) ( p (cid:62) C (cid:62) ⊗ I d ) + C (cid:62) ( I d ⊗ p (cid:62) D (cid:62) )] (cid:98) r + 1(1 (cid:62) mn p ) tr( (cid:98) RCpp (cid:62) D (cid:62) )1 mn } (7)where (cid:98) s and (cid:98) R satisfy Φ( p, (cid:98) s, (cid:98) R ) = min s,R Φ( p, s, R ) . The optimal (cid:98) s and (cid:98) R can be obtainedby the method described previously. Here the vector (cid:98) r (cid:44) vec( (cid:98) R (cid:62) ) and I d denotes the d × d identity matrix. The PF algorithm [4] is used to optimize Φ by constructing an interpolation function betweena convex function (cid:107) p (cid:107) and the concave function Φ, E λ = (1 − λ ) (cid:107) p (cid:107) + λ Φ( p )and gradually increasing λ from 0 to 1 so that E λ gradually transitions from the convex function (cid:107) p (cid:107) to the concave function Φ. With each value of λ , E λ is locally minimized. We refer thereader to [5] for detail. We compare our method with state-of-the-art methods including RPM-PF [5], RPM [1], Go-ICP [10], CPD [6] and gmmreg [11]. To ensure fairness, for RPM-PF, transformation is notregularized. We implement all the methods in MATLAB on a PC with a 3.3 GHz CPU and16 G RAM. For methods only outputting point correspondence, affine transformation is usedto warp the model point set. For our method, we set parameters s = 0 . s = 1 . . , .
5] are applied to the prototype shape when generating the scene point set.The average matching accuracies (fraction of correct matches) by different methods arepresented in Fig. 3. One can see that our method performs considerably better than othermethods. This demonstrates our method’s robustness to various types of disturbances. Exam-ples of matching results by different methods in the coexisting outlier test are shown in Fig. 4.The average running times (in second) by different methods are 3.6054 for our method, 4.1525for RPM-PF, 1.8160 for RPM, 5.7997 for Go-ICP, 0.0612 for CPD and 0.2865 for gmmreg. It’sclear our method is efficient.
We proposed a PF based point matching method in this letter by adapting the method of[5] to use the similarity transformation. Due to nonlinearity of 3D similarity transformation,this is a nontrivial extension of the method of [5]. Experimental results demonstrate betterrobustness of the proposed method over state-of-the-art methods.
References [1] Chui, H., Rangarajan, A.: ‘A new point matching algorithm for non-rigid registration’.
Computer Vision and Image Understanding , 2003, , pp. 114-141[2] Lian, W., Zhang, L.: ‘Robust point matching revisited: a concave optimization approach’. European conference on computer vision , 20126 .02 0.04 0.06
Degree of deformation A cc u r a cy Deformation test
Noise level A cc u r a cy Noise test
Outlier to data ratio A cc u r a cy Outlier test
Occlusion to data ratio A cc u r a cy Occlusion test
Outlier to data ratio A cc u r a cy Coexisting outlier test oursRPM-PFRPMGo-ICPCPDgmmreg
Figure 3: Average matching accuracies by different methods in the 5 categories of tests. Theerror bars indicate standard deviations of the methods over 100 random trials.Figure 4: Examples of matching results by (from left to right and from top to bottom) ourmethod, RPM-PF, RPM, Go-ICP, CPD and gmmreg in the coexisting outlier test.[3] Lian, W., Zhang, L.: ‘Point matching in the presence of outliers in both point sets: Aconcave optimization approach’,
IEEE Conf. Computer Vision and Pattern Recognition ,2014, pp. 352-359[4] Liu, Z. Y. and Qiao, H.: ‘Gnccp-graduated nonconvexity and concavity procedure’.
IEEE rans. Pattern Analysis and Machine Intelligence , 2014, , pp. 1258-1267[5] Lian, W.: ‘A path-following algorithm for robust point matching’. IEEE Signal ProcessingLetters , 2015, , pp. 89-93[6] Myronenko, A., Song, X.: ‘Point set registration: Coherent point drift’. IEEE Transactionson Pattern Analysis and Machine Intelligence , 2010, , pp. 2262-2275[7] Papadimitriou, C.H., Steiglitz, K.: ‘Combinatorial optimization: algorithms and complex-ity’, Dover Publications, INC. Mineola. New York , 1998[8] Jiang, H., Drew, M. S., Li, Z. N.: ‘Matching by linear programming and successive convex-ification’.
IEEE Trans. Pattern Analysis and Machine Intelligence , 2007, , pp. 959-975[9] Bertsekas, D. P.: ‘convex analysis and optimization’. Athena Scientific, Belmont, Mas-sachusetts, 2003[10] Yang, J., Li, H., Jia, Y.: ‘Go-icp: Solving 3d registration efficiently and globally opti-mally’. IEEE International Conference on Computer Vision , 2013[11] Jian, B., Vemuri, B. C.: ‘Robust point set registration using gaussian mixture models’.
IEEE Trans. Pattern Analysis and Machine Intelligence , 2011, , pp. 1633-1645[12] Zheng, Y., Doermann, D.: ‘Robust point matching for nonrigid shapes by preserving localneighborhood structures’. IEEE Trans. Pattern Analysis and Machine Intelligence , 2006,28