Robust Nonrigid ICP Using Outlier-Sparsity Regularization
Robust Nonrigid ICP Using Outlier-Sparsity Regularization Hidekata Hontani
Takamiti Matsuno Nagoya Institute of Technology Aichi, JAPAN
Yoshihide Sawada
[email protected]
Abstract
is registered to target data, which are extracted from the images. The problem here is the existence of nuisances: False target points would be extracted from anatomical structures that are not targeted. In addition, some local but large deformation would occur because of lesions near the target surface, and some true target points would be missed because of, e.g., the contact with neighboring organs. In such situations, it is known that statistical shape model (SSM) is required for registering the template robustly[10]. The active shape model (ASM) is widely used for the registration[6][7]. In the ASM, the template surface is represented by a point distribution model (PDM), which represents a surface by using a set of n points on it. Let wi denote the 3D coordinates of the point Pi (i = 1, 2, · · · , n), and let w ≡ [wT1 , wT2 , · · · , wTn ]T denote the global shape of the template. Given a set of m training surfaces, {w1 , w2 , · · · , wm }, one constructs a linear shape model of w for ASM by applying PCA to the covariance matrix of the training surfaces:
We show how to incorporate a statistical shape model into the nonrigid ICP framework, and propose a robust nonrigid ICP algorithm. In the nonrigid ICP framework, a template surface is represented by a set of points, and the shape of the template is parametrized by a transformation matrix per one template point. In the proposed method, the statistics of the matrices are estimated based on a set of training surfaces, and the statistical shape model is incorporated into the nonrigid ICP framework by modifying the representation of the stiffness of the template. The statistical shape model and a noise model make it possible to discriminate outliers from inliers in given targets. Our proposed method detects the outliers, which are not represented by the models appropriately, based on their sparseness. The detected outliers are automatically excluded from the target to be registered, and the template is deformed to fit the inliers only. As the result, the accuracy of the registration is improved. The performance of the proposed method is evaluated qualitatively and quantitatively using synthetic data and clinical CT images.
¯ + P θ, w=w
(1)
¯ is the mean of the training surfaces, P is a 3n × k where w (k < 3n) matrix each of which column is the eigen vector of the covariance matrix, and θ is a k-vector of the deformation parameters. The ASM iteratively updates the estimate of θ that minimizes the sum of the squared residuals between the point Pi on the template and its corresponding point on a given target. It is known that the least square approach is optimal only when the residuals obey the Gaussian distribution. This Gaussian model, though, does not always describe the distribution of the residuals appropriately, because of the nuisances: They often become the outliers of the Gaussian model. Many methods, hence, are proposed for making the ASM robust: Outliers and inliers are automatically discriminated, and the template is deformed to fit the inlier targets in ways of consistent with the SSM[10]. Duta and Sonka[9], for example, proposed a robust ASM that detects outliers based on the variance information from the PDM. Lekadir et al.[12] proposed another robust method that de-
1. Introduction The registration of 3D surfaces is an important task in computer vision. It has broad applications, such as segmentation of anatomical structures in 3D medical images, face morphing, 3D shape retrieval and object recognition[5]. Registering a model surface (template) to a given measured surface (target) means finding a mapping between the template surface and the target surface that describes the position of semantically corresponding points[1]. Numerous algorithms for the model-to-data registration have been proposed, and it is difficult to focus on some specific drawbacks and advantages of those algorithms without limiting their applications. The main application of the proposed method is the segmentation of the anatomical structures in given 3D medical images: The surface template of each anatomical structure 1
tects outliers based on the ratio of landmark distances. The techniques of the statistical sampling or of the robust point matching are also used for making ASM more robust[13]. The iterative closest point (ICP)[3] is also the algorithm that registers the PDM to given measured targets, and many modifications to the original rigid ICP are proposed to improve the robustness against outliers. Phillips et al.[14], for example, proposed the Outlier Robust ICP, which registers the template by using a subset of the template points with the smallest residual distances. One of the key points of this method is that the fraction, f ∈ [0, 1], of the subset can be adaptively updated in the registration process. Horaud et al.[11], proposed another robust ICP, which can explicitly detect outliers by means of the EM algorithm. A latent variable is introduced for each template point in order to describe if the point corresponds to the outlier or not, so that the template is registered to the inlier target points while the outliers are excluded. These are, though, the rigid registration method and no SSM is used. The nonrigid ICP (NICP)[1] is the another modification to the original rigid ICP. The deformation of the template is represented by an affine transformation matrix per one vertex, and the template is registered by estimating the optimal matrices that minimizes a cost function. The details will be described below. One of the advantages of the NICP is that the optimal transformation matrices can be computed directly when the corresponding points on the target are fixed, because the cost function is a quadratic system of the unknown matrices. The cost function includes a regularization term that requires the template to be stiff, but no SSM is used for restricting the shape of the template. If the aim is to register the template to anatomical structures, the specification of the template[8] is not enough when its shape is restricted only by the stiffness. In this paper, we incorporate the SSM of the template surface to the NICP. The SSM does not directly represent the statistics of the coordinates of the template points but of the statistics of the transformation matrices. Given the SSM and a noise model, we propose a robust NICP algorithm, which detects outliers based on their sparseness and which registers the template only to the inliers. This framework can be found in the optical flow computation proposed by Ayvaci et al.[2]: The occlusion detection and the optical flow computation are simultaneously realized based on the sparse distribution of the occlusions, at which the optical flow constraint is not satisfied. Before describing our proposition, the outline of the NICP is explained in the next subsection.
1.1. Nonrigid ICP Algorithm The template S = (V, E) is given as a set of n vertices, V, and a set of edges, E. Let the homogeneous coordinates of the vertex i be denoted by v i = [xi , yi , zi , 1]T
(i = 1, 2, · · · , n). The deformation of the template is parametrized by one affine 3 × 4 transformation matrix, Xi , per template vertex, v i . The unknowns are organized in a 4n × 3 matrix: X ≡ [X1 X2 · · · Xn ]T .
(2)
The template is registered to the given target surface, T , by minimizing the following cost function: E(X) ≡ Ed (X) + αEs (X) + βEl (X).
(3)
Here, Ed (X) is a sum of squared distances between the target surface, T , and the template vertices, Ti v i , such that ∑ Ed (X) ≡ wi dist2 (T , Xi v i ), (4) v i ∈V where wi is a positive weight coefficient and dist(T , Xi v i ) is the distance between Xi v i and its closest point on T . Es (X) in (3) is a stiffness term, which regularizes the deformation. The weighted difference of the transformations of neighbouring vertices is penalized as follows: ∑ Es (X) ≡ k(Xi − Xj )Gk2F , (5) (i,j)∈E
where k · kF denotes the Frobenius norm and G is a weighting 4×4 matrix, which is needed for balancing the scales of the elements of Xi : G ≡ diag(1, 1, 1, s). The coefficient, α, in (3) controls the stiffness of the template S, which becomes more stiff when α is larger. El (X) in (3) is a landmark term, which is used for the initialization and the guidance of the template. Let a set of landmarks be denoted by L ≡ {(v i1 , l1 ), (v i2 , l2 ), · · · , (v il , ll )}. The distance between the landmark v i on S and the corresponding one ls on T is penalized as follows: ∑ kXi v i − li k2 . (6) El (X) ≡ (v i ,li )∈L When the closest points on T are fixed, one can directly compute the exact minimizer of the cost function E(X) in (3). Let ui denote the point on T that is closest to Xi v i . Then, Ed (X) in (3) becomes ∑ ¯d (X) ≡ E wi kXi v i − ui k2 (7) v i ∈V 2
= kW (DX − U )kF , √ √ where W ≡ diag( w1 , · · · , wn ), D is an n × 4n sparse matrix: T v1 v T2 D≡ (8) , .. . v Tn
and U ≡ [u1 u2 · · · un ]T . It should be noted that Eq. (7) is a sparse quadratic system of the unknown parameters ¯d , Es and El in (3) can be represented in X. Similar to E quadratic forms. The stiffness term, Es , becomes as follows: Es (X) = k(M ⊗ G)Xk2F , (9) where M is the node-arc matrix. When two vertices, i and j, are connected by an edge, r, then the nonzero entries of M in row r are (M )ri = −1 and (M )rj = 1. The landmark term, El , becomes El (X) = kDL X − UL k2F ,
(10)
where DL is a matrix obtained by taking the rows of D corresponding to the landmark vertices, and UL = [l1 , · · · , ll ]T . As the result, the cost function E(X) in (3) can be represented by a quadratic function of X:
2
WU
√ WD
¯ α · M ⊗ G X − 0 E(X) =
√ √
β · DL β · U L F (11) = kAX − Bk2F . ¯ One can directly minimize E(X) in (11) by computing T −1 T ˆ X = (A A) A B. The NICP algorithm can be described as follows: • Initialize X 0 . • For each αj ∈ {α1 , α2 , · · · , αK }, αj−1 > αj . ˆj − X ˆ j−1 kF < . • Until kX ˆ j−1 v i . • Find corresponding point uj−1 for each X i i j ˆ that minimizes E(X). ¯ • Compute X In the registration process, the stiffness coefficient parameter, α, is controlled: Starting with a more stiff model, the global alignment is recovered firstly, and then successively lower stiffness coefficient is used[1].
2. Method As mentioned above, a statistical model of the template shape is incorporated to the NICP method in this article. The statistical model is constructed based on a set of training surfaces. This model is useful for distinguishing true target points from false ones in given data. In addition, the model provides a criteria for detecting the outliers.
2.1. Nonrigid ICP with Statistical Shape Model Same with the NICP, the proposed method represents the template surface by using a set of n points, {Pi |i = 1, 2, · · · , n}, and its deformation is parametrized by one affine 3 × 4 transformation matrix, Xi , per template vertex. We firstly explain how to construct the probabilistic model of {Xi } from a given set of training surfaces, and show the registration algorithm. Then, we propose the robust registration method, which automatically detects outliers based on their sparsity.
2.1.1
Construction of Statistical Shape Model
Let S k (k = 1, 2, · · · , m) denote the kth training surface. For constructing the probabilistic model of {Xi }, the average shape of {S k } is firstly computed, and then each of the training surfaces is registered to the average shape by means of the NICP. Let Xik denote the resultant matrix obtained by the NICP. Then, the probabilistic model of Xi is constructed based on {Xik |k = 1, 2, · · · , m}. ¯ a set of n correFor computing the average shape, S, k sponding points, {Pi |i = 1, 2, · · · , n}, is generated on each surface, S k . Let wki denote the 3D coordinates of ¯ is the polygonal surface, Pik . The average surface, S, ¯ whose vertices, Pi , are determined by computing the mean ¯i = coordinates for ∑ each of the corresponding points: w ∑ k k ¯ k w i /m = ( k Xi /m)v i = Xi v i . Then, registering S k to S¯ by means of the NICP, we obtain a set of transformation matrices, {Xik |i = 1, 2, · · · , n}. Let the transformation matrix, Xi , be vectorized to a 12vector, xi by vertically concatenating the columns of XiT , and let x ≡ [xT1 , xT2 , · · · , xTn ]T . In the proposed method, two probability distributions, p(xi ) and p(xi −xj ), are estimated based on the transformation matrices, {Xik }. We represent these two distributions by using the 12-dimensional multivariate Gaussian distributions, N (·|µ, Σ), where µ denotes the mean and Σ denotes the covariance. Let p(xi ) ∼ N (xi |µi , Σi ) and let p(xi − xj ) = N (xi − xj |µij , Σij ). These parameters, µi , µij , Σi , and Σij are estimated based on {Xik }. Here, p(xi − xj ) is estimated for the pair of Pi and Pj , if there exists an edge between them in the polygo¯ nal surface, S.
2.1.2
Optimal Step for Nonrigid ICP with Statistical Shape Model
The probabilistic models, p(xi ) and p(xi − xj ), are incor¯ porated to the NICP by modifying the cost function E(X) in (11) as follows. ¯d (X) + αE ¯s (X) + βEl (X) + γEp (X). (12) E0 (x) = E ¯d (X) Here, the first and the third terms of the right side, E and El (X), are shown in (7) and (6), respectively. Using x ¯d (X) as instead of X, one can rewrite E ¯d (x) = kW 0 (D0 x − u)k2 , E
(13)
where W 0 is a 3n × 3n matrix and is defined as √ W0 ≡
w1 · I3 ..
.
√
, wn · I3
(14)
in which Id denotes the d × d identity matrix. D0 in (13) denotes a 3n × 12n matrix and is defined as T v 1 ⊗ I3 .. (15) D0 ≡ . . v Tn ⊗ I3 u in (13) denotes a 3n-vector: u = [uT1 , uT2 , · · · , uTn ]T . ¯d , El (X) in (12) can be written as follows: Analogous to E 0 El (x) = kDL x − uL k2 ,
(16)
0 where DL and uL are obtained by taking the elements that corresponds to the landmarks from D0 and u in (13), respectively. ¯s (X) in (12), is a modified version of The second term, E the stiffness term, which penalizes local deformations that are statistically rare, and is defined as follows: ∑ ¯s (X) ≡ E (xi − xj − µij )T Σ−1 ij (xi − xj − µij ) (i,j)∈E
=
∑
kGij x − Gij µk2
(i,j)∈E
= kGx − Gµk2 .
(17) −1/2
⊗ Mij , Here, µ ≡ [µT1 , µT2 , · · · , µTn ]T and Gij ≡ Σij where Mij is an n-dimensional row vector, in which ith element is one, jth element is minus one, and the others are zero. Concatenating all Gij vertically, we obtain a 12|E| × 12n matrix G in (17). The fourth term, Ep (X) in (12), is introduced for penalizing the distance between the ¯ i , as transformed location, Xi v i , and the mean location, w follows. ∑ Ep (X) ≡ (xi − µi )T Σ−1 i (xi − µi ) i
= kΛx − Λµk2 , where
Λ≡
(18)
−1/2
Σ1
..
. −1/2 Σn
.
(19)
¯i vi = w ¯ i denotes ¯ i , where X It should be reminded that X the transformation matrix that corresponds to µi . The corresponding parameter, γ, controls the weight of the prior knowledge on the location of Pi . Each point on the template, S, becomes more difficult to move away from its average location, when γ is larger. The quadratic form of the cost function, E0 (X) in (12) is represented as follows:
2
W 0u W 0 D0
√
√
√ α · G0 x − √α · GµE (20) E0 (x) =
β · D
β · uL L
√ √
γ·Λ γ · Λµ
= kA0 x − bk2 . One can directly minimize E0 (x) in (20) by computing ˆ = (A0T A0 )−1 A0T b. Converting x ˆ to a set of the transforx ˆ i }, we obtain the optimal parameters of mation matrix {X the deformation under the condition that the corresponding points, {ui }, are fixed. The outline of the registration algorithm is same with that of the NICP described in sec.1.1. There are two loops in the process: In the inner loop, the corresponding points, ˆ i , are alternatively ui , and the transformation matrices, X updated until the matrices are converged. In the outer loop, not only the stiffness parameter, α, but also γ are successively decreased to αK and γ K , respectively. When γ is large, each point prefers to be located around its average location, which often helps to find appropriate corresponding points in the early stage. The NICP that uses the SSM is abbreviated as NICP-SSM in the followings. ¯i. • Initialize Xi0 = X • For each pair of the parameters (αj , γ j ) ∈ {(α1 , γ 1 ), (α2 , γ 2 ), · · · , (αK , γ K )}, where αj−1 > αj and γ j−1 > γ j • Until kˆ ˆ j−1 k < . xj − x ˆ j−1 v i . • Find corresponding point uj−1 for each X i i • Compute x ˆ j that minimizes E0 (x) in (20).
2.2. Robust Nonrigid ICP In the measurements, the points on the target surface could be perturbed from their true locations by measurement noises. We assume that the perturbation obeys the Gaussian distribution. The problem here is that the combination of the SSM and the model of the measurement noise cannot always describe all points of the given target. For example, measurements of the target surface would include some false points, which do not come from the true target surface. It can also happen that some parts of the target surface are missed by the measurements. As mentioned in the previous section, the closest point, ui , in the target is selected for each template point, Pi . The point in the target, ui , is called outlier when it is not described by the combination of the SSM and the measurement noise model. When there exist outliers, the performance of registration is degraded, hence it is required to automatically detect outliers and exclude them for robust registration. In the proposed method, the outliers are detected based on the assumption that the number of outliers are small compared to the total number of the template points, n. Let Ω denote a set of template points, Pi , which correspond to outliers. The relation between ui and Xi v i is given by { Xi v i + ri di , Pi ∈ S\Ω, ui = (21) ρi , Pi ∈ Ω,
where ri is the size of the deviation between ui and Xi v i , di denotes a unit three-vector defined as di ≡ (ui − Xi v i )/kui − Xi v i k and ρi denotes the location of ui that corresponds to Pi ∈ Ω. As mentioned, we assume that the distribution of ri obeys the Gaussian: p(ri ) = N (ri |0, λ). On the other hand, the distribution of ρi cannot be defined. Let ei di ≡ ui − Xi v i for all points in the template, S. Then, we obtain { e1i , Pi ∈ Ω, ei = kui − Xi v i k = (22) e2i , Pi ∈ S\Ω, where e2i = ri obeys the Gaussian if Pi ∈ S\Ω, and ρi = ui = Xi v i + e1i di ,
(23)
if Pi ∈ Ω. The distribution of e1i is unknown. We set e1i = 0 in S\Ω and e2i = 0 in Ω to rewrite Eq. (21) as follows ui = Xi v i + (e1i + e2i )di .
(24)
Let eq ≡ [eq1 , eq2 , · · · , eqn ]T (q = 1, 2). It should be noted that e1 is large but sparse, while e2 is small but dense. The penalty, φdev , corresponding to the deviation is a sum of the data fidelity term that penalizes the number of nonzero elements in e1 and the negative log-likelihood of e2 : φdev (X, e1 ) ≡ ke1 kL0 +
1 2 2 ke k λ
(25)
1∑ kui − Xi v i − e1i di k2 + ke1 kL0 . λ i=1
Using x, we rewrite Eq. (25) as φ0dev (x, e1 ) = ku − D0 x − Je1 k2 + λke1 kL0 , where J is a 3n × n matrix: d1 d2 J =
(26)
..
.
.
(27)
dn ¯d (x) in (12), we define the cost Using φ0dev instead of E 1 function, E1 (x, e ) for robust registration as follows: E1 (x, e1 )
= ku − D0 x − Je1 k2 + λke1 kL0 ¯s (x) + βEl (x) + γEp (x). (28) + αE
It is NP-hard to minimize this cost function, E1 , hence, as customary, we relax the term of L0-norm in (28) to L1 norm of e1 , such that ¯1 (x, e1 ) E
This is a standard convex relaxation of (28). When all corresponding points, ui , are fixed, the globally optimal solution of (29) can be reached efficiently[15]. The proposed algorithm of the robust NICP is abbreviated as RNICP. The algorithm of RNICP is almost same with NICP-SSM shown in the previous section. The only difference is that RNICP minimizes (29) instead of (20).
3. Experimental Results In the first experiments, using a set of range images, we compared the performance of three methods qualitatively: the nonrigid ICP (NICP), the nonrigid ICP with the SSM (NICP-SSM), and the robust nonrigid ICP (RNICP). Then, we applied these methods for segmenting the livers in given clinical X-ray CT images, and compared the accuracy quantitatively. In the experiments, we did not use any landmark. The El terms in (3), (12), and (29) were omitted in the experiments reported in this section.
3.1. Simulation Experiment
n
=
Figure 1. Left: Regions deformed when the training surfaces were generated. Others: Examples of training surfaces.
= ku − D0 x − Je1 k2 + λke1 kL1 ¯s (x) + βEl (x) + γEp (x). (29) + αE
Deforming the Stanford bunny by using the thin plate spline (TPS), we constructed a set of training data and a set of test data. A set of the control points of the splines was put around the bunny. Let the control point be denoted by ci (i = 1, 2, · · · ), and the set of all of the control points be denoted by C = {ci }. For generating the training surfaces, the original bunny was deformed by moving the positions of some control points randomly. The control points were selected randomly from a set of specific predetermined control points. Let the set of the predetermined control points be denoted by CT ⊂ C. The number of the control points was ten (|CT | = 10) in this experiment. The regions that can be deformed by the control points in CT are shown in Fig.1. We generated m = 1, 000 training surfaces. Some of the examples are also shown in Fig.1. The test data (targets) were generated as follows. The control points that were included in CT were firstly moved randomly and independently to deform the original bunny. Then, dents, holes, or/and false points were added to the deformed surfaces. The dents were also made by applying TPS to the surface, and the moved control point was not included in cj ∈ / CT . The top column of Fig 2 shows some examples of the targets. In the figure, red arrows indicate the places deformed by moving cj ∈ CT . The blue arrows
indicate the outliers: the dents, holes, and false points. The SSM was constructed from the training surfaces and the template was registered to the targets. Three registration algorithms, NICP, NICP-SSM, and RNICP, were applied to each target. We used the identical SSM in NICP-SSM and RNICP, which was constructed from the set of the training surfaces. The series of (αi , γ i ) was also common to those two algorithms. The different series of {αi } was used in NICP. We experimentally found that the value of αK , which controlled the final stiffness of the template, had a larger influence on the registration results. Hence, αK was tuned up more carefully. The initial shape of the template was set to the mean shape for all of the three algorithms. Figure 2 shows examples of experimental results. The top row shows the given target data. The results obtained by NICP, NICP-SSM, and RNICP are shown in the second, the third, and the bottom row, respectively. In this section, the red color density indicates the size of the residual between the resultant template and the inlier target. When e1i 6= 0, then it means that the point is judged as the outlier and |e1i | is indicated by the blue color. The column (A) shows the results obtained from the target that has two dents: one is an inlier dent (indicated by the red arrow) and the other is an outlier one (indicated by the blue arrow). The former dent was generated by moving cj ∈ CT and the latter was generated by moving cj ∈ / CT . NICP deformed the template to fit both of the two dents. NICP-SSM and RNICP, on the other hand, deformed the template to fit only the inlier dent (indicated by the red arrow). At the outlier dent, the shape of the original Stanford bunny was recovered, as shown in the figure. As the result, the residual between the registered template and the target became larger than that obtained by NICP. When such the recovery of the original shape is required as in the case of the anatomical structure segmentation, NICP-SSM outperforms NICP. As shown in the figure, RNICP explicitly detected the outlier region: e1i became nonzero only at the outlier dent. The target shown in the column (B) was made by putting some false points in the target data shown in the column (A). The false points were located above the outlier dent. NICP deformed the template to fit the false points, and the resultant shape changed from that obtained when the false points were not put. NICP-SSM and RNICP recovered the original shape more correctly. Figure3 shows the closeups of the registered templates even when the false points existed. The given target data are shown in Fig.3 (A). The original shape before the dent was generated is shown in (B). The results obtained by NICP, NICP-SSM and RNCIP are shown in (C), (D), and (E), respectively. The largest distance between the original bunny and the recovered template was 0.24, 0.03 and 0.009 when NICP, NICP-SSM, and RNCIP were used, respectively. RNICP recovered the orig-
(A)
(B)
(C)
Given Target
NICP
NICP-SSM
RNICP Figure 2. Examples of registration results. Top row: Original target data. Registration results obtained from the original target data are shown in each column. Red color density indicates the residual size between the corresponding points in inlier regions. Blue color indicates |e1i | in outlier regions.
inal shape more accurately. In addition, RNICP detected the false points as outliers. These results would change if one changes the value of the stiffness parameter αK . One of the advantages of SSM is that the tune-up of the parameter is much easier: The stiffness of the template is homogeneous in NICP, but is inhomogeneous in NICP-SSM and in RNICP reflecting how deformable each position of the template is. The column (C) of Fig.2 shows the results obtained from the target data, in which the left ear of the bunny was cut. NICP and NICP-SSM covered the hole of the left ear with flat patch, as shown in the figures. NICP-SSM could not recover the shape of the ear because the residual would have been too large if the regular shape of the ear had been recovered. RNICP recovered the shape of the ear. It was because the template points in the recovered region were judged as outliers and were automatically excluded from the sum of squared residuals: The shape of the recovered region was ¯s and Ep , determined based on the regularization terms, E in (29). Figure 4 shows the closeups of the cut ear. The
(A) Target
(B) Original
Figure 5. Left: An example of the X-CT image. Right: The corresponding target data. The nuisances are indicated by the blue arrows.
(C) NICP
(D) NICP-SSM
RNICP
Figure 3. Closeups of the registered templates shown in Fig.2 (B)
(A) Target
(C) NICP-SSM
(A) Gold standard
(B) NICP
(C) NICP-SSM
(D) RNICP
(B) NICP
(D) RNICP
Figure 4. Closeups of the registered templates shown in Fig.2 (C). The pale blue shows the shape of the original ear.
largest distance between the registered template and the ear that was cut was 2.22, 1.83, and 0.013 when NICP, NICPSSM, and RNCIP were used, respectively.
3.2. Clinical 3D Image Experiment We applied the three algorithms for registering the template of the liver to enhanced CT images of abdomen. The image size was 512 × 512 × 200, and each voxel size was 0.625 × 0.625 × 1 mm3 . Given 36 CT images, we normalized the size and the shape of the body in each image based on the costal bones. For the leave-one-out cross validation, we constructed the SSM by using 35 of the given images, and registered the resultant template to the one image unused for the SSM construction. In order to construct the training surfaces of the liver, we manually extracted the surface of the liver in each image and generated a set of n = 500 corresponding points, Pik , on each of the surfaces. There are many methods for generating corresponding points on given surfaces[10]. We employed a method
Figure 6. The templates registered to the target data shown in Fig.5
based on the GMDS[4]. For obtaining target data, we applied an edge detector to the test image. Figure 5 shows a slice image of the target. There are tumors in the liver. An example of the target data is shown in Fig. 5. The surface of the liver has two large dents corresponding to the tumors (indicated by the arrow (A) and (B)). In addition, many false points were detected near the target surface. The lug indicated by the arrow (C) were generated by a bone contacting with the liver. Figure 6 shows the results of the registration. NICP deformed the template to fit the dents generated by the tumors and to false points extracted from the surfaces of the anatomical structures near the liver, as shown in Fig.6 (B). NICP-SSM and RNICP interpolated the dents of the tumors and excluded the false points successfully, as shown in Fig.6 (C) and (D), respectively. The latter one explicitly detected the two tumors and the false lug as shown in Fig.6 (D). Table 1 shows the average residuals between the registered template and the gold standard surface extracted manually. RNICP registered the template more accurate than NICP and NICP-SSM. The resultant locations of some of the template points are
Table 1. Comparison of the residuals
Averaged residual [mm]
NICP 2.61
NICP-SSM 2.06
RNICP 1.94
http://graphics.stanford.edu/data/3Dscanrep. This work was funded by Grant-in-Aid for Scientific Research on Innovative Areas: Computational Anatomy for Computeraided Diagnosis and Therapy, MEXT, Japan.
References
(A) NICP
(B) NICP-SSM
(C) RNICP
Figure 7. Locations of the registered template points. The red dots indicate the inlier ones and the blue ones indicate the outliers.
demonstrated in Fig.7. Only the template points that located in the slice images are shown in the figure. The template was deformed largely by NICP at the tumor, as shown in Fig.7 (A). NICP would have interpolated the large dent if the template had been more stiff. We used the soft template, though, because the registration accuracy was better when the template was soft. NICP-SSM located many template points at the boundary of the dent as shown in Fig.7 (B), while RNICP located some outlier points at the lid of the dent as shown in (C).
4. Conclusion and Future Work We proposed the robust nonrigid ICP algorithm (RNICP). The statistical shape model (SSM) of the template is incorporated to the nonrigid ICP (NICP), and outliers are detected based on their sparseness. The shape of the template is parametrized by the transformation matrices of the template points, and the optimal matrices that minimize the cost function can be computed directly, under the condition that the corresponding points on the target are fixed. The algorithm updates the transformation matrices and the corresponding points alternately, while the stiffness coefficients are decreased. Experimental results showed that RNICP improved the registration accuracy by automatically detecting the outliers. Some false targets, though, could not be detected as outliers, when their local shapes were included in the representation of the SSM. It has reported that the soft assignment[11] also improves the robustness against the nuisances. It is included in the future works to improve the way to assign corresponding points to the template points.
5. Acknowledgement The bunny Stanford 3D
model is Scanning
courtesy of Repository
the at
[1] B. Amberg, S. Romdhani, and T. Vetter. Optimal step nonrigid icp algorithms for surface registration. In Computer Vision and Pattern Recognition, 2007. CVPR’07. IEEE Conference on, pages 1–8. IEEE, 2007. 1, 2, 3 [2] A. Ayvaci, M. Raptis, and S. Soatto. Occlusion detection and motion estimation with convex optimization. Advances in Neural Information Processing Systems, 2010. 2 [3] P. Besl and N. McKay. A method for registration of 3-d shapes. IEEE Transactions on pattern analysis and machine intelligence, 14(2):239–256, 1992. 2 [4] A. Bronstein, M. Bronstein, and R. Kimmel. Generalized multidimensional scaling: a framework for isometryinvariant partial surface matching. Proceedings of the National Academy of Sciences of the United States of America, 103(5):1168, 2006. 7 [5] R. Campbell and P. Flynn. A survey of free-form object representation and recognition techniques. Computer Vision and Image Understanding, 81(2):166–210, 2001. 1 [6] T. Cootes, A. Hill, C. Taylor, and J. Haslam. Use of active shape models for locating structures in medical images. Image and vision computing, 12(6):355–365, 1994. 1 [7] T. Cootes, C. Taylor, D. Cooper, J. Graham, et al. Active shape models-their training and application. Computer vision and image understanding, 61(1):38–59, 1995. 1 [8] R. Davies, C. Twining, T. Cootes, and C. Taylor. Building 3-d statistical shape models by direct optimization. Medical Imaging, IEEE Transactions on, 29(4):961–981, 2010. 2 [9] N. Duta and M. Sonka. Segmentation and interpretation of mr brain images. an improved active shape model. Medical Imaging, IEEE Transactions on, 17(6):1049–1062, 1998. 1 [10] T. Heimann and H. Meinzer. Statistical shape models for 3d medical image segmentation: A review. Medical image analysis, 13(4):543–563, 2009. 1, 7 [11] R. Horaud, F. Forbes, M. Yguel, G. Dewaele, and J. Zhang. Rigid and articulated point registration with expectation conditional maximization. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 33(3):587 – 602, 2009. 2, 8 [12] K. Lekadir, R. Merrifield, and G. Yang. Outlier detection and handling for robust 3-d active shape models search. Medical Imaging, IEEE Transactions on, 26(2):212–222, 2007. 1 [13] H. Li and O. Chutatape. Automated feature extraction in color retinal images by a model based approach. Biomedical Engineering, IEEE Transactions on, 51(2):246–254, 2004. 2 [14] J. Phillips, R. Liu, and C. Tomasi. Outlier robust icp for minimizing fractional rmsd. In 3-D Digital Imaging and Modeling, 2007. 3DIM’07. Sixth International Conference on, pages 427–434. IEEE, 2007. 2 [15] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288, 1996. 5
Takamiti Matsuno Nagoya Institute of Technology Aichi, JAPAN
Yoshihide Sawada
[email protected]
Abstract
is registered to target data, which are extracted from the images. The problem here is the existence of nuisances: False target points would be extracted from anatomical structures that are not targeted. In addition, some local but large deformation would occur because of lesions near the target surface, and some true target points would be missed because of, e.g., the contact with neighboring organs. In such situations, it is known that statistical shape model (SSM) is required for registering the template robustly[10]. The active shape model (ASM) is widely used for the registration[6][7]. In the ASM, the template surface is represented by a point distribution model (PDM), which represents a surface by using a set of n points on it. Let wi denote the 3D coordinates of the point Pi (i = 1, 2, · · · , n), and let w ≡ [wT1 , wT2 , · · · , wTn ]T denote the global shape of the template. Given a set of m training surfaces, {w1 , w2 , · · · , wm }, one constructs a linear shape model of w for ASM by applying PCA to the covariance matrix of the training surfaces:
We show how to incorporate a statistical shape model into the nonrigid ICP framework, and propose a robust nonrigid ICP algorithm. In the nonrigid ICP framework, a template surface is represented by a set of points, and the shape of the template is parametrized by a transformation matrix per one template point. In the proposed method, the statistics of the matrices are estimated based on a set of training surfaces, and the statistical shape model is incorporated into the nonrigid ICP framework by modifying the representation of the stiffness of the template. The statistical shape model and a noise model make it possible to discriminate outliers from inliers in given targets. Our proposed method detects the outliers, which are not represented by the models appropriately, based on their sparseness. The detected outliers are automatically excluded from the target to be registered, and the template is deformed to fit the inliers only. As the result, the accuracy of the registration is improved. The performance of the proposed method is evaluated qualitatively and quantitatively using synthetic data and clinical CT images.
¯ + P θ, w=w
(1)
¯ is the mean of the training surfaces, P is a 3n × k where w (k < 3n) matrix each of which column is the eigen vector of the covariance matrix, and θ is a k-vector of the deformation parameters. The ASM iteratively updates the estimate of θ that minimizes the sum of the squared residuals between the point Pi on the template and its corresponding point on a given target. It is known that the least square approach is optimal only when the residuals obey the Gaussian distribution. This Gaussian model, though, does not always describe the distribution of the residuals appropriately, because of the nuisances: They often become the outliers of the Gaussian model. Many methods, hence, are proposed for making the ASM robust: Outliers and inliers are automatically discriminated, and the template is deformed to fit the inlier targets in ways of consistent with the SSM[10]. Duta and Sonka[9], for example, proposed a robust ASM that detects outliers based on the variance information from the PDM. Lekadir et al.[12] proposed another robust method that de-
1. Introduction The registration of 3D surfaces is an important task in computer vision. It has broad applications, such as segmentation of anatomical structures in 3D medical images, face morphing, 3D shape retrieval and object recognition[5]. Registering a model surface (template) to a given measured surface (target) means finding a mapping between the template surface and the target surface that describes the position of semantically corresponding points[1]. Numerous algorithms for the model-to-data registration have been proposed, and it is difficult to focus on some specific drawbacks and advantages of those algorithms without limiting their applications. The main application of the proposed method is the segmentation of the anatomical structures in given 3D medical images: The surface template of each anatomical structure 1
tects outliers based on the ratio of landmark distances. The techniques of the statistical sampling or of the robust point matching are also used for making ASM more robust[13]. The iterative closest point (ICP)[3] is also the algorithm that registers the PDM to given measured targets, and many modifications to the original rigid ICP are proposed to improve the robustness against outliers. Phillips et al.[14], for example, proposed the Outlier Robust ICP, which registers the template by using a subset of the template points with the smallest residual distances. One of the key points of this method is that the fraction, f ∈ [0, 1], of the subset can be adaptively updated in the registration process. Horaud et al.[11], proposed another robust ICP, which can explicitly detect outliers by means of the EM algorithm. A latent variable is introduced for each template point in order to describe if the point corresponds to the outlier or not, so that the template is registered to the inlier target points while the outliers are excluded. These are, though, the rigid registration method and no SSM is used. The nonrigid ICP (NICP)[1] is the another modification to the original rigid ICP. The deformation of the template is represented by an affine transformation matrix per one vertex, and the template is registered by estimating the optimal matrices that minimizes a cost function. The details will be described below. One of the advantages of the NICP is that the optimal transformation matrices can be computed directly when the corresponding points on the target are fixed, because the cost function is a quadratic system of the unknown matrices. The cost function includes a regularization term that requires the template to be stiff, but no SSM is used for restricting the shape of the template. If the aim is to register the template to anatomical structures, the specification of the template[8] is not enough when its shape is restricted only by the stiffness. In this paper, we incorporate the SSM of the template surface to the NICP. The SSM does not directly represent the statistics of the coordinates of the template points but of the statistics of the transformation matrices. Given the SSM and a noise model, we propose a robust NICP algorithm, which detects outliers based on their sparseness and which registers the template only to the inliers. This framework can be found in the optical flow computation proposed by Ayvaci et al.[2]: The occlusion detection and the optical flow computation are simultaneously realized based on the sparse distribution of the occlusions, at which the optical flow constraint is not satisfied. Before describing our proposition, the outline of the NICP is explained in the next subsection.
1.1. Nonrigid ICP Algorithm The template S = (V, E) is given as a set of n vertices, V, and a set of edges, E. Let the homogeneous coordinates of the vertex i be denoted by v i = [xi , yi , zi , 1]T
(i = 1, 2, · · · , n). The deformation of the template is parametrized by one affine 3 × 4 transformation matrix, Xi , per template vertex, v i . The unknowns are organized in a 4n × 3 matrix: X ≡ [X1 X2 · · · Xn ]T .
(2)
The template is registered to the given target surface, T , by minimizing the following cost function: E(X) ≡ Ed (X) + αEs (X) + βEl (X).
(3)
Here, Ed (X) is a sum of squared distances between the target surface, T , and the template vertices, Ti v i , such that ∑ Ed (X) ≡ wi dist2 (T , Xi v i ), (4) v i ∈V where wi is a positive weight coefficient and dist(T , Xi v i ) is the distance between Xi v i and its closest point on T . Es (X) in (3) is a stiffness term, which regularizes the deformation. The weighted difference of the transformations of neighbouring vertices is penalized as follows: ∑ Es (X) ≡ k(Xi − Xj )Gk2F , (5) (i,j)∈E
where k · kF denotes the Frobenius norm and G is a weighting 4×4 matrix, which is needed for balancing the scales of the elements of Xi : G ≡ diag(1, 1, 1, s). The coefficient, α, in (3) controls the stiffness of the template S, which becomes more stiff when α is larger. El (X) in (3) is a landmark term, which is used for the initialization and the guidance of the template. Let a set of landmarks be denoted by L ≡ {(v i1 , l1 ), (v i2 , l2 ), · · · , (v il , ll )}. The distance between the landmark v i on S and the corresponding one ls on T is penalized as follows: ∑ kXi v i − li k2 . (6) El (X) ≡ (v i ,li )∈L When the closest points on T are fixed, one can directly compute the exact minimizer of the cost function E(X) in (3). Let ui denote the point on T that is closest to Xi v i . Then, Ed (X) in (3) becomes ∑ ¯d (X) ≡ E wi kXi v i − ui k2 (7) v i ∈V 2
= kW (DX − U )kF , √ √ where W ≡ diag( w1 , · · · , wn ), D is an n × 4n sparse matrix: T v1 v T2 D≡ (8) , .. . v Tn
and U ≡ [u1 u2 · · · un ]T . It should be noted that Eq. (7) is a sparse quadratic system of the unknown parameters ¯d , Es and El in (3) can be represented in X. Similar to E quadratic forms. The stiffness term, Es , becomes as follows: Es (X) = k(M ⊗ G)Xk2F , (9) where M is the node-arc matrix. When two vertices, i and j, are connected by an edge, r, then the nonzero entries of M in row r are (M )ri = −1 and (M )rj = 1. The landmark term, El , becomes El (X) = kDL X − UL k2F ,
(10)
where DL is a matrix obtained by taking the rows of D corresponding to the landmark vertices, and UL = [l1 , · · · , ll ]T . As the result, the cost function E(X) in (3) can be represented by a quadratic function of X:
2
WU
√ WD
¯ α · M ⊗ G X − 0 E(X) =
√ √
β · DL β · U L F (11) = kAX − Bk2F . ¯ One can directly minimize E(X) in (11) by computing T −1 T ˆ X = (A A) A B. The NICP algorithm can be described as follows: • Initialize X 0 . • For each αj ∈ {α1 , α2 , · · · , αK }, αj−1 > αj . ˆj − X ˆ j−1 kF < . • Until kX ˆ j−1 v i . • Find corresponding point uj−1 for each X i i j ˆ that minimizes E(X). ¯ • Compute X In the registration process, the stiffness coefficient parameter, α, is controlled: Starting with a more stiff model, the global alignment is recovered firstly, and then successively lower stiffness coefficient is used[1].
2. Method As mentioned above, a statistical model of the template shape is incorporated to the NICP method in this article. The statistical model is constructed based on a set of training surfaces. This model is useful for distinguishing true target points from false ones in given data. In addition, the model provides a criteria for detecting the outliers.
2.1. Nonrigid ICP with Statistical Shape Model Same with the NICP, the proposed method represents the template surface by using a set of n points, {Pi |i = 1, 2, · · · , n}, and its deformation is parametrized by one affine 3 × 4 transformation matrix, Xi , per template vertex. We firstly explain how to construct the probabilistic model of {Xi } from a given set of training surfaces, and show the registration algorithm. Then, we propose the robust registration method, which automatically detects outliers based on their sparsity.
2.1.1
Construction of Statistical Shape Model
Let S k (k = 1, 2, · · · , m) denote the kth training surface. For constructing the probabilistic model of {Xi }, the average shape of {S k } is firstly computed, and then each of the training surfaces is registered to the average shape by means of the NICP. Let Xik denote the resultant matrix obtained by the NICP. Then, the probabilistic model of Xi is constructed based on {Xik |k = 1, 2, · · · , m}. ¯ a set of n correFor computing the average shape, S, k sponding points, {Pi |i = 1, 2, · · · , n}, is generated on each surface, S k . Let wki denote the 3D coordinates of ¯ is the polygonal surface, Pik . The average surface, S, ¯ whose vertices, Pi , are determined by computing the mean ¯i = coordinates for ∑ each of the corresponding points: w ∑ k k ¯ k w i /m = ( k Xi /m)v i = Xi v i . Then, registering S k to S¯ by means of the NICP, we obtain a set of transformation matrices, {Xik |i = 1, 2, · · · , n}. Let the transformation matrix, Xi , be vectorized to a 12vector, xi by vertically concatenating the columns of XiT , and let x ≡ [xT1 , xT2 , · · · , xTn ]T . In the proposed method, two probability distributions, p(xi ) and p(xi −xj ), are estimated based on the transformation matrices, {Xik }. We represent these two distributions by using the 12-dimensional multivariate Gaussian distributions, N (·|µ, Σ), where µ denotes the mean and Σ denotes the covariance. Let p(xi ) ∼ N (xi |µi , Σi ) and let p(xi − xj ) = N (xi − xj |µij , Σij ). These parameters, µi , µij , Σi , and Σij are estimated based on {Xik }. Here, p(xi − xj ) is estimated for the pair of Pi and Pj , if there exists an edge between them in the polygo¯ nal surface, S.
2.1.2
Optimal Step for Nonrigid ICP with Statistical Shape Model
The probabilistic models, p(xi ) and p(xi − xj ), are incor¯ porated to the NICP by modifying the cost function E(X) in (11) as follows. ¯d (X) + αE ¯s (X) + βEl (X) + γEp (X). (12) E0 (x) = E ¯d (X) Here, the first and the third terms of the right side, E and El (X), are shown in (7) and (6), respectively. Using x ¯d (X) as instead of X, one can rewrite E ¯d (x) = kW 0 (D0 x − u)k2 , E
(13)
where W 0 is a 3n × 3n matrix and is defined as √ W0 ≡
w1 · I3 ..
.
√
, wn · I3
(14)
in which Id denotes the d × d identity matrix. D0 in (13) denotes a 3n × 12n matrix and is defined as T v 1 ⊗ I3 .. (15) D0 ≡ . . v Tn ⊗ I3 u in (13) denotes a 3n-vector: u = [uT1 , uT2 , · · · , uTn ]T . ¯d , El (X) in (12) can be written as follows: Analogous to E 0 El (x) = kDL x − uL k2 ,
(16)
0 where DL and uL are obtained by taking the elements that corresponds to the landmarks from D0 and u in (13), respectively. ¯s (X) in (12), is a modified version of The second term, E the stiffness term, which penalizes local deformations that are statistically rare, and is defined as follows: ∑ ¯s (X) ≡ E (xi − xj − µij )T Σ−1 ij (xi − xj − µij ) (i,j)∈E
=
∑
kGij x − Gij µk2
(i,j)∈E
= kGx − Gµk2 .
(17) −1/2
⊗ Mij , Here, µ ≡ [µT1 , µT2 , · · · , µTn ]T and Gij ≡ Σij where Mij is an n-dimensional row vector, in which ith element is one, jth element is minus one, and the others are zero. Concatenating all Gij vertically, we obtain a 12|E| × 12n matrix G in (17). The fourth term, Ep (X) in (12), is introduced for penalizing the distance between the ¯ i , as transformed location, Xi v i , and the mean location, w follows. ∑ Ep (X) ≡ (xi − µi )T Σ−1 i (xi − µi ) i
= kΛx − Λµk2 , where
Λ≡
(18)
−1/2
Σ1
..
. −1/2 Σn
.
(19)
¯i vi = w ¯ i denotes ¯ i , where X It should be reminded that X the transformation matrix that corresponds to µi . The corresponding parameter, γ, controls the weight of the prior knowledge on the location of Pi . Each point on the template, S, becomes more difficult to move away from its average location, when γ is larger. The quadratic form of the cost function, E0 (X) in (12) is represented as follows:
2
W 0u W 0 D0
√
√
√ α · G0 x − √α · GµE (20) E0 (x) =
β · D
β · uL L
√ √
γ·Λ γ · Λµ
= kA0 x − bk2 . One can directly minimize E0 (x) in (20) by computing ˆ = (A0T A0 )−1 A0T b. Converting x ˆ to a set of the transforx ˆ i }, we obtain the optimal parameters of mation matrix {X the deformation under the condition that the corresponding points, {ui }, are fixed. The outline of the registration algorithm is same with that of the NICP described in sec.1.1. There are two loops in the process: In the inner loop, the corresponding points, ˆ i , are alternatively ui , and the transformation matrices, X updated until the matrices are converged. In the outer loop, not only the stiffness parameter, α, but also γ are successively decreased to αK and γ K , respectively. When γ is large, each point prefers to be located around its average location, which often helps to find appropriate corresponding points in the early stage. The NICP that uses the SSM is abbreviated as NICP-SSM in the followings. ¯i. • Initialize Xi0 = X • For each pair of the parameters (αj , γ j ) ∈ {(α1 , γ 1 ), (α2 , γ 2 ), · · · , (αK , γ K )}, where αj−1 > αj and γ j−1 > γ j • Until kˆ ˆ j−1 k < . xj − x ˆ j−1 v i . • Find corresponding point uj−1 for each X i i • Compute x ˆ j that minimizes E0 (x) in (20).
2.2. Robust Nonrigid ICP In the measurements, the points on the target surface could be perturbed from their true locations by measurement noises. We assume that the perturbation obeys the Gaussian distribution. The problem here is that the combination of the SSM and the model of the measurement noise cannot always describe all points of the given target. For example, measurements of the target surface would include some false points, which do not come from the true target surface. It can also happen that some parts of the target surface are missed by the measurements. As mentioned in the previous section, the closest point, ui , in the target is selected for each template point, Pi . The point in the target, ui , is called outlier when it is not described by the combination of the SSM and the measurement noise model. When there exist outliers, the performance of registration is degraded, hence it is required to automatically detect outliers and exclude them for robust registration. In the proposed method, the outliers are detected based on the assumption that the number of outliers are small compared to the total number of the template points, n. Let Ω denote a set of template points, Pi , which correspond to outliers. The relation between ui and Xi v i is given by { Xi v i + ri di , Pi ∈ S\Ω, ui = (21) ρi , Pi ∈ Ω,
where ri is the size of the deviation between ui and Xi v i , di denotes a unit three-vector defined as di ≡ (ui − Xi v i )/kui − Xi v i k and ρi denotes the location of ui that corresponds to Pi ∈ Ω. As mentioned, we assume that the distribution of ri obeys the Gaussian: p(ri ) = N (ri |0, λ). On the other hand, the distribution of ρi cannot be defined. Let ei di ≡ ui − Xi v i for all points in the template, S. Then, we obtain { e1i , Pi ∈ Ω, ei = kui − Xi v i k = (22) e2i , Pi ∈ S\Ω, where e2i = ri obeys the Gaussian if Pi ∈ S\Ω, and ρi = ui = Xi v i + e1i di ,
(23)
if Pi ∈ Ω. The distribution of e1i is unknown. We set e1i = 0 in S\Ω and e2i = 0 in Ω to rewrite Eq. (21) as follows ui = Xi v i + (e1i + e2i )di .
(24)
Let eq ≡ [eq1 , eq2 , · · · , eqn ]T (q = 1, 2). It should be noted that e1 is large but sparse, while e2 is small but dense. The penalty, φdev , corresponding to the deviation is a sum of the data fidelity term that penalizes the number of nonzero elements in e1 and the negative log-likelihood of e2 : φdev (X, e1 ) ≡ ke1 kL0 +
1 2 2 ke k λ
(25)
1∑ kui − Xi v i − e1i di k2 + ke1 kL0 . λ i=1
Using x, we rewrite Eq. (25) as φ0dev (x, e1 ) = ku − D0 x − Je1 k2 + λke1 kL0 , where J is a 3n × n matrix: d1 d2 J =
(26)
..
.
.
(27)
dn ¯d (x) in (12), we define the cost Using φ0dev instead of E 1 function, E1 (x, e ) for robust registration as follows: E1 (x, e1 )
= ku − D0 x − Je1 k2 + λke1 kL0 ¯s (x) + βEl (x) + γEp (x). (28) + αE
It is NP-hard to minimize this cost function, E1 , hence, as customary, we relax the term of L0-norm in (28) to L1 norm of e1 , such that ¯1 (x, e1 ) E
This is a standard convex relaxation of (28). When all corresponding points, ui , are fixed, the globally optimal solution of (29) can be reached efficiently[15]. The proposed algorithm of the robust NICP is abbreviated as RNICP. The algorithm of RNICP is almost same with NICP-SSM shown in the previous section. The only difference is that RNICP minimizes (29) instead of (20).
3. Experimental Results In the first experiments, using a set of range images, we compared the performance of three methods qualitatively: the nonrigid ICP (NICP), the nonrigid ICP with the SSM (NICP-SSM), and the robust nonrigid ICP (RNICP). Then, we applied these methods for segmenting the livers in given clinical X-ray CT images, and compared the accuracy quantitatively. In the experiments, we did not use any landmark. The El terms in (3), (12), and (29) were omitted in the experiments reported in this section.
3.1. Simulation Experiment
n
=
Figure 1. Left: Regions deformed when the training surfaces were generated. Others: Examples of training surfaces.
= ku − D0 x − Je1 k2 + λke1 kL1 ¯s (x) + βEl (x) + γEp (x). (29) + αE
Deforming the Stanford bunny by using the thin plate spline (TPS), we constructed a set of training data and a set of test data. A set of the control points of the splines was put around the bunny. Let the control point be denoted by ci (i = 1, 2, · · · ), and the set of all of the control points be denoted by C = {ci }. For generating the training surfaces, the original bunny was deformed by moving the positions of some control points randomly. The control points were selected randomly from a set of specific predetermined control points. Let the set of the predetermined control points be denoted by CT ⊂ C. The number of the control points was ten (|CT | = 10) in this experiment. The regions that can be deformed by the control points in CT are shown in Fig.1. We generated m = 1, 000 training surfaces. Some of the examples are also shown in Fig.1. The test data (targets) were generated as follows. The control points that were included in CT were firstly moved randomly and independently to deform the original bunny. Then, dents, holes, or/and false points were added to the deformed surfaces. The dents were also made by applying TPS to the surface, and the moved control point was not included in cj ∈ / CT . The top column of Fig 2 shows some examples of the targets. In the figure, red arrows indicate the places deformed by moving cj ∈ CT . The blue arrows
indicate the outliers: the dents, holes, and false points. The SSM was constructed from the training surfaces and the template was registered to the targets. Three registration algorithms, NICP, NICP-SSM, and RNICP, were applied to each target. We used the identical SSM in NICP-SSM and RNICP, which was constructed from the set of the training surfaces. The series of (αi , γ i ) was also common to those two algorithms. The different series of {αi } was used in NICP. We experimentally found that the value of αK , which controlled the final stiffness of the template, had a larger influence on the registration results. Hence, αK was tuned up more carefully. The initial shape of the template was set to the mean shape for all of the three algorithms. Figure 2 shows examples of experimental results. The top row shows the given target data. The results obtained by NICP, NICP-SSM, and RNICP are shown in the second, the third, and the bottom row, respectively. In this section, the red color density indicates the size of the residual between the resultant template and the inlier target. When e1i 6= 0, then it means that the point is judged as the outlier and |e1i | is indicated by the blue color. The column (A) shows the results obtained from the target that has two dents: one is an inlier dent (indicated by the red arrow) and the other is an outlier one (indicated by the blue arrow). The former dent was generated by moving cj ∈ CT and the latter was generated by moving cj ∈ / CT . NICP deformed the template to fit both of the two dents. NICP-SSM and RNICP, on the other hand, deformed the template to fit only the inlier dent (indicated by the red arrow). At the outlier dent, the shape of the original Stanford bunny was recovered, as shown in the figure. As the result, the residual between the registered template and the target became larger than that obtained by NICP. When such the recovery of the original shape is required as in the case of the anatomical structure segmentation, NICP-SSM outperforms NICP. As shown in the figure, RNICP explicitly detected the outlier region: e1i became nonzero only at the outlier dent. The target shown in the column (B) was made by putting some false points in the target data shown in the column (A). The false points were located above the outlier dent. NICP deformed the template to fit the false points, and the resultant shape changed from that obtained when the false points were not put. NICP-SSM and RNICP recovered the original shape more correctly. Figure3 shows the closeups of the registered templates even when the false points existed. The given target data are shown in Fig.3 (A). The original shape before the dent was generated is shown in (B). The results obtained by NICP, NICP-SSM and RNCIP are shown in (C), (D), and (E), respectively. The largest distance between the original bunny and the recovered template was 0.24, 0.03 and 0.009 when NICP, NICP-SSM, and RNCIP were used, respectively. RNICP recovered the orig-
(A)
(B)
(C)
Given Target
NICP
NICP-SSM
RNICP Figure 2. Examples of registration results. Top row: Original target data. Registration results obtained from the original target data are shown in each column. Red color density indicates the residual size between the corresponding points in inlier regions. Blue color indicates |e1i | in outlier regions.
inal shape more accurately. In addition, RNICP detected the false points as outliers. These results would change if one changes the value of the stiffness parameter αK . One of the advantages of SSM is that the tune-up of the parameter is much easier: The stiffness of the template is homogeneous in NICP, but is inhomogeneous in NICP-SSM and in RNICP reflecting how deformable each position of the template is. The column (C) of Fig.2 shows the results obtained from the target data, in which the left ear of the bunny was cut. NICP and NICP-SSM covered the hole of the left ear with flat patch, as shown in the figures. NICP-SSM could not recover the shape of the ear because the residual would have been too large if the regular shape of the ear had been recovered. RNICP recovered the shape of the ear. It was because the template points in the recovered region were judged as outliers and were automatically excluded from the sum of squared residuals: The shape of the recovered region was ¯s and Ep , determined based on the regularization terms, E in (29). Figure 4 shows the closeups of the cut ear. The
(A) Target
(B) Original
Figure 5. Left: An example of the X-CT image. Right: The corresponding target data. The nuisances are indicated by the blue arrows.
(C) NICP
(D) NICP-SSM
RNICP
Figure 3. Closeups of the registered templates shown in Fig.2 (B)
(A) Target
(C) NICP-SSM
(A) Gold standard
(B) NICP
(C) NICP-SSM
(D) RNICP
(B) NICP
(D) RNICP
Figure 4. Closeups of the registered templates shown in Fig.2 (C). The pale blue shows the shape of the original ear.
largest distance between the registered template and the ear that was cut was 2.22, 1.83, and 0.013 when NICP, NICPSSM, and RNCIP were used, respectively.
3.2. Clinical 3D Image Experiment We applied the three algorithms for registering the template of the liver to enhanced CT images of abdomen. The image size was 512 × 512 × 200, and each voxel size was 0.625 × 0.625 × 1 mm3 . Given 36 CT images, we normalized the size and the shape of the body in each image based on the costal bones. For the leave-one-out cross validation, we constructed the SSM by using 35 of the given images, and registered the resultant template to the one image unused for the SSM construction. In order to construct the training surfaces of the liver, we manually extracted the surface of the liver in each image and generated a set of n = 500 corresponding points, Pik , on each of the surfaces. There are many methods for generating corresponding points on given surfaces[10]. We employed a method
Figure 6. The templates registered to the target data shown in Fig.5
based on the GMDS[4]. For obtaining target data, we applied an edge detector to the test image. Figure 5 shows a slice image of the target. There are tumors in the liver. An example of the target data is shown in Fig. 5. The surface of the liver has two large dents corresponding to the tumors (indicated by the arrow (A) and (B)). In addition, many false points were detected near the target surface. The lug indicated by the arrow (C) were generated by a bone contacting with the liver. Figure 6 shows the results of the registration. NICP deformed the template to fit the dents generated by the tumors and to false points extracted from the surfaces of the anatomical structures near the liver, as shown in Fig.6 (B). NICP-SSM and RNICP interpolated the dents of the tumors and excluded the false points successfully, as shown in Fig.6 (C) and (D), respectively. The latter one explicitly detected the two tumors and the false lug as shown in Fig.6 (D). Table 1 shows the average residuals between the registered template and the gold standard surface extracted manually. RNICP registered the template more accurate than NICP and NICP-SSM. The resultant locations of some of the template points are
Table 1. Comparison of the residuals
Averaged residual [mm]
NICP 2.61
NICP-SSM 2.06
RNICP 1.94
http://graphics.stanford.edu/data/3Dscanrep. This work was funded by Grant-in-Aid for Scientific Research on Innovative Areas: Computational Anatomy for Computeraided Diagnosis and Therapy, MEXT, Japan.
References
(A) NICP
(B) NICP-SSM
(C) RNICP
Figure 7. Locations of the registered template points. The red dots indicate the inlier ones and the blue ones indicate the outliers.
demonstrated in Fig.7. Only the template points that located in the slice images are shown in the figure. The template was deformed largely by NICP at the tumor, as shown in Fig.7 (A). NICP would have interpolated the large dent if the template had been more stiff. We used the soft template, though, because the registration accuracy was better when the template was soft. NICP-SSM located many template points at the boundary of the dent as shown in Fig.7 (B), while RNICP located some outlier points at the lid of the dent as shown in (C).
4. Conclusion and Future Work We proposed the robust nonrigid ICP algorithm (RNICP). The statistical shape model (SSM) of the template is incorporated to the nonrigid ICP (NICP), and outliers are detected based on their sparseness. The shape of the template is parametrized by the transformation matrices of the template points, and the optimal matrices that minimize the cost function can be computed directly, under the condition that the corresponding points on the target are fixed. The algorithm updates the transformation matrices and the corresponding points alternately, while the stiffness coefficients are decreased. Experimental results showed that RNICP improved the registration accuracy by automatically detecting the outliers. Some false targets, though, could not be detected as outliers, when their local shapes were included in the representation of the SSM. It has reported that the soft assignment[11] also improves the robustness against the nuisances. It is included in the future works to improve the way to assign corresponding points to the template points.
5. Acknowledgement The bunny Stanford 3D
model is Scanning
courtesy of Repository
the at
[1] B. Amberg, S. Romdhani, and T. Vetter. Optimal step nonrigid icp algorithms for surface registration. In Computer Vision and Pattern Recognition, 2007. CVPR’07. IEEE Conference on, pages 1–8. IEEE, 2007. 1, 2, 3 [2] A. Ayvaci, M. Raptis, and S. Soatto. Occlusion detection and motion estimation with convex optimization. Advances in Neural Information Processing Systems, 2010. 2 [3] P. Besl and N. McKay. A method for registration of 3-d shapes. IEEE Transactions on pattern analysis and machine intelligence, 14(2):239–256, 1992. 2 [4] A. Bronstein, M. Bronstein, and R. Kimmel. Generalized multidimensional scaling: a framework for isometryinvariant partial surface matching. Proceedings of the National Academy of Sciences of the United States of America, 103(5):1168, 2006. 7 [5] R. Campbell and P. Flynn. A survey of free-form object representation and recognition techniques. Computer Vision and Image Understanding, 81(2):166–210, 2001. 1 [6] T. Cootes, A. Hill, C. Taylor, and J. Haslam. Use of active shape models for locating structures in medical images. Image and vision computing, 12(6):355–365, 1994. 1 [7] T. Cootes, C. Taylor, D. Cooper, J. Graham, et al. Active shape models-their training and application. Computer vision and image understanding, 61(1):38–59, 1995. 1 [8] R. Davies, C. Twining, T. Cootes, and C. Taylor. Building 3-d statistical shape models by direct optimization. Medical Imaging, IEEE Transactions on, 29(4):961–981, 2010. 2 [9] N. Duta and M. Sonka. Segmentation and interpretation of mr brain images. an improved active shape model. Medical Imaging, IEEE Transactions on, 17(6):1049–1062, 1998. 1 [10] T. Heimann and H. Meinzer. Statistical shape models for 3d medical image segmentation: A review. Medical image analysis, 13(4):543–563, 2009. 1, 7 [11] R. Horaud, F. Forbes, M. Yguel, G. Dewaele, and J. Zhang. Rigid and articulated point registration with expectation conditional maximization. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 33(3):587 – 602, 2009. 2, 8 [12] K. Lekadir, R. Merrifield, and G. Yang. Outlier detection and handling for robust 3-d active shape models search. Medical Imaging, IEEE Transactions on, 26(2):212–222, 2007. 1 [13] H. Li and O. Chutatape. Automated feature extraction in color retinal images by a model based approach. Biomedical Engineering, IEEE Transactions on, 51(2):246–254, 2004. 2 [14] J. Phillips, R. Liu, and C. Tomasi. Outlier robust icp for minimizing fractional rmsd. In 3-D Digital Imaging and Modeling, 2007. 3DIM’07. Sixth International Conference on, pages 427–434. IEEE, 2007. 2 [15] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288, 1996. 5
