IROS2019国际学术会议论文集 1673

Robust Non Rigid Point Set Registration Algorithm Considering Anisotropic Uncertainties Based on Coherent Point Drift Zhe Min Jin Pan Ang Zhang Max Q H Meng Fellow IEEE Abstract Non rigid point set registration PSR is an out standing and fundamental problem in fi elds of robotics com puter vision medical image analysis and image guided surgery The aim of a non rigid registration problem is to align together two point sets that have been deformed We have derived and presented a novel registration algorithm that non rigidly registers two point sets together The assumption of isotropic localization error is shared in the previous non rigid registra tion algorithms In this paper the position localization error is generalized to the anisotropic cases which means that the error distribution is not the same in different spatial directions The motivation of considering the anisotropic characteristic is that the point localization error is actually different in three spatial directions in real applications Mathematically the diffi culty in dealing with the anisotropic error case comes from the change from a standard deviation that is a scalar to a covariance matrix The formulas for updating the parameters in both expectation and maximization steps are derived In the expectation step we compute the posterior probabilities that represent the correspondences between points in two PSs In the maximization step given the current posteriors the covariance matrix of the position localization error and the non rigid transformation are updated To facilitate the proposed algorithm the low rank approximation variation of our method is also presented We have demonstrated through experiments that the proposed algorithm outperforms the state of the art ones in terms of registration and accuracy and robustness to noise More specifi cally most of the experimental results have passed the statistical tests at the 5 signifi cance level Index Terms Anisotropic Noise Image Guided Surgery Non Rigid Point Set Registration Medical Imaging I INTRODUCTION Point set registration PSR is an important and fundamen tal problem in the fi eld of robotics medical image analysis computer vision and image guided surgery 1 2 Broadly speaking PSR can be divided into rigid PSR and non rigid PSR In terms of the rigid point set registration the parameters are a few and the problem is rather simple 3 In our recent work the normal vectors are adopted to assist the rigid registration and signifi cantly improve the registration s performance 3 4 Non rigid point set registration is still a challenging problem because the transformation between the two point sets PSs is non rigid and unknown All the authors are with the Robotics Perception and Artifi cial Intel ligence Lab in the Department of Electronic Engineering The Chinese University of Hong Kong N T Hong Kong SAR China Max Q H Meng is also with the Shenzhen Research Institute of the Chinese University of Hong Kong in Shenzhen China ThisprojectispartiallysupportedbytheHongKongRGC GRF grants 14210117 RGC NSFC RGC Joint Research Scheme N CUHK448 17 and the shenzhen Science and Technology Innovation projects JCYJ20170413161616163 awarded to Max Q H Meng Corresponding Author max meng cuhk edu hk Non rigid registration has many applications in many different research fi elds In robotics the aim of the dynamic fusion task is to reconstruct and track the dynamic non rigid scenes in real time 5 One critical step in this task is to non rigidly warp the canonical volume into the live frame In this step the canonical volume s zero level set is non rigidly registered with the live frame s depth map In the framework proposed by Newcombe et al the dense non rigid ICP algorithm is utilized to fi nish this step 5 In the image guided liver surgery for example the pre operative image space has to be non rigidly aligned with the intra operative space 6 In the neurosurgery procedures such as the tumor resection the so called brain deformation or brain shift has a signifi cant impact on the navigation accuracy To compensate for this deformation in neurosurgery Bayer et al 7 have proposed a novel method to register the vessel centerlines derived from pre operative C Arm cone beam CT CBCT images to the intra operative ones Besides non rigid registration is very essential in constructing the statistical shape models SSMs In terms of the rigid registration we have already consid ered the anisotropic characteristic of the point localization error in the generalized point set registration i e the normal vectors that can be extracted from the point set are also used in the registration 3 4 In 4 the positional localization error is assumed to be isotropic The positional localization error is generalized to be anisotropic in 3 The above two algorithms are also based on the coherent point drift CPD method 8 Under the ICP framework there are several methods that obey the assumption of anisotropic localiza tion error 9 10 In terms of the non rigid registration algorithms as far as we know the point localization error in the related literature is all assumed to be isotropic The most classic non rigid method is the coherent point drift CPD method All the points are assumed to move coherently in the CPD method There exists plenty of work to improve the CPD method in various aspects 11 In this paper on the other hand we signifi cantly improve the CPD method of the robustness to the anisotropic noise in the non rigid registration We believe that the proposed algorithm is essential since the positional noise is actually anisotropic in many biomedical applications In addition in both expectation and maximization steps the matrix forms associated with the solutions are developed The advantage of the matrix forms is that the algorithm will be rather easy and faster to be implemented in a practical sense 2019 IEEE RSJ International Conference on Intelligent Robots and Systems IROS Macau China November 4 8 2019 978 1 7281 4003 2 19 31 00 2019 IEEE7897 A Organizations of This Paper The rest of this paper is organized as follows Section II describes the relevant work Section III describes the problem formulation More specifi cally the non rigid registration is formally re defi ned in a probabilistic manner considering the anisotropic positional error Section IV fi rst presents an overview of the expectation maximization framework and then introduces the expectation E step maximiza tion covariance M Cov step and maximization weighting matrix M W step in details A low rank approximation of our method is also described Section V describes the experiments and results Section VI gives a brief discussion including the main limitations of our method and possible extensions Section VII concludes the paper II RELATEDWORK Iterative closest point ICP method is one of the well known registration algorithms 12 Binary correspondences between points in two point sets PSs are computed us ing the nearest point strategy The spatial transformation between two PSs is then computed based on the current correspondences It is also well known that ICP method has the following drawbacks 3 1 easily converges to the local optimum 2 needs a good initial transformation 3 is not robust to noise and outliers Chui and Rangarajan developed a general framework for estimating correspondences and transformation to recover a non rigid transformation 13 In this category of registration methods the PS is represented by a density function and the registration is then regarded as a density estimation problem In 13 thin plate splines TPS are utilized to model the non rigid transformation Soft assignment and deterministic annealing are used to determine the correspondences and transformation iteratively With the drawback of higher com putational cost TPS RPM is more robust than ICP when some degree of degradations exist in the PSs 13 In the Kernel Correlation method the correlation between two kernel densities representing two PSs is maximized to re cover the non rigid transformation based on the M estimator 14 In the KC method the non rigid transformation is also modeled as a thin plate spline The most related work with our method is the coherent point drift CPD method presented in 8 Two steps are involved in CPD expectation E step and maximization M step In the E step the correspondences between points in two PSs are computed In the M step the non rigid transformation is updated under the current correspondences In CPD the displacement fi eld is regularized following the motion coherence theory MCT 8 Instead of thin plate splines in CPD Gauss radial basis function GRBF is used to model the transformation With the help of the variational calculus the authors prove that the optimal displacement fi eld has an elegant kernel form in multiple dimensions In CPD one PS is considered as the model while the other is assumed to be the data The data PS is considered to be transformed samples of the model PS In the non rigid CPD framework the model PS is deformed to be aligned with the data PS We remark that the CPD method is proved to be more robust and accurate than TPS RPM 8 There exists plenty of work improving the CPD method in many aspects particularly the robustness to outliers To perform the registration they introduced a L2E registration algorithm 11 The L2E algorithm is robust to a very large proportion of outliers even up to 90 In 15 both global and local structures among PSs are exploited to improve the algorithm s robustness to deformation noise outliers rotation and occlusion More recently Ma et al developed an effi cient approach based on maintaining the local neigh borhood structures of potential true matches 16 It provides a closed form solution that can accomplish the mismatch removal from thousands of putative correspondences in only a few milliseconds We should note that all the above registration methods share one common strong assumption may not be realistic of isotropic localization error associating with the observed PS As we have emphasized above this paper focuses on signifi cantly improving the CPD method s robustness to anisotropic noise A Motivations of This Paper In this work we are particularly interested in enhanc ing the non rigid registration algorithm s robustness to the anisotropic noise The point localization error is assumed to be anisotropic which means that the error distribution is different in the two 2D or three 3D directions 17 For example in the stereo camera system that is often used to reconstruct the surface of an organ the localization error is larger in the viewing direction of the camera than the other two directions 6 Mathematically the original standard deviation is replaced with the symmetric covariance matrix when we construct the mixture models B Contributions of This Paper Our contributions in this paper are summarized as follows Firstly as far as we know the point localization error is assumed to be anisotropic in the non rigid point set registration for the fi rst time We validate the robustness of our proposed algorithm with varying percentages of noise anisotropy Extensive experimental results demonstrate that our algorithm is signifi cantly more robust to anisotropic noise than the state of the art registration methods Secondly to speed up the computation time we provide the matrix form of the expression of the updated covariance matrix and the weighting matrix in the M Step In addition the low rank approximation of the algorithm is provided We also present various ways to further improve our algorithm in various aspects III PROBLEMFORMULATION The problem of seeking the non rigid transformation ma trix between two point sets PSs is mathematically indeed an optimization problem Generally speaking every optimiza tion problem owns its explicit or implicit objective function We have to minimize the value of a certain objective function 7898 to seek a specifi c non rigid transformation matrix In this paper we are warping the model shape Y y1 yM Rd Mto the observed shape X x1 xN Rd N M and N are the number of points in Y and X respectively and d is the dimension of the PSs We consider the points in Y as the Gaussian Mixture Model GMM centroids while the points in X as observed points generated by the GMM The warped model PS is defi ned as the following T Y Y v Y 1 where v Y represents the motion model of points in Y More specifi cally we are using the v i e displacement function to represent the non rigid transformation The probability density function PDF of a observed point xn being generated by the m th model point ym is defi ned as the following p xn zn m 1 2 3 2 1 2 e 1 2 zmn T 1 zmn 2 where xn Rd 1and ym Rd 1are the column vectors representing the n th and m th points locations in X and Y respectively zmn xn ym v ym represents the distance vector between the observed point xnand the m th warped model point ym v ym Sd dis the covariance matrix associated with the point localization error and Tdenote the determinant and the transpose of a matrix respectively The GMM s PDF p xn Y is then formulated as the following p xn Y w 1 N 1 w M X m 1 1 M p xn zn m 3 where w R represents the weight of the uniform distribu tion i e p xn zn M 1 1 N that is added to account for the outliers p xn zn m is defi ned in 2 v includes the model parameters Z z1 zN are the latent variables indicating the correspondences By assuming all the observed points xnare independent with each other the likelihood L defi ned as the following is maximized to estimate L n Y i 1 p xn Y 4 To effectively deal with the point set registration with more complex transformation one regularization term that is de fi ned as the following is added p v e 2 v 5 where v v 2 H Pv 2is the norm of v in the Re producing Kernel Hibert Space RKHS is a positive real number and the operator P extracted the high frequency content part Usually instead of maximizing the likelihood the negative log likelihood is minimized to estimate Q ln p v N X n 1 lnp xn Y 2 v 2 H N X n 1 lnp xn Y 6 where p xn Y is defi ned in 3 Algorithm 1 Anisotropic Point Set Registration Based on Coherent Point Drift Anisotropic CPD 1Inputs Point Sets X Y 2Outputs v the displacement function 3Initialization Initialize 0 v0 0 4Initialization Initialize W and 5while not converged do 1 E step Compute posterior probability matrix P in 9 given current 2 M Cov step Update using 11 3 M W Step Solve the following equation for W diag P1 GW W PXT diag P1 YT 6end while 7return The warped model PS as T Y WTGT IV EXPECTATIONMAXIMIZATIONALGORITHM Several ways exist to estimate the parameters of mixture models EM algorithm gradient descent and the variational inference 18 Since we are dealing with an optimization problem with missing data i e the correspondence prob abilities we cannot solve the optimization directly The maximum likelihood ML function is minimized under the EM framework More specifi cally in the expectation step the correspondence probabilities are computed In the maxi mization step the non rigid transformation matrix is updated given the current probabilities The non rigid transformation matrix is actually computed by updating the weighting matrix W The procedures of our proposed method are summarized in Algorithm 1 By expanding 6 using 2 and 3 and omitting the terms independent of the complete data log likelihood becomes Q 1 2 N X n 1 M X m 1 pmnzT mn 1zmn 1 2NPlog 2 Pv 2 7 where pmn Pold zn m xn is the posterior probability NP PN n 1 PM m 1pmn The fi rst two terms in 7 are the negative log likelihood function represents the weight of the regularization term 1 A E Step In this step we compute the posterior probabilities pmn between points in X and Y We should note that the warped model point can be directly acquired in the M W step which updates the weighting matrix The posterior probability is computed by applying the Bayes rule P zn m Y old pmn P m p xn zn m p xn Y old 8 1Or the trade off parameter between the negative log likelihood value and the regularization term 7899 where olddenotes the model parameters in the last EM step the prior probability is defi ned as P m 1 M By sub stituting 2 and 3 into 8 we get the explicit expression of the posteriors pmn 1 w e 1 2 zmn T 1 zmn wM N 2 3 2 1 2 1 w PM m 1e 1 2 zmn T 1 zmn 9 The results of pmnare stored in a matrix P RM Nwith P m n pmn On the other hand the posterior probability of xnbeing an outlier point is computed as the following pM 1 n 1 PM m 1pmn B M Cov Step The covariance matrix associated with the positional localization error is updated by solving Q 0 Q in 7 The expression of is presented as the following PN n 1 PM m 1pmn xn ym W xn ym W T NP 10 where ym W ym v ym represents the warped model point using the current weighting matrix W RM d W w1 wM Tdenotes the kernel matrix wm 1 PN n 1pmn 1 xn ym v ym By substituting the expression of v ym in 22 into 10 and with some further manipulations the matrix form of in 10 is as the following Xdiag PT1 XT XPTYT XPTGW XPTYT XPTGW T Ydiag P1 YT Ydiag P1 GW Ydiag P1 GW TNP WTGTdiag P1 GW NP 11 where G RM Mis a kernel matrix or Gram matrix with elements gij G yi yj e 1 2 yi yj 2 is the kernel bandwidth that controls the local structure of the PSs v ym is derived in Appendix I C M W Step In this step we retain the terms that are related with v in 7 Q v in 7 then becomes 1 2 N X n 1 M X m 1 pmn zmn T 1zmn 2 tr WTGW 12 where the term 2 Pv 2 is equal to 2tr W TGW is presented in Appendix III By expanding 12 and ignoring the terms that are irrelevant with v ym Q v includes the following with six terms 1 2 N X n 1 M X m 1 pmnxT n 1v ym z term 1 pmnyT m 1v ym z term 2 pmnv ym T 1xn z term 3 pmnv ym T 1ym z term 4 pmnv ym T 1v ym z term 5 2 tr WTGW z term 6 13 By substituting the expression of v ym in 22 into 12 the matrix form of the expression in 13 is the following Q W tr WTGTPXT 1 z term 1 term 3 tr WTGTdiag P1 YT 1 z term 2 term 4 tr WTGTdiag P1 GW z term 5 2 tr WTGW z term 6 14 The d