# Joint Segmentation and Groupwise Registration of Cardiac Perfusion Images Using Temporal Information

Authors: Dwarikanath Mahapatra
Journal: ISRN Machine Vision
 (2012)

---

## Abstract

We propose a joint segmentation and groupwise registration method for cardiac perfusion images by using temporal information. The nature of perfusion images makes groupwise registration especially attractive as the temporal information from the entire image sequence can be used. Registration aims to maximize the smoothness of the intensity signal, while segmentation minimizes a pixel’s dissimilarity with other pixels having the same segmentation label. The cost function is optimized in an iterative fashion using B-splines. Tests on real patient datasets show that compared to two other methods, our method shows lower registration error and higher segmentation accuracy. This is attributed to the use of temporal information for groupwise registration and mutually complementary registration and segmentation information in one framework, while other methods solve the two problems separately.


---

## Body

## 1. Introduction

Dynamic contrast-enhanced (DCE) magnetic resonance (MR) images (or perfusion MRI) have developed as a popular noninvasive tool for the functional analysis of internal organs. Contrast agent is injected intravenously into the patient, and a series of MR images are acquired over a period of time. As the contrast agent flows through the blood stream, the intensity of corresponding regions increases. Since the image acquisition process can take up to 20 minutes, patient movement is inevitable. Additionally, elastic deformation of cardiac tissues due to patient breathing needs to be corrected. Perfusion images are characterized by rapid intensity change over time, low spatial resolution, and noise and make registration challenging. Previous techniques mostly employed a pairwise registration approach, that is, all images of a sequence are individually registered to a fixed reference image [1–3]. The success of such approaches depends upon the robustness of the cost function to intensity change. Although intensity change due to contrast agent flow poses challenges, it also provides important temporal information for registration and segmentation of the perfusion image sequence. In this work, we propose a method which makes use of the temporal dynamics of contrast agent flow to achieve joint segmentation and groupwise registration of an MR cardiac image sequence.The changing intensity due to contrast agent flow provides valuable segmentation by highlighting the organ of interest. It is interesting to note that different regions of a cardiac image sequence have different intensity-time characteristics (Figure4). Temporal flow information was used for registration [4] and segmentation [5] of perfusion images. In [6], the registration framework constrained the deformation field such that different regions follow a particular intensity-time profile. But it is an accepted fact that improved registration leads to accurate segmentation and vice versa. In previous works [7, 8], we have approached the problem of joint registration and segmentation of cardiac MRI. However, we did not make use of the temporal information from the image sequence, rather focusing on pairwise registration of images. In this paper, we make use of temporal information to achieve registration and segmentation of the entire image sequence. Instead of pairwise registration, we solve the problem using groupwise registration as it allows us to impose constraints based on temporal information.Of late, groupwise registration methods have gained popularity because of a need to register large number of datasets for atlas construction [9, 10]. Groupwise registration is generally approached by two techniques. The first approach uses pairwise registration between a template and all other images in the population as in [11–13]. Yang et al. in [13] utilize the voxelwise geometry and orientation information for registration. Different measurements and rotation invariant features are used in a deformable matching mechanism for image alignment. However, pairwise registration has two limitations. First, selecting a fixed template may not accurately represent the population. Secondly, pairwise registration is not very effective when registering two images with significant anatomical differences. Only those subjects that are close to the template are registered properly.To overcome the above limitations, there are many methods that achieve registration using all images from the population. This approach is more faithful to the term “groupwise registration.” The goal is to warp all subjects in a population towards a hidden common space simultaneously within a single framework [14, 15]. The groupwise registration problem is formulated as one of optimization, with a global cost function defined on all aligned images [16, 17]. The cost function in [16, 17] is defined as the stack entropy or the entropy of corresponding voxels in different volumes. The justification is that if images are properly aligned, intensity values at corresponding voxels from all volumes will form a low entropy distribution. With this constraint, groupwise registration is achieved within a B-spline-based free-form deformation framework.When across-subject variation is very large, it is generally difficult to achieve good registration by simply registering each image to a template image. Therefore, many approaches make us of intermediate templates for registration [18–21]. In [18, 21], an intermediate template which does not belong to the original dataset is created to aid registration between two images. Tang et al. in [21] warp the template image with a set of simulated deformation fields learnt using principal component analysis (PCA) on a set of training deformation fields. Such approaches do not guarantee that the intermediate template is realistic which may affect registration results. Jia et al. in [19] construct a tree structure where each node of the tree is represented by an image, and similar images are represented by connected nodes. Each image is registered with the help of intermediate templates determined along its own path with respect to the final template.Groupwise registration methods are particularly suitable for perfusion DCE-MRI. Each region of the scanned organ is characterized by a different intensity profile over time. For example, in cardiac perfusion MRI, the contrast agent first flows into the right ventricle (RV) and then into the left ventricle (LV) before being flushed out of the cardiovascular system. Thus, a pixel within the RV shows a peak intensity magnitude early in the scanning sequence, while for the LV blood pool, the intensity peak occurs later (Figure4). With this available information, we can formulate the cost function such that, after registration, pixels from certain regions follow a particular intensity time profile. This is achieved by joint segmentation and groupwise registration of DCE-MRI. There are not many works dealing with joint segmentation and groupwise registration of cardiac perfusion MRI, although Zhang et al. in [22] describe a method for the rigid registration of brain perfusion images. The cost function is derived from the total quadratic variation of image intensity. Metz et al. in [23] propose a method for groupwise registration of dynamic lung data using both spatial and temporal constraints where groupwise optimization of B-splines is used.The temporal intensity patterns of pixels also determine their segmentation labels (i.e., RV, LV blood pool, myocardium background, etc.). Thus, the image sequence can be segmented along with groupwise registration. It is a well-known fact that registration and segmentation are mutually complementary approaches. Many works have combined them in a joint registration and segmentation framework [24–26]. Including segmentation information into the cost function reduces registration error and improves segmentation accuracy. In this paper, we propose a joint segmentation and groupwise registration (JSGR) approach for cardiac perfusion MRI. Our method combines intensity information from the entire image sequence (for registration) and maximizes the similarity between a pixel and other pixels belonging to the same class (for segmentation). We describe our method in Section 2, present experimental results in Section 3, and conclude with Section 4.


## 2. Theory

JSGR aims to find the transformation for every image to minimize a cost function. We do not have an explicitly defined reference image but constrain the transformations such that the registered images approach a common image space (which is the reference image) and are approximated to be the center of the images being registered. Wu et al. in [27] have highlighted the importance of a sharp mean image for accurate groupwise registration. Their observations are derived from constructing atlases for a large population of brain images. Such datasets show a lot of variability, and a fixed reference image is sure to introduce bias. Our method is used to achieve registration and segmentation using images from the same dataset, where each patient has multiple frames from different time intervals.The general objective function for JSGR consists of two terms, that is,(1)E=Edata+Esmooth,
where Edata is the data cost and Esmooth is the smoothness cost. The data cost depends upon the type of images being registered, and the smoothness cost depends upon the optimization framework.Perfusion images are characterized by rapid intensity change over time. Instead of relying on low-level information, we aim to exploit the temporal information for groupwise registration. First, we give a brief description of B-splines and their optimization. Then we explain the formulation of our data cost (Edata).


### 2.1. B-Spline-Based Registration

A B-spline-based freeform deformation (FFD) transformation model was presented in [28] for the elastic registration of breast images. The basic idea of FFDs is to deform an object by manipulating an underlying mesh of control points. The resulting deformation controls the shape of the 3D (or 2D) object and produces a smooth and continuous transformation. The transformation field consists of a global and local component and is defined as
(2)T(x)=Tlocal(Tglobal(x)),
where Tglobal is an affine transform obtained using [1], and Tlocal is the deformation based on B-splines. The DCE-MR images are first affinely registered to a chosen reference image. Note that this reference image is only for the purpose of rigid registration and is usually the scan showing all tissues without ambiguity. For elastic registration, there is no explicitly defined reference image. Further discussion is restricted to Tlocal.We define an initialnx×ny×nz grid of control points denoted as Φ. The grid points are denoted as Φi,j,k and have uniform spacing. The free-form deformation can be written as the 3D tensor product of 1D cubic B-splines,
(3)Tlocal(x)=x+∑l=03∑m=03∑n=03Bl(u)Bm(v)Bn(w)Φi+l,j+m,k+n,
where x=(x1,x2,x3) is the displacement vector, i=⌊x/nx⌋-1, j=⌊y/ny⌋-1, k=⌊z/nz⌋-1, u=x/nx-⌊x/nx⌋, v=y/ny-⌊y/ny⌋, w=z/nz-⌊z/nz⌋, and Bl is the lth cubic B-spline basis function given by the following equations:
(4)B0(u)=(1-u)36,B1(u)=(3u3-6u2+4)6,B2(u)=(-3u3+3u2+3u+1)6,B3(u)=(u3)6.Since our datasets are in 2D, the corresponding equations are(5)Tlocal(x)=x+∑l=02∑m=02Bl(u)Bm(v)Φi+l,j+m.B-splines are locally controlled, which makes them computationally efficient even for a large number of control points. In particular, the basis function of cubic B-splines have a limited support, that is, changing control pointΦi,j affects the transformation only in the local neighborhood of that control point.The local deformation of the cardiac tissues should be characterized by a smooth transformation. To constrain the spline-based FFD transformation to be smooth, one can introduce a penalty term which regularizes the transformation. The general form of such a transformation in 3D takes the following form:(6)Esmooth=1V∫0X∫0Y∫0Z[(∂2T∂x2)2+(∂2T∂y2)2+(∂2T∂z2)2+2(∂2T∂xy)2+2(∂2T∂xz)2+2(∂2T∂yz)2]dxdydz,
where V denotes the volume of the image domain. Our registration algorithm deals with 2D images in the spatial domain. Therefore, the smoothness cost can be formulated as
(7)Esmooth=1A∫0X∫0Y[(∂2T∂x2)2+(∂2T∂y2)2+2(∂2T∂xy)2]dxdy,
where A is the area of the image domain.


### 2.2. Similarity Measure

The images are first affinely aligned with respect to a reference image using [1]. Note that this reference image is only for rigid alignment and is not used for groupwise registration. Seed points belonging to RV, LV blood pool, myocardium, and background are identified (as shown by red arrows in the first image of the first row of Figure 3), and the labels of other pixels are determined using graph cuts. The initial labeling is used to calculate the cost functions in the first round of iteration. After the B-spline grid of each image is updated, the images are transformed and the segmentation labels are also updated based on the transformed images. We define the data cost as a combination of two terms which, individually, exploit the different characteristics of the perfusion datasets. Thus, Edata is defined as
(8)Edata=w1×EW+w2×EQ,
where EQ calculates the total quadratic variation of the dataset, EW calculates the within-class distance of each pixel, and w1,w2 are weights that determine the relative contribution of each term. w1=.4 and w2=1. Two weights are used in order to examine the relative contribution of each term to the results (discussed in Section 3.4). Below we explain each term in greater detail.


### 2.3. Total Quadratic Variation (EQ)

After registration of DCE-MRI, we expect the intensity time variation of pixels to be smooth after motion correction. Thetotal quadratic variation, EQ, measures the smoothness of the intensity signal by a combination of its first- and second-order derivatives. The sum of first derivatives over the entire sequence contributes to a smooth signal, while the sum of second derivatives favours a piecewise linear model of the signal. During the precontrast and postcontrast stages, or for regions without contrast enhancement, the first-order derivative is relevant because we expect the intensity of the same tissue to remain constant (i.e., the first-order derivative of a constant signal is zero). During the wash-in and wash-out stages, we expect the intensity of the same tissue to increase or decrease approximately with a constant rate. In an unregistered image sequence, the derivatives during the wash-in and wash-out stage will alternate between positive and negative values. However, as the image sequence is registered, the intensity changes are gradual and the sum of derivative values is minimal. We also use the second-order derivative of the intensity vector to encourage a piecewise linear intensity signal, as the second derivative of a linear signal is zero. Except for position of peak enhancement, the second derivative is zero at other time points.LetIn denote the image at the nth time point (or the nth frame in the dynamic sequence. The intensity of its ith pixel at time t is given by In,i(t), where t=1,…,T. Thus EQ is given by
(9)EQ=∑n=1N∑i∑t=0T(In,i′(t)+In,i′′(t))dt.
Here, I′ and I′′ are, respectively, the first and second derivatives of the intensity signal, and N is the total number of pixels in each image. Note that both the first and second derivatives are used for the entire image sequence.


#### 2.3.1. Within-Class Distance (EW)

The within class distance,EW, integrates segmentation information into the JSGR framework. EW ensures that pixels of the same class have similar intensity time profiles. In other words, pixels with the same label are made similar to a representative signal from that class. We use the mean intensity vector for representing each class of pixels. For a pixel i with known label l, its intensity vector provides greater information about the class labels. The within-class distance is calculated as the difference with respect to the mean intensity vector of class l(I¯l) and is given by
(10)EW=∑i‖Ii-I¯l‖.
Note that if pixel i belongs to class l (as determined from the current labels), then its difference only with respect to the mean intensity vector of class l is calculated; ∥·∥ denotes the Euclidean distance of the vectors. If a pixel has been correctly labeled as LV blood pool (or RV), then the residual error from the mean of LV blood pool (or RV) will be low. On the other hand, if the labeling is wrong, then the corresponding error is high. Initially, due to many unregistered pixels, the mean intensity vector may not be a very accurate representation of the class. But after every iteration, the mean intensity vector becomes smooth with the update of segmentation labels and starts to truly represent the particular class. In an iterative method, these constraints (EQ and EW) ensure that the labels converge correctly. Although the blood pool shows a lot of intensity change, it does not affect our method. By combining the contributions of the registration and segmentation terms, we overcome the effects of intensity change.Here, we need a representative intensity vector for a particular class. Although the mean vector is not necessarily the best representative vector at the beginning of registration, it is the best choice for a balance between registration accuracy and computational complexity. PCA could be a more accurate choice for the representative vector but significantly increases computation complexity. Moreover, we observe that with increasing number of iterations, the mean vector does converge to the representative vector of the class. This is supported by the fact that the registration and segmentation accuracies are the same when using the mean vector (after convergence) and principal components. But by using PCA, the computation time increases by more than 1.5 times.


### 2.4. Optimization

B-splines are used to optimize the cost function in (1). A uniform grid is initialized for all images. The first image’s grid coefficients are updated based on the present value of the energy function. The image is transformed, and the segmentation labels are immediately updated. The next image’s grid coefficients are then updated followed by image transformation and update of segmentation labels. This is repeated for all images of the sequence corresponding to different time points. This constitutes one iteration for JSGR. The segmentation labels are updated after the iteration to be used as the starting point for the next iteration. We repeat the process until the cost function does not decrease further.The advantage of dynamic update of segmentation labels is that it reflects the value of the energy function of the updated image sequence based on the latest transformations. Thus, if the energy function value does not change while updating the grid coefficients of three consecutive images, then JSGR can be immediately terminated. Once the process has been terminated, the final segmentation labels can be immediately obtained using graph cuts.In order that the images are registered to a common image space, the average deformation of a pixel is constrained to be zero. This is achieved by making the sum of coefficients for the corresponding grid point over all images to be zero. After the grid of every image is updated, the images are transformed, and the segmentation labels are dynamically updated, thus ensuring that the change in grid coefficients reflects immediately in the registration process. It also leads to more accurate registration than updating the segmentation labels at the end of each iteration. The algorithm is summarized in Algorithm1.Algorithm 1:Joint segmentation and groupwise registration framework.
(1) Determine the initial labelingLf of the image sequence.(2) Sett=1 for first registration round.(3) Initialize uniform B-spline grid for all images.(4) Calculate the value of energy function using the data cost formulation ((9) and (10)) and the smoothness cost (7).(5) Update the B-spline grid for the image, transform the image, and update the segmentation labels.(6) Repeat the above steps for all images of the sequence.(7) Repeat for all images until convergence.


## 2.1. B-Spline-Based Registration

A B-spline-based freeform deformation (FFD) transformation model was presented in [28] for the elastic registration of breast images. The basic idea of FFDs is to deform an object by manipulating an underlying mesh of control points. The resulting deformation controls the shape of the 3D (or 2D) object and produces a smooth and continuous transformation. The transformation field consists of a global and local component and is defined as
(2)T(x)=Tlocal(Tglobal(x)),
where Tglobal is an affine transform obtained using [1], and Tlocal is the deformation based on B-splines. The DCE-MR images are first affinely registered to a chosen reference image. Note that this reference image is only for the purpose of rigid registration and is usually the scan showing all tissues without ambiguity. For elastic registration, there is no explicitly defined reference image. Further discussion is restricted to Tlocal.We define an initialnx×ny×nz grid of control points denoted as Φ. The grid points are denoted as Φi,j,k and have uniform spacing. The free-form deformation can be written as the 3D tensor product of 1D cubic B-splines,
(3)Tlocal(x)=x+∑l=03∑m=03∑n=03Bl(u)Bm(v)Bn(w)Φi+l,j+m,k+n,
where x=(x1,x2,x3) is the displacement vector, i=⌊x/nx⌋-1, j=⌊y/ny⌋-1, k=⌊z/nz⌋-1, u=x/nx-⌊x/nx⌋, v=y/ny-⌊y/ny⌋, w=z/nz-⌊z/nz⌋, and Bl is the lth cubic B-spline basis function given by the following equations:
(4)B0(u)=(1-u)36,B1(u)=(3u3-6u2+4)6,B2(u)=(-3u3+3u2+3u+1)6,B3(u)=(u3)6.Since our datasets are in 2D, the corresponding equations are(5)Tlocal(x)=x+∑l=02∑m=02Bl(u)Bm(v)Φi+l,j+m.B-splines are locally controlled, which makes them computationally efficient even for a large number of control points. In particular, the basis function of cubic B-splines have a limited support, that is, changing control pointΦi,j affects the transformation only in the local neighborhood of that control point.The local deformation of the cardiac tissues should be characterized by a smooth transformation. To constrain the spline-based FFD transformation to be smooth, one can introduce a penalty term which regularizes the transformation. The general form of such a transformation in 3D takes the following form:(6)Esmooth=1V∫0X∫0Y∫0Z[(∂2T∂x2)2+(∂2T∂y2)2+(∂2T∂z2)2+2(∂2T∂xy)2+2(∂2T∂xz)2+2(∂2T∂yz)2]dxdydz,
where V denotes the volume of the image domain. Our registration algorithm deals with 2D images in the spatial domain. Therefore, the smoothness cost can be formulated as
(7)Esmooth=1A∫0X∫0Y[(∂2T∂x2)2+(∂2T∂y2)2+2(∂2T∂xy)2]dxdy,
where A is the area of the image domain.


## 2.2. Similarity Measure

The images are first affinely aligned with respect to a reference image using [1]. Note that this reference image is only for rigid alignment and is not used for groupwise registration. Seed points belonging to RV, LV blood pool, myocardium, and background are identified (as shown by red arrows in the first image of the first row of Figure 3), and the labels of other pixels are determined using graph cuts. The initial labeling is used to calculate the cost functions in the first round of iteration. After the B-spline grid of each image is updated, the images are transformed and the segmentation labels are also updated based on the transformed images. We define the data cost as a combination of two terms which, individually, exploit the different characteristics of the perfusion datasets. Thus, Edata is defined as
(8)Edata=w1×EW+w2×EQ,
where EQ calculates the total quadratic variation of the dataset, EW calculates the within-class distance of each pixel, and w1,w2 are weights that determine the relative contribution of each term. w1=.4 and w2=1. Two weights are used in order to examine the relative contribution of each term to the results (discussed in Section 3.4). Below we explain each term in greater detail.


## 2.3. Total Quadratic Variation (EQ)

After registration of DCE-MRI, we expect the intensity time variation of pixels to be smooth after motion correction. Thetotal quadratic variation, EQ, measures the smoothness of the intensity signal by a combination of its first- and second-order derivatives. The sum of first derivatives over the entire sequence contributes to a smooth signal, while the sum of second derivatives favours a piecewise linear model of the signal. During the precontrast and postcontrast stages, or for regions without contrast enhancement, the first-order derivative is relevant because we expect the intensity of the same tissue to remain constant (i.e., the first-order derivative of a constant signal is zero). During the wash-in and wash-out stages, we expect the intensity of the same tissue to increase or decrease approximately with a constant rate. In an unregistered image sequence, the derivatives during the wash-in and wash-out stage will alternate between positive and negative values. However, as the image sequence is registered, the intensity changes are gradual and the sum of derivative values is minimal. We also use the second-order derivative of the intensity vector to encourage a piecewise linear intensity signal, as the second derivative of a linear signal is zero. Except for position of peak enhancement, the second derivative is zero at other time points.LetIn denote the image at the nth time point (or the nth frame in the dynamic sequence. The intensity of its ith pixel at time t is given by In,i(t), where t=1,…,T. Thus EQ is given by
(9)EQ=∑n=1N∑i∑t=0T(In,i′(t)+In,i′′(t))dt.
Here, I′ and I′′ are, respectively, the first and second derivatives of the intensity signal, and N is the total number of pixels in each image. Note that both the first and second derivatives are used for the entire image sequence.


### 2.3.1. Within-Class Distance (EW)

The within class distance,EW, integrates segmentation information into the JSGR framework. EW ensures that pixels of the same class have similar intensity time profiles. In other words, pixels with the same label are made similar to a representative signal from that class. We use the mean intensity vector for representing each class of pixels. For a pixel i with known label l, its intensity vector provides greater information about the class labels. The within-class distance is calculated as the difference with respect to the mean intensity vector of class l(I¯l) and is given by
(10)EW=∑i‖Ii-I¯l‖.
Note that if pixel i belongs to class l (as determined from the current labels), then its difference only with respect to the mean intensity vector of class l is calculated; ∥·∥ denotes the Euclidean distance of the vectors. If a pixel has been correctly labeled as LV blood pool (or RV), then the residual error from the mean of LV blood pool (or RV) will be low. On the other hand, if the labeling is wrong, then the corresponding error is high. Initially, due to many unregistered pixels, the mean intensity vector may not be a very accurate representation of the class. But after every iteration, the mean intensity vector becomes smooth with the update of segmentation labels and starts to truly represent the particular class. In an iterative method, these constraints (EQ and EW) ensure that the labels converge correctly. Although the blood pool shows a lot of intensity change, it does not affect our method. By combining the contributions of the registration and segmentation terms, we overcome the effects of intensity change.Here, we need a representative intensity vector for a particular class. Although the mean vector is not necessarily the best representative vector at the beginning of registration, it is the best choice for a balance between registration accuracy and computational complexity. PCA could be a more accurate choice for the representative vector but significantly increases computation complexity. Moreover, we observe that with increasing number of iterations, the mean vector does converge to the representative vector of the class. This is supported by the fact that the registration and segmentation accuracies are the same when using the mean vector (after convergence) and principal components. But by using PCA, the computation time increases by more than 1.5 times.


## 2.3.1. Within-Class Distance (EW)

The within class distance,EW, integrates segmentation information into the JSGR framework. EW ensures that pixels of the same class have similar intensity time profiles. In other words, pixels with the same label are made similar to a representative signal from that class. We use the mean intensity vector for representing each class of pixels. For a pixel i with known label l, its intensity vector provides greater information about the class labels. The within-class distance is calculated as the difference with respect to the mean intensity vector of class l(I¯l) and is given by
(10)EW=∑i‖Ii-I¯l‖.
Note that if pixel i belongs to class l (as determined from the current labels), then its difference only with respect to the mean intensity vector of class l is calculated; ∥·∥ denotes the Euclidean distance of the vectors. If a pixel has been correctly labeled as LV blood pool (or RV), then the residual error from the mean of LV blood pool (or RV) will be low. On the other hand, if the labeling is wrong, then the corresponding error is high. Initially, due to many unregistered pixels, the mean intensity vector may not be a very accurate representation of the class. But after every iteration, the mean intensity vector becomes smooth with the update of segmentation labels and starts to truly represent the particular class. In an iterative method, these constraints (EQ and EW) ensure that the labels converge correctly. Although the blood pool shows a lot of intensity change, it does not affect our method. By combining the contributions of the registration and segmentation terms, we overcome the effects of intensity change.Here, we need a representative intensity vector for a particular class. Although the mean vector is not necessarily the best representative vector at the beginning of registration, it is the best choice for a balance between registration accuracy and computational complexity. PCA could be a more accurate choice for the representative vector but significantly increases computation complexity. Moreover, we observe that with increasing number of iterations, the mean vector does converge to the representative vector of the class. This is supported by the fact that the registration and segmentation accuracies are the same when using the mean vector (after convergence) and principal components. But by using PCA, the computation time increases by more than 1.5 times.


## 2.4. Optimization

B-splines are used to optimize the cost function in (1). A uniform grid is initialized for all images. The first image’s grid coefficients are updated based on the present value of the energy function. The image is transformed, and the segmentation labels are immediately updated. The next image’s grid coefficients are then updated followed by image transformation and update of segmentation labels. This is repeated for all images of the sequence corresponding to different time points. This constitutes one iteration for JSGR. The segmentation labels are updated after the iteration to be used as the starting point for the next iteration. We repeat the process until the cost function does not decrease further.The advantage of dynamic update of segmentation labels is that it reflects the value of the energy function of the updated image sequence based on the latest transformations. Thus, if the energy function value does not change while updating the grid coefficients of three consecutive images, then JSGR can be immediately terminated. Once the process has been terminated, the final segmentation labels can be immediately obtained using graph cuts.In order that the images are registered to a common image space, the average deformation of a pixel is constrained to be zero. This is achieved by making the sum of coefficients for the corresponding grid point over all images to be zero. After the grid of every image is updated, the images are transformed, and the segmentation labels are dynamically updated, thus ensuring that the change in grid coefficients reflects immediately in the registration process. It also leads to more accurate registration than updating the segmentation labels at the end of each iteration. The algorithm is summarized in Algorithm1.Algorithm 1:Joint segmentation and groupwise registration framework.
(1) Determine the initial labelingLf of the image sequence.(2) Sett=1 for first registration round.(3) Initialize uniform B-spline grid for all images.(4) Calculate the value of energy function using the data cost formulation ((9) and (10)) and the smoothness cost (7).(5) Update the B-spline grid for the image, transform the image, and update the segmentation labels.(6) Repeat the above steps for all images of the sequence.(7) Repeat for all images until convergence.


## 3. Experiments and Results

Cardiac images were acquired on a 1.5T Siemens Sonata MR scanner following bolus injection of Gd-DTPA contrast agent. There are 12 sequences acquired from 12 patients in whom it was important to look at myocardial perfusion. Each sequence comprised of 60 frames with a total of 720 frames. The datasets were acquired with electrocardiographic (ECG) gating such that the images were acquired during the same phase of the cardiac cycle. This minimized cardiac motion, but some deformations were still observed due to patient breathing. The pixel spacing ranges from(1.5×1.5) to (2.8×2.8) mm2. The acquired images were from the same midcavity slice. The images were corrected for rotation and translation motion before segmentation. The initial labeling was obtained using graph cuts [5]. The B-spline grid was of size 10×10 with the spacing between grid control points varying from 7 to 9 pixels.


### 3.1. B-Spline Grids

An example of B-spline-based illustration is shown in Figure1. A reference image, floating image, difference images before and after registration, and the deformed grid obtained after registration are shown. The initial grid is uniform and without any deformations. This example has been used for illustration purpose in which a reference image was fixed and all other images were registered using JRGS. From the difference images and deformed grid, we can conclude that the deformations of the floating image are captured by B-spline registration.B-spline registration results for cardiac perfusion images. (a) Reference image; (b) floating image; difference image; (c) before registration; (d) after registration usingJRGS; (e) deformed grids obtained from JRGS.
(a)(b)(c)(d)(e)


### 3.2. Registration Results

We compare the registration error fromJSGR with [16] (Met1), [28] (Met2), and our pairwise joint registration and segmentation method in [7] (JRS). Met1 does not define an explicit reference image and uses the entropy of the pixel stack for the cost function. Met2 uses pairwise registration with normalized mutual information (NMI) as the similarity measure. It has been used before for registering contrast-enhanced breast MRI. JRS is a joint registration and segmentation method, but it only uses information from a pair of images and not from the entire image sequence. Our algorithm converged after 6 iterations. The threshold cost difference above which registration continues is set at 0.1. We present qualitative and quantitative results of our method (JSGR) in terms of registration error and segmentation accuracy. For Met1, the number of iterations required for convergence was 5, while on an average, Met2 needed 4 iterations for convergence while registering a pair of images. JRS is a one-time optimization method and is not iterative. The reported results are after convergence of all algorithms.The outline of the LV (epicardium and endocardium) and RV is manually identified by expert observers in the original image sequence. These contours are denoted asCorg. The transformation of each image is used to map the contours to the registered image space. Let these contours be denoted as Ctrans. Since it is practically impossible to have ground truth value for elastic registration, we take a different approach to calculate registration error. In a perfusion image sequence, the shape of the heart should be same for ideal registration of all images. To calculate the registration error, we fix one image as reference (usually the first image of the sequence) and calculate the contour distance between the LV (or any object of interest) in all other images and the reference image. At the end of our JSGR method, not only has the sequence been registered but the segmentation labels of each frame have also been determined. Thus, we just need to find the contour distance between regions with the same labels. By contour, we refer to the outline of the segmented region. Note that the same procedure as above can be used to calculate the registration error before registration. The contour distance is calculated using the following steps:(1)
let the reference contour for the registered image sequence be denoted asCtransref,(2)
let theith point on Ctrans be denoted by Ctrans(i) and let the jth point on Ctransref be denoted by Ctransref(j),(3)
for everyCtrans(i) find the point on Ctransref(j) such that the distance between Ctrans(i) and Ctransref(j) is minimum(11)d(i,j)=min⁡j‖Ctrans⁡(i)-Ctrans⁡ref(j)‖,(4)
for each pointCtrans(i), the corresponding d(i,j) is calculated. The contour distance (CD) is the average distance and is defined as(12)CD=∑id(i,j)∑i.The average error measures for 12 datasets are given in Table1. Lower the registration error better is the method’s performance. We observe that JSGR has the lowest registration, while Met2 has the highest error. Met2 registers all images to a fixed reference image. Intensity change is common between two images of the dataset. Consequently, registration is prone to error while using a simple NMI-based similarity measure. Groupwise registration has the advantage that information from the entire sequence can be exploited for registration which is particularly important for DCE-MRI. The importance of dynamic information is highlighted by the fact that JRGS shows slightly higher accuracy than JRS which is a pairwise joint registration and segmentation method. JRS’s performance is better than Met1 although Met1 is a groupwise registration method. This indicates that inclusion of segmentation information improves registration more than the inclusion of dynamic information. Between Met1 and JSGR, the latter performs better because it combines segmentation and registration information in the cost function. Subsequently, the final registration errors for JSGR are lower than Met1. In Figure 2(a), we show the error measures for the LV for each of the 12 datasets.Table 1
Summary of registration and segmentation performance on cardiac perfusion datasets. The values indicate average and standard deviations for all datasets.
Registration error (mm)Segmentation resultBeforeAfter registrationDM (%)RegistrationJSGRMet1Met2JRSJSGRMet1Met2JRSLV epicardium2.2±1.20.7±0.11.0±0.21.1±0.10.8±0.493.6±0.790.1±0.390.6±0.692.1±1.1LV endocardium2.8±1.00.5±0.20.9±0.10.9±0.20.7±0.392.8±0.688.9±0.188.4±0.392.1±0.8LV overall2.6±1.10.5±0.20.9±0.31.0±0.10.7±0.493.1±0.589.4±0.288.8±0.892.1±0.9Registration and segmentation results of the LV for each of the 12 datasets. (a) Registration error (in pixels); (b) segmentation accuracy in%.
(a)(b)Contours of segmented regions overlaid on images from the dataset. First column shows results forJSGR, second column for Met1, third column for Met2, and fourth column for JRS. Each row corresponds to a different dataset.Intensity change with time for pixels on epicardium, RV, and endocardium. (a) Before groupwise registration; (b) after groupwise registration.
(a)(b)The time taken for registering one full dataset (including 60 images) is 1 hr 33 minutes usingJSGR, 58 mins using Met1, 45 mins using Met2, and 55 mins using JRS. JSGR and JRS were implemented using MATLAB 7.5 on a PC having a Pentium 4, 3 GhZ processor. Met1 was implemented by the authors using ITK and thus has low execution time. Met2 was also implemented in MATLAB for pairwise registration. Since JSGR is a joint segmentation and groupwise registration method, it takes more than twice the time compared to Met2.


### 3.3. Segmentation Results

Segmentation accuracy is calculated based on Dice metric (DM) values between manual segmentation and automatic segmentation for different methods. Manual segmentations for each slice are obtained after registration by each method. After registration is complete for each method, the segmentation labels are obtained by applying graph cuts on the intensity vectors. Now, the segmentation labels for corresponding pixels on different slices will be the same. These labels are compared with the manual segmentations to get the segmentation accuracy using DM. For JSGR, the labels are already obtained after the registration process.The averageDM values for the LV over all 12 datasets are shown in Table 1. Figure 2(b) shows the average DM values of the LV for each dataset over all frames of the sequence. The DM values are highest for JSGR thus indicating maximum accuracy amongst the three methods.Figure3 shows the segmented contours for LV endocardium overlaid on a representative image of the database. The representative image is chosen such that the blood pool is visible without any ambiguities. The manual segmentations are shown in red, while the automatic segmentations are shown in green. The first row also shows the initial segmentation in yellow. The initial segmentation is a result of applying graphcut to the unregistered image sequence and using the intensity vector of each pixel. Note that the automatic segmentations are the average contours over all frames of the sequence and the manual segmentations are also the average of the manually drawn contours. Since the segmentation labels are calculated from the intensity vectors, the labels will be the same for corresponding pixels in all frames. The first column shows results for JSGR, second column shows results for Met1, the third column shows results for Met2, and the fourth column shows results for JRS. Each row shows results for different datasets. Again we observe that JSGR shows the best agreement with manual segmentations due to the combination of registration and segmentation information. The other two methods being solely focused on registration perform inadequately for cardiac perfusion images. The accuracy measures highlight the importance of integrating registration and segmentation information. This combination is particularly important when images have low contrast, and also when segmentation information is available to be exploited for registration. Although the DM and registration error for JRGS and JRS are similar, a t-test gives P<0.032, thus indicating statistically different results, and hence improvement in results by using dynamic information.Figure4 shows the intensity variations with time for pixels on the LV blood pool, RV, and endocardium. Edge pixels are chosen for epicardium and RV, and therefore, the intensity change before registration (Figure 4(a)) is noisy. After registration using JSGR, the intensity variation is smoother (Figure 4(b)) and highlights the success of our method for time-varying data.


### 3.4. Importance ofEQand EW

It is important to look at the contribution ofEQ and EW to the overall registration procedure. EQ can be termed as the registration energy, while EW is the segmentation energy. We examine the improvement brought about by EW to registration accuracy, and also the improvement in segmentation accuracy due to EQ. We vary w1 (8) from 0 to 5 (keeping w2=1 fixed) and calculate the registration and segmentation accuracy values for the LV (shown in Table 2). It is observed that as w1 increases from 0, the registration and segmentation performance both improve. In fact, when w1≤0.4, the registration accuracy is less than or comparable to Met1 but improves with greater contribution of EW. However, if w1>2, the registration accuracy starts to degrade (evident from higher registration error), and the DM values also decrease due to unbalancing of each terms contribution.Table 2
Average registration and segmentation accuracy of LV with change inw1.
w100.10.20.411.522.5345Registration error (pix)10.90.70.60.60.70.70.80.91.21.5DM values (%)89.390.691.993.893.793.192.99189.989.388.8Similarly, when we increase the value ofw2 (with w1=0.4) from 0 to 5, we observe low DM values when w2<1. The best segmentation accuracy is obtained for 1≤w2≤3. However, when w2>3, the DM values start to decrease and the registration error also increases. Table 3 shows the change in registration error and DM values with change in w2. The best results are obtained for w2=1 and w1=0.4.Table 3
Average registration and segmentation accuracy of LV with change inw2.
w200.10.20.411.522.5345Registration error (pix)1.41.210.80.60.60.70.750.81.11.4DM values (%)88.289.590.491.793.893.592.792.191.889.288.3


## 3.1. B-Spline Grids

An example of B-spline-based illustration is shown in Figure1. A reference image, floating image, difference images before and after registration, and the deformed grid obtained after registration are shown. The initial grid is uniform and without any deformations. This example has been used for illustration purpose in which a reference image was fixed and all other images were registered using JRGS. From the difference images and deformed grid, we can conclude that the deformations of the floating image are captured by B-spline registration.B-spline registration results for cardiac perfusion images. (a) Reference image; (b) floating image; difference image; (c) before registration; (d) after registration usingJRGS; (e) deformed grids obtained from JRGS.
(a)(b)(c)(d)(e)


## 3.2. Registration Results

We compare the registration error fromJSGR with [16] (Met1), [28] (Met2), and our pairwise joint registration and segmentation method in [7] (JRS). Met1 does not define an explicit reference image and uses the entropy of the pixel stack for the cost function. Met2 uses pairwise registration with normalized mutual information (NMI) as the similarity measure. It has been used before for registering contrast-enhanced breast MRI. JRS is a joint registration and segmentation method, but it only uses information from a pair of images and not from the entire image sequence. Our algorithm converged after 6 iterations. The threshold cost difference above which registration continues is set at 0.1. We present qualitative and quantitative results of our method (JSGR) in terms of registration error and segmentation accuracy. For Met1, the number of iterations required for convergence was 5, while on an average, Met2 needed 4 iterations for convergence while registering a pair of images. JRS is a one-time optimization method and is not iterative. The reported results are after convergence of all algorithms.The outline of the LV (epicardium and endocardium) and RV is manually identified by expert observers in the original image sequence. These contours are denoted asCorg. The transformation of each image is used to map the contours to the registered image space. Let these contours be denoted as Ctrans. Since it is practically impossible to have ground truth value for elastic registration, we take a different approach to calculate registration error. In a perfusion image sequence, the shape of the heart should be same for ideal registration of all images. To calculate the registration error, we fix one image as reference (usually the first image of the sequence) and calculate the contour distance between the LV (or any object of interest) in all other images and the reference image. At the end of our JSGR method, not only has the sequence been registered but the segmentation labels of each frame have also been determined. Thus, we just need to find the contour distance between regions with the same labels. By contour, we refer to the outline of the segmented region. Note that the same procedure as above can be used to calculate the registration error before registration. The contour distance is calculated using the following steps:(1)
let the reference contour for the registered image sequence be denoted asCtransref,(2)
let theith point on Ctrans be denoted by Ctrans(i) and let the jth point on Ctransref be denoted by Ctransref(j),(3)
for everyCtrans(i) find the point on Ctransref(j) such that the distance between Ctrans(i) and Ctransref(j) is minimum(11)d(i,j)=min⁡j‖Ctrans⁡(i)-Ctrans⁡ref(j)‖,(4)
for each pointCtrans(i), the corresponding d(i,j) is calculated. The contour distance (CD) is the average distance and is defined as(12)CD=∑id(i,j)∑i.The average error measures for 12 datasets are given in Table1. Lower the registration error better is the method’s performance. We observe that JSGR has the lowest registration, while Met2 has the highest error. Met2 registers all images to a fixed reference image. Intensity change is common between two images of the dataset. Consequently, registration is prone to error while using a simple NMI-based similarity measure. Groupwise registration has the advantage that information from the entire sequence can be exploited for registration which is particularly important for DCE-MRI. The importance of dynamic information is highlighted by the fact that JRGS shows slightly higher accuracy than JRS which is a pairwise joint registration and segmentation method. JRS’s performance is better than Met1 although Met1 is a groupwise registration method. This indicates that inclusion of segmentation information improves registration more than the inclusion of dynamic information. Between Met1 and JSGR, the latter performs better because it combines segmentation and registration information in the cost function. Subsequently, the final registration errors for JSGR are lower than Met1. In Figure 2(a), we show the error measures for the LV for each of the 12 datasets.Table 1
Summary of registration and segmentation performance on cardiac perfusion datasets. The values indicate average and standard deviations for all datasets.
Registration error (mm)Segmentation resultBeforeAfter registrationDM (%)RegistrationJSGRMet1Met2JRSJSGRMet1Met2JRSLV epicardium2.2±1.20.7±0.11.0±0.21.1±0.10.8±0.493.6±0.790.1±0.390.6±0.692.1±1.1LV endocardium2.8±1.00.5±0.20.9±0.10.9±0.20.7±0.392.8±0.688.9±0.188.4±0.392.1±0.8LV overall2.6±1.10.5±0.20.9±0.31.0±0.10.7±0.493.1±0.589.4±0.288.8±0.892.1±0.9Registration and segmentation results of the LV for each of the 12 datasets. (a) Registration error (in pixels); (b) segmentation accuracy in%.
(a)(b)Contours of segmented regions overlaid on images from the dataset. First column shows results forJSGR, second column for Met1, third column for Met2, and fourth column for JRS. Each row corresponds to a different dataset.Intensity change with time for pixels on epicardium, RV, and endocardium. (a) Before groupwise registration; (b) after groupwise registration.
(a)(b)The time taken for registering one full dataset (including 60 images) is 1 hr 33 minutes usingJSGR, 58 mins using Met1, 45 mins using Met2, and 55 mins using JRS. JSGR and JRS were implemented using MATLAB 7.5 on a PC having a Pentium 4, 3 GhZ processor. Met1 was implemented by the authors using ITK and thus has low execution time. Met2 was also implemented in MATLAB for pairwise registration. Since JSGR is a joint segmentation and groupwise registration method, it takes more than twice the time compared to Met2.


## 3.3. Segmentation Results

Segmentation accuracy is calculated based on Dice metric (DM) values between manual segmentation and automatic segmentation for different methods. Manual segmentations for each slice are obtained after registration by each method. After registration is complete for each method, the segmentation labels are obtained by applying graph cuts on the intensity vectors. Now, the segmentation labels for corresponding pixels on different slices will be the same. These labels are compared with the manual segmentations to get the segmentation accuracy using DM. For JSGR, the labels are already obtained after the registration process.The averageDM values for the LV over all 12 datasets are shown in Table 1. Figure 2(b) shows the average DM values of the LV for each dataset over all frames of the sequence. The DM values are highest for JSGR thus indicating maximum accuracy amongst the three methods.Figure3 shows the segmented contours for LV endocardium overlaid on a representative image of the database. The representative image is chosen such that the blood pool is visible without any ambiguities. The manual segmentations are shown in red, while the automatic segmentations are shown in green. The first row also shows the initial segmentation in yellow. The initial segmentation is a result of applying graphcut to the unregistered image sequence and using the intensity vector of each pixel. Note that the automatic segmentations are the average contours over all frames of the sequence and the manual segmentations are also the average of the manually drawn contours. Since the segmentation labels are calculated from the intensity vectors, the labels will be the same for corresponding pixels in all frames. The first column shows results for JSGR, second column shows results for Met1, the third column shows results for Met2, and the fourth column shows results for JRS. Each row shows results for different datasets. Again we observe that JSGR shows the best agreement with manual segmentations due to the combination of registration and segmentation information. The other two methods being solely focused on registration perform inadequately for cardiac perfusion images. The accuracy measures highlight the importance of integrating registration and segmentation information. This combination is particularly important when images have low contrast, and also when segmentation information is available to be exploited for registration. Although the DM and registration error for JRGS and JRS are similar, a t-test gives P<0.032, thus indicating statistically different results, and hence improvement in results by using dynamic information.Figure4 shows the intensity variations with time for pixels on the LV blood pool, RV, and endocardium. Edge pixels are chosen for epicardium and RV, and therefore, the intensity change before registration (Figure 4(a)) is noisy. After registration using JSGR, the intensity variation is smoother (Figure 4(b)) and highlights the success of our method for time-varying data.


## 3.4. Importance ofEQand EW

It is important to look at the contribution ofEQ and EW to the overall registration procedure. EQ can be termed as the registration energy, while EW is the segmentation energy. We examine the improvement brought about by EW to registration accuracy, and also the improvement in segmentation accuracy due to EQ. We vary w1 (8) from 0 to 5 (keeping w2=1 fixed) and calculate the registration and segmentation accuracy values for the LV (shown in Table 2). It is observed that as w1 increases from 0, the registration and segmentation performance both improve. In fact, when w1≤0.4, the registration accuracy is less than or comparable to Met1 but improves with greater contribution of EW. However, if w1>2, the registration accuracy starts to degrade (evident from higher registration error), and the DM values also decrease due to unbalancing of each terms contribution.Table 2
Average registration and segmentation accuracy of LV with change inw1.
w100.10.20.411.522.5345Registration error (pix)10.90.70.60.60.70.70.80.91.21.5DM values (%)89.390.691.993.893.793.192.99189.989.388.8Similarly, when we increase the value ofw2 (with w1=0.4) from 0 to 5, we observe low DM values when w2<1. The best segmentation accuracy is obtained for 1≤w2≤3. However, when w2>3, the DM values start to decrease and the registration error also increases. Table 3 shows the change in registration error and DM values with change in w2. The best results are obtained for w2=1 and w1=0.4.Table 3
Average registration and segmentation accuracy of LV with change inw2.
w200.10.20.411.522.5345Registration error (pix)1.41.210.80.60.60.70.750.81.11.4DM values (%)88.289.590.491.793.893.592.792.191.889.288.3


## 4. Conclusion

We have proposed a novel method for joint segmentation and groupwise registration of the LV in cardiac perfusion images. By maximizing the smoothness of the temporal intensity signal, our method uses available temporal information from the entire image sequence. This helps to overcome the effects of intensity change. Segmentation information is incorporated by minimizing the error between a pixel’s intensity vector and mean intensity vector of the same class. Compared to manual segmentations, our method gives higher segmentation accuracy than other methods. The registration errors for our method are the lowest from amongst the three methods. Our being ajoint segmentation and groupwise registration approach, both registration and segmentation performance are better than conventional methods. This is because of two factors: (1) exploiting temporal information from DCE-MRI sequence in a groupwise registration framework and (2) use of mutually complementary registration and segmentation information, while most other methods solve registration and segmentation separately. Our method has the potential to be used for other data types having time-varying characteristics, and in future, we aim to use our method on other dynamic datasets.



---

*Source: 705853-2012-06-11.xml*