A similarity measure for nonrigid volume registration using known joint

0 downloads 0 Views 862KB Size Report
bDepartment of Radiology, Osaka University Graduate School of Medicine, Suita, ..... (Osaka University Hospital, Japan), the position of the VOI .... in which d(x; F) represents Free-form deformation (FFD) .... (d) Color-coded vector field obtained using our measure. ..... registration of medical images using a statistical atlas.
Medical Image Analysis 7 (2003) 553–564 www.elsevier.com / locate / media

A similarity measure for nonrigid volume registration using known joint distribution of targeted tissue: Application to dynamic CT data of the liver Jun Masumoto a , *, Yoshinobu Sato a , Masatoshi Hori b , Takamichi Murakami b , Takeshi Johkoh c , Hironobu Nakamura b , Shinichi Tamura a a

Division of Interdisciplinary Image Analysis, Osaka University Graduate School of Medicine, Suita, Osaka 565 -0871, Japan b Department of Radiology, Osaka University Graduate School of Medicine, Suita, Osaka 565 -0871, Japan c School of Allied Health Sciences, Faculty of Medicine, Osaka University, Suita, Osaka 565 -0871, Japan

Abstract A similarity measure for nonrigid volume registration with known joint distribution of a targeted tissue is developed to process tissue slide at the boundaries between the targeted and non-targeted tissues. Pre-segmentation of the targeted tissue is unnecessary. This measure is applied to registering volumes acquired at different time-phases in dynamic CT scans of the liver using contrast materials and can be derived for the case where only the joint distribution of the targeted tissue is known. The similarity measure is formulated as a likelihood by introducing a concept termed ‘exclusivity condition’ and embedded into a cost function for nonrigid registration to be combined with the smoothness term. In addition, a practical method for estimating the joint distribution of the liver from unregistered clinical CT data is described. We demonstrate experimentally that tissue slide is effectively processed by this proposed measure using simulated dynamic CT data generated from a software phantom and clinical CT data of eight patients.  2003 Elsevier B.V. All rights reserved.

1. Introduction Dynamic contrast-enhanced CT scans are effective means for hepatic disease diagnosis and surgical planning. In a dynamic CT scan operation, CT volumes are typically acquired at different time-phases but not necessarily within a single breath-hold. Hence, these volumes are not always registered between different time-phases due to the respiratory movements. Their registration by post-processing is however highly desirable to (1) perform accurate correlation between different time-phases, (2) register more accurately in 3D rendering of the liver, the portal / hepatic veins and tumors enhanced at different phases, and (3) estimate time–density curves at every voxel, which should eventually permit automatic cancer characterization (Carrillo et al., 2000). The present paper addresses non-rigid registration between volumes acquired at different time-phases of dynamic CT scans of the liver. An important issue in hepatic *Corresponding author. E-mail address: [email protected] (J. Masumoto). 1361-8415 / 03 / $ – see front matter  2003 Elsevier B.V. All rights reserved. doi:10.1016 / S1361-8415(03)00036-7

registration is tissue slide that occurs along the boundaries between the liver and other tissues. Tissue slide results in discontinuities in the 3D vector field, designated non-rigid deformation (Rueckert et al., 1999; Lester and Aridge, 1999). Previous attempts to deal with this problem require either pre-segmentation of the hepatic region (Chen et al., 1999) or specification for possible locations of tissue slide before registration (Wang and Staib, 2000). However, segmentation of the liver from CT data is a rather complicated task (Schenk et al., 2000), direct registration between raw CT volumes without segmentation is therefore desirable in clinical practice. To process tissue slide without pre-segmentation, we have developed a similarity measure for volume registration using known joint distribution. In the dynamic CT, the tissue contrast during scans at different time-phases varies with the particular tissue involved. Thus, unlike crosscorrelation measure (Lemieux et al., 1998), a similarity measure should also be able to cope with differences in contrast between volumes to be registered. Although mutual information (Bell and Sejnowski, 1995; Collignon et al., 1995; Maes et al., 1997; Wells et al., 1996) (or the

J. Masumoto et al. / Medical Image Analysis 7 (2003) 553–564

554

entropy correlation coefficient: ECC (Pluim et al., 2000)) is known to be useful as a similarity measure in such a case (Holden et al., 2000; Roche et al., 1999), we employ known joint distribution of a targeted tissue, originally suggested by Leventon et al. (1998) and extended by Chung et al. (2002). In our present study, we utilize the known joint distribution of the targeted tissue only rather than that of the entire volume in order to register solely the targeted tissue, for example the liver, while ignore nontargeted tissues. And we effectively overcome the tissue slide issue that is inevitable in registration of the abdominal domain. A notable feature of our measure is the simultaneous execution of both rough segmentation and registration of the targeted tissue using its known joint distribution, therefore any pre-segmentation or manual interaction is not needed to deal with tissue slide, an improvement from the reported methods (Chen et al., 1999; Schenk et al., 2000). The following of the paper describes the formulation of our similarity measure as a likelihood and its embedding into a nonrigid registration procedure to be combined with the smoothness constraint (Section 2). The effectiveness of the similarity measure for tissue slide is experimentally demonstrated by using simulated dynamic CT data generated from a software phantom and clinical CT data of patients (Section 3). The significance, advantages and limitations of the similarity measure are discussed (Section 4), followed by future directions (Section 5).

Fig. 1. Example of joint distribution.

distribution P(i, j) estimated from training data sets, the maximum log likelihood transformation T ML is defined as T ML 5argmax T

O logsPsI(x), JsT(x)ddd.

(2)

x

In previous work by Leventon et al. (1998), aimed at rigid registration of the brain in two different MR scans, the whole joint distribution with various anatomical structures in the brain is modeled by mixture of Gaussian and Parzen window models and estimated from correctly aligned training data sets.

2.2. Similarity measure as likelihood 2. Methods

2.1. Registration based on joint intensity model We assume that a pair of images represented as I(x) and T J(x) is aligned correctly, where x 5 [x, y, z] . We define the set of anatomical structures G 5 hg1 , g2 , . . . ,gn j, and assume that anatomical structure g (g [ G) appears with the intensity value i in the image I and j in the image J with a joint probability P(i, jug ). We also define P(g ) be the prior probability of g. Fig. 1 illustrates an example of the joint distribution. In the case of an unregistered pair of images I(x) and J(x), which is assumed to be modeled by the joint distribution described above if aligned correctly, given a hypothesis of registration transformation T(x) to align J with I, the likelihood can be calculated by

P P(I(x), J(T(x))) 5 P O P(g ) ? P(I(x), J(T(x))ug ),

P(I, JuT ) 5

x [I

(1)

x [I g [G

where og [G P(g ) 5 1, and we make the assumption that voxels are independent samples from the distribution (and relative positions of voxels are ignored). Given the joint

The present study is aimed at the application of nonrigid registration of the liver at two different time-phases of dynamic contrast-enhanced abdominal CT data. The difficulties of dealing with the liver in CT data as compared to the brain in MR data are summarized as follows. • In formulating a procedure of nonrigid registration, maximization of image similarity, which is defined as a likelihood in the present study, acts as image constraint, while the so-called ‘smoothness constraint’ is typically combined with it to stabilize the estimation of the deformation field. Although the spatial variations of the deformation vector field within the liver can be assumed to be smooth, those at the boundaries between the liver and other tissues are often discontinuous due to tissue slide. The smoothness constraint shows unwanted behaviors when it is applied across the boundaries along which discontinuities occur. • There are large variations of the joint distribution among different patients in contrast-enhanced CT data (Fig. 2) and thus it is difficult to construct a general model for the whole joint distribution with various anatomical structures in the abdomen from training data sets. In addition, correctly aligned training data sets are not easy to prepare due to the difficulty of manual nonrigid registration. In order to overcome the above difficulties, we model the

J. Masumoto et al. / Medical Image Analysis 7 (2003) 553–564

555

Fig. 2. Joint histograms of two different phases of dynamic contrast-enhanced CT data from eight different patients. The joint histograms were calculated from the volume of interest described in Section 2.2.2. Significant variations among patients were observed. The horizontal and vertical axes correspond to the intensity values of pre- and post-contrast images, respectively. Brightness represents frequency.

joint distributions and formulate a similarity measure as a likelihood in the following ways: • The patient-dependent joint distribution of the targeted tissue, i.e. the liver, is estimated from the unregistered data sets of the patient. • The joint distribution of other tissues is basically modeled as a uniform distribution in order not to affect likelihood (i.e. similarity) computation. Based on the joint distribution modeling described above, the likelihood maximization process tries to register within the targeted tissue regions to make the joint histogram resemble the distribution model of the targeted tissue estimated beforehand while it does not force to register for other non-targeted regions since any image constraint for them is not provided. As a result, the unwanted behavior of the smoothness constraint should be avoided because contradicting image constraints in the non-targeted tissue side of the boundaries are not forced. We assume that the set of anatomical structures G consists of two tissues, the targeted tissue, i.e. the liver (L) and non-hepatic tissues (O), where O represents all the tissues except the targeted tissue. In this case, Eq. (1) is rewritten as P(I, J) 5

P (P(L) ? P(I(x), J(T(x))uL)

;(i, j) h(i, j)uPsi, jugsd . e j,

O and OO O Pst, jug d ]]]]] . r, O O Pst, jug d Psi, tugsd t ]]]]] .r Psi, tugkd gk [G

t

s

t

x[I

1 P(O) ? P(I(x), J(T(x))uO)).

in one image to identified boundaries in the other image will become extremely difficult. To overcome the difficulty, we make an assumption described below. We define Psi, jugsd as being exclusive if Psiugsd and Ps jugsd, which are projections of Psi, jugsd onto the i- and j-axes, respectively, have only negligible overlap with Psiugtd and Ps jugtd for all t (where t ± s), respectively. Fig. 3 illustrates the exclusivity condition. Tissue g3 is exclusive if the distributions of all other tissues gt (where t ± 3) do not overlap with the shaded area shown in Fig. 3. Tissue g3 is not regarded as exclusive in Fig. 3(a) because the distributions of other tissues (g4 , g5 , g6 g8 ) are overlapped with the shaded area in Fig. 3(a), while tissue g3 is regarded as exclusive in Fig. 3(b). Mathematically, Psi, jugsd is regarded as exclusive ifPsi, jugsd satisfies the following conditions:

(3)

The similarity measure is derived by defining P(i, juL), P(i, juO), P(L) and P(O) in the above likelihood.

2.2.1. Exclusivity condition We introduce the exclusivity condition in order to ensure that boundaries between the targeted tissue and nontargeted tissues in local regions are identified in both images to be registered. If the boundaries in either of these images are not identified, finding the corresponding parts

gk [G

t

(4)

k

where e and r are sufficiently close to 0 and 1, respectively. We assume that the joint distribution of the targeted tissue, i.e. the liver (L), satisfies the exclusivity condition. Although this condition may appear to be too rigorous to be satisfied in real term, when nonrigid registration is considered, it is unnecessary to satisfy the exclusivity condition in the whole image since the likelihood calculation (joint histogram evaluation) for nonrigid registration is

556

J. Masumoto et al. / Medical Image Analysis 7 (2003) 553–564

performed in each local area rather than in the whole image. Thus, the condition just needs to be satisfied in a local area, which is practically reasonable. Fig. 4 illustrates the local exclusivity condition. Assuming the joint distribution of the whole image shown in Fig. 4(b), the expected joint histograms of local areas specified in Fig. 4(a) are shown in Fig. 4(c)–(g). Although g3 is not exclusive in the joint distribution of the whole image, g3 can be regarded as exclusive in each local area except Fig. 4(d). The exclusivity condition guarantees the intensity distribution of the targeted tissue (liver) regions be wellseparated from those of non-targeted tissues in each single image of I(x) and J(x). This means that image boundaries between the targeted and non-targeted tissue regions are detectable in both I(x) and J(x). The robustness against the violation of the exclusivity exclusive condition is experimentally examined in Section 3. The methods for estimating P(i, juL), P(i, juO), P(L) and P(O) in Eq. (3) are as follows. Firstly, P(i, juL) is estimated from an unregistered pair of images for which nonrigid registration is performed. Secondly, P(i, juO) is derived so that P(i, juL) satisfies the exclusivity condition. Finally, P(L) and P(O) are estimated based on P(i, juL) and P(i, juO).

Fig. 3. Exclusivity condition. (a) Tissue g3 is not regarded as exclusive. (b) Tissue g3 is regarded as exclusive.

2.2.2. Estimation of P( i, ju L) The joint distribution of a targeted tissue P(i, juL) needs to be obtained before the maximum likelihood registration. As described above, there are large variations in P(i, juL) among patient data sets. Thus, a practical method is

Fig. 4. Local exclusivity condition. (a) Original image. (b) Joint distribution in the whole image. The horizontal and vertical axes correspond to the intensity values of image I and image J, respectively. (c)–(g) Joint distribution in local areas: Area 1, Area 2, Area 3, Area 4 and Area 5.

J. Masumoto et al. / Medical Image Analysis 7 (2003) 553–564

necessary for estimating the joint distribution from two unregistered images to be registered. We assume that the joint distribution P(i, juL) is well approximated by the Gaussian function given by ] ] 1 2( 1 / 2)(v 2v) T S 21 (v 2v) P(vuL) 5 ]]] , ]] e œ2p uSu

(5)

] ] ] ] where v 5 (i, j), ]v 5 (i, j), (i, j) are the average values, and S is covariance matrix. A method for estimating the joint distribution of the liver region from unregistered two images at different time-phases in dynamic contrast-enhanced CT scans of the abdomen is established as follows. The field of view (FOV) for abdominal CT scans is usually set based on the position of the spine. We set the volume of interest (VOI) so that it would be mostly occupied by hepatic tissue (Fig. 5(a)). In abdominal CT data acquired in our hospital (Osaka University Hospital, Japan), the position of the VOI can be fixed for each patient because the position of the liver relative to the spine is not greatly different in each ] ] case. We estimate the averages (i, j) and covariance matrix S of the joint probability distribution Psi, juLd of Eq. (3) by analyzing the joint histogram obtained from the VOI of the ] ] two volumes. (i, j) and S are estimated from the histogram region whose center is the mode of the joint histogram and whose horizontal and vertical widths are three times the full width half maximum (FWHM) values of 1D histograms projected onto the i- and j-axes, respectively. Although the two volumes are not registered at this point, it still provide a good approximation of the joint distribution. Figs. 5(b) and (c) show the joint histogram obtained from the VOI and the estimated joint distribution of the liver, respectively.

2.2.3. Modeling of P( i, ju O) The joint distribution of non-targeted tissues P(i, juO) is basically modeled as a uniform distribution in order not to provide any image constraint on the condition that P(i, juL) being exclusive should be satisfied. P(i, juO) should be modeled as uniform except for the shaded area shown in Fig. 3. We model the distribution of the shaded

557

area shown in Fig. 3 based on the projections P(iuL) and P( juL) of P(i, juL) onto the i-axis and j-axis, respectively. This is realized by subtracting back-projections of P(iuL) and P( juL) from the uniform distribution, after normalization of P(iuL) and P( juL). ] Let GsiuLd and Gs juLd be exp h 2 (i 2i9)2 / 2si 9 j and ] 2 ] ] exph 2 ( j 2j9) / 2sj 9 j, where i9 and j9 are average values of the projections of Psi, juLd onto the i-axis and j-axis, respectively, and si 9 and sj 9 are their standard deviations, respectively. Psi, juOd is modeled as 1 Psi, juOd 5 ] h1 2 max hGsiuLd, Gs juLdjj, S

(6)

where S5

O h1 2 maxhGsiuLd, Gs juLdjj.

(7)

i, j

2.2.4. Estimation of P( L) and P( O) The prior probabilities of tissue L and O basically depend on the ratio of their volumes. In nonrigid registration, the prior probabilities should depend on the ratio of their volumes in the local area where the likelihood is calculated. Because it is not easy to precisely estimate this ratio, we use the following empirical method to determine the prior probabilities. Let P(L) and P(O) be a (0 < a < 1) and 1 2 a, respectively. The value of a is not so sensitive to the registration results. When P(L) ? P(i, juL) 4 P(O) ? P(i, juO), however, the registration processes tend to register the target tissue in one image to the peripheral tissue in the other image if the peripheral tissue have voxel values somewhat similar to the target tissue in one image. When P(L) ? P(i, juL) < P(O) ? P(i, juO), the image constraint derived from the known distribution of the target tissue does not become effective. In order to avoid the these situations, we assume that the maximum probability of P(L) ? P(i, juL) is equal to that of P(O) ? P(i, juO). Based on the assumption, we obtain the constraint on a given by max P(L) ? P(i, juL) 5max P(O) ? P(i, juO). i, j

i, j

(8)

From the above constraint, a is derived as

Fig. 5. Estimation of liver distribution. (a) Volume of interest (VOI). (b) Joint histogram in VOI. (c) Estimated P(i, juL) (in Eq. (3)). (d) Estimated P(i, juO) (in Eq. (3)).

J. Masumoto et al. / Medical Image Analysis 7 (2003) 553–564

558

1 a 5 ]]]]]. max P(i, juL) i, j 1 1 ]]]] max P(i, juO)

(9)

? Update the gradient vector =C ? Move up to finer hierarchical grid (F , → F , 11 ) ? until Finest hierarchical grid

i, j

2.3. Embedding similarity measure into nonrigid registration procedure We embed the similarity measure as a likelihood described above into a cost function defined as C(F) 5

O h2C

similarity

f [F

(f ) 1 lCsmooth (f ) j,

(10)

where Csimilarity (f ) denotes the similarity measure term defined as the likelihood described above, Csmooth (f ) denotes the smoothness term, and l is a weight parameter balancing the two terms. F denotes the parameters describing the registration transform given by T(x) 5 x 1 d (x; F),

(11)

in which d(x; F) represents Free-form deformation (FFD) (Rueckert et al., 1999) described by B-splines, F denotes the whole set of the B-spline control points, and f (f [ F) denotes a subset of the control points involved in each local area. Csimilarity (f ) and Csmooth (f ) are the similarity and smoothness terms (Wahba, 1990) in each local area, which are given by Csimilarity (f ) 5

O logsPsI(x), J(x 1 d(x; f ))dd,

(12)

x[V

and

S D S D S D S D S D S D

≠ 2f ]] Csmooth (f ) 5 ≠x 2

2

2

≠ 2f ]] 1 ≠y 2

≠ f 1 2 ]] ≠xy

2

2

≠ 2f ]] 1 ≠z 2

2

≠ f 1 2 ]] ≠xz

2

2

2

≠ f 1 2 ]] ≠yz

2

, (13)

respectively, where V denotes a local area in which f involves. In order to minimize Eq. (10), we use a steepest descend algorithm (Maes et al., 1999) using a hierarchical grid. We successively refine the deformation field by using the deformation field estimated on a coarse grid as an initial field of the estimation procedure on its next finer grid. This coarse–fine procedure is given by ? Initialize the control points F , ( , 5 0) ? repeat ? Calculate the gradient vector of Eq. (10): ≠C(F , ) =C 5 ]]] ≠F , ? while i=Ciu . e do ? Update the control points:

=C F 5 F 1 m ]] i=Ci

where F , denotes the B-spline control point grid at , th hierarchical level. The method proposed by (Forsey and Bartels, 1988) is employed for the operation F , → F , 11 . In our experiments, the hierarchical grid consisted of three levels, and the grid intervals were 42, 21 and 10.5 mm.

3. Experiments

3.1. Simulation experiments 3.1.1. Synthesized data sets and evaluation methods We evaluated the above described method using synthesized CT images generated from a software phantom shown in Fig. 6(a). The software phantom was designed to simulate the dynamic contrast-enhanced CT scans of the abdomen. The CT values of organs and tissues in the phantom were set to be similar to those in real dynamic CT image. Gaussian noise was added to the phantom data to simulate noise in CT image. Parameters related to image intensities such as the CT values and the standard deviation (S.D.) of Gaussian noise at two difference time-phases (pre- and post-contrast) were summarized in Table 1. It should be noted that the CT values of the liver and the vessel were different between two time-phases to simulate dynamic contrast enhancement. The experiments were specially aimed at evaluating effectiveness for tissue slide and robustness against the violation of the exclusivity condition. The CT values shown in D9 of Table 1 were used in the experiments for evaluating the robustness against the violation of the exclusivity condition (shown in Fig. 10) only. The CT value shown in D was used for other experiments (shown in Figs. 7–9). Synthesized CT images used for the experiments are shown in Fig. 6(b). Geometric parameters used to generate the phantom are summarized in Table 2. We compared the proposed similarity measure with the entropy correlation coefficient (ECC), which is equivalent to normalized mutual information (NMI) (Pluim et al., 2000) and Pearson product–moment cross correlation (NCC). For the comparison, the similarity term, Csimilarity (f ), in Eq. (10) was replaced by ECC and NCC. In each of our measure, ECC and NCC, registration accuracy was evaluated for various values of l in Eq. (10). The following two criteria for the evaluations were used: (1) Difference of region shape (Dice, 1945) (the difference between the liver region shape in I(x) and that in J(T(x))): Let VI be the liver region pre-contrast image I(x),

J. Masumoto et al. / Medical Image Analysis 7 (2003) 553–564

559

Fig. 6. Software phantom and synthesized CT images used in simulation experiments. (a) Software phantom (A: air, B: rib, C: fat, D: liver, E: vessel). Its xy-plane (left), xz-plane (middle) and 3D visualization (right) are shown. (b) Synthesized (pre-contrast) CT images generated from software phantom (xy-plane). CT values of the liver region (D) was 80 (left), 60 (middle) and 50 (right). The middle and right images were used for evaluating robustness against violation of the exclusivity condition.

and VJ be that in registered post-contrast image J(T(x)). We defined the difference of region shape Eregion as VI VJ Eregion 5 ]]]]] 3 100 (%). VI