PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Nucl Sci. Author manuscript; available in PMC 2010 June 8.
Published in final edited form as:
IEEE Trans Nucl Sci. 2009 December 8; 56(6): 3662–3671.
doi:  10.1109/TNS.2009.2031642
PMCID: PMC2858412
NIHMSID: NIHMS184710

Estimation of Cardiac Respiratory-Motion by Semi-Automatic Segmentation and Registration of Non-Contrast-Enhanced 4D-CT Cardiac Datasets

Joyoni Dey, Member, IEEE,1 Tinsu Pan,2 David J. Choi,3 Dennis Robotis,4 Mark S. Smyczynski,1 P. Hendrik Pretorius, Associate Member, IEEE,1 and Michael A. King, Senior Member, IEEE1

Abstract

The goal of this work is to investigate, for a large set of patients, the motion of the heart with respiration during free-breathing supine medical imaging. For this purpose we analyzed the motion of the heart in 32 non-contrast enhanced respiratory-gated 4D-CT datasets acquired during quiet unconstrained breathing. The respiratory-gated CT images covered the cardiac region and were acquired at each of 10 stages of the respiratory cycle, with the first stage being end-inspiration. We devised a 3-D semi-automated segmentation algorithm that segments the heart in the 4D-CT datasets acquired without contrast enhancement for use in estimating respiratory motion of the heart. Our semi-automated segmentation results were compared against interactive hand segmentations of the coronal slices by a cardiologist and a radiologist. The pairwise difference in segmentation among the algorithm and the physicians was on the average 11% and 10% of the total average segmented volume across the patient, with a couple of patients as outliers above the 95% agreement limit. The mean difference among the two physicians was 8% with an outlier above the 95% agreement limit. The 3-D segmentation was an order of magnitude faster than the Physicians’ manual segmentation and represents significant reduction of Physicians’ time. The segmented first stages of respiration were used in 12 degree-of-freedom (DOF) affine registration to estimate the motion at each subsequent stage of respiration. The registration results from the 32 patients indicate that the translation in the superior-inferior direction was the largest component motion, with a maximum of 10.7 mm, mean of 6.4 mm, and standard deviation of 2.2 mm. Translation in the anterior-posterior direction was the next largest component of motion, with a maximum of 4.0 mm, mean of 1.7 mm, and standard deviation of 1.0 mm. Rotation about the right-left axis was on average the largest component of rotation observed, with a maximum of 4.6 degrees, mean of 1.6 degrees, and standard deviation of 2.1 degrees. The other rotation and shear parameters were all close to zero on average indicting the motion could be reasonably well approximated by rigid-body motion. However, the product of the three scale factors averaged about 0.97 indicating the possibility of a small decrease in heart volume with expiration. The motion results were similar whether we used the Physician’s segmentation or the 3-D algorithm.

Index Terms: Image segmentation, imaging applications, PET/CT, SPECT/CT

I. Introduction

The nature and magnitude of respiratory motion of the heart has been previously reported [1]–[11]. In [1] the predominant motion of the heart due to respiration was found to be along the superior-inferior (or foot-head) direction. Manke et al. [2] built a patient specific affine-model for cardiac respiratory motion for correction in coronary MR angiography. In [4] they applied multiple navigators to capture the variability across the breathing cycle. McLeish et al. [5] performed breath-hold cardiac MRI at different stages during the respiratory cycle for 8 volunteers and 9 patients to study motion and deformation of the heart due to respiration. In their methodology, one stage of respiration was hand segmented for each patient, and then registered with the remaining datasets. The largest extent of motion of the heart was found in the superior-inferior (SI) direction with a mean translation of 8.25 mm, observed in the 9 patients. This value was larger in the volunteers. Shechter et al. [6] analyzed cardiac motion in biplane coronary angiograms from 10 patients while normally breathing. Again the largest component of motion was SI translation, with a smaller mean value of 4.9 mm and maximum of 8.0 mm. Jahnke et al. [7] considered a total of 32 patients undergoing routine cardiac examination and applied different respiration motion correction techniques with increasing degree of freedom (just foot-head translation, all three translations and a fully affine-transformation). All motion compensation models showed significantly higher cross-correlations when compared to “no compensation”. In particular, the affine transformation algorithm achieved the highest cross-correlation values with a significant increase compared to a fixed 1-D translation (foot-head direction). The registration kernel was chosen to be a 3-D ellipsoid around the entire heart. A. P. King et al. [8] used a patient-specific model for correction of respiration in x-ray guided catheterization using a MRI roadmap. They build a affine-model of respiration from a quick pre-procedure scan registered to the high-resolution “road-map” scan. The structures of interest (such as 4 chambers and/or major vessels) were segmented in the high-resolution scan by using Analyze software. Boucher et al. [9] measured the amplitude of motion of the left ventricular apex as seen in respiratory-gated PET studies of 10 subjects. They observed a mean translation of 6.7 mm in the coronal plane. Livieratos et al. [10] measured the displacement of the centroid of the cardiac blood pool in 10 patients undergoing C15O PET imaging. Again the greatest motion was translational in the SI direction, with an average of 8.5 mm. Klein et al.. [11] investigated left ventricular motion in 9 subjects breathing normally during respiratory-gated PET imaging. They segmented the left-ventricle from one stage by thresholding, and then registered the remaining stages using a 12-degree-of-freedom (DOF) affine motion model. The calculated mean motion in the SI direction over 9 patients was 11.5 mm (leaving out one patient with exaggerated forced breathing).

We also have presented results employing preliminary algorithms for semi-automated segmentation and registration of the heart in non-contrast-enhanced CT studies, and investigated the effect of respiratory motion of the heart in several patients [12], [13]. Herein, we report an improved algorithm on motion estimation in a set of 32 patients imaged during normal tidal breathing by non-contrast-enhanced respiratory-gated 4D-CT. The 4D-CT datasets from each of these patients consist of 3-D data covering the heart region for 10 stages during the respiratory cycle. Our strategy for motion estimation was to first semi-automatically segment the heart region in the first of the 10 stages, and then estimation cardiac motion in the subsequent stages by registering the segmented stage to each of them. The output transformation parameters of the registration allowed us to estimate the gross 12-DOF affine motion of the heart at each stage due to respiration. It is emphasized that segmentation step was needed because the heart is adjoined by the lungs, liver, ribs, and other structures, each of which move independently from the heart and to different extents. We devised and report on a method of semi-automatic methods of cardiac segmentation for use in non-contrast enhanced CT studies. Our segmentation algorithm does not require the elaborate step of constructing an atlas from many training datasets, as used in [14], [15]. Instead, we use simple geometric shape-priors to initiate segmentation which is followed by a 3-D level-set based semi-automatic algorithm. Our segmentation methodology was compared against interactive hand segmentation performed by two physicians. The physicians’ segmentations were also used for independent respiratory motion estimation and the three motion estimates were compared. In the following we describe the acquisition of the 4DCT dataset employed, our segmentation methodology, the use of the segmentation cardiac regions to estimate motion by registration, validation of our segmentation methodology by comparison to the hand-segmentations of two physicians, and finally the estimated cardiac motion determined in this set of patients.

II. Respiratory Gated CT Acquisition

4D-CT datasets from 32 patients acquired using a GE Light-speed 8-slice CT system were obtained with IRB approval from the University of Texas, M.D. Anderson Cancer Center. The studies were not acquired for the purposes of this investigation. They were retrospectively selected as those which matched the needs of this investigation from those acquired clinically. The acquisition protocol has previously been described in detail [16]. Briefly, the x-ray collimation was 8 times 2.5 mm. CT scanning was conducted using cine-mode acquisition. Gantry rotation was 0.5 sec. The duration of data acquisition at each table position was the duration of the patient’s average respiratory cycle plus 1 sec. The reconstruction interval between two consecutive images was 0.45 sec or less so that there were over 15 samples per respiratory cycle. The RPM optical monitoring device from Varian Medical Systems, Palo Alto, CA was used to record the height of a reflector placed on the abdomen of the patients as a signal related to respiration. After image reconstruction, the cine CT images were correlated with the corresponding respiratory waveform to create 10 phases of 3-D CT images over a respiratory cycle, thus making the 4D-CT images by use of phase binning. All the patient data were stripped of patient identifiers in compliance with the Health Insurance Portability and Accountability Act (HIPAA). The volumes were 512 × 512 × S voxels, with the z-dimension (S) varying with different patients. The voxel dimensions were 0.98 mm × 0.98 mm × 2.5 mm. Before applying the segmentation algorithm, the CT-numbers of the datasets were windowed for each patient, to provide better visualization for interactive segmentation by physicians which was used to validate our semi-automatic segmentation. The process involved windowing such that CT-numbers in the range of −400 to 200 were mapped to 0–255. The windowing parameters were kept the same for all patients. We also cropped the datasets to reduce the data size, yet cover the heart adequately. The cropped dimensions were 201–220 voxels in transverse (xy) directions, and 32–81 voxels in axial (z) direction. One caveat is that the 4D-CT acquisitions were acquired using phase as opposed to amplitude binning, but the latter has been shown to result in fewer motion artifacts in CT slices [17]. The 32 patients included in this study were specifically selected based on their RPM tracings to have had a regular respiratory pattern during CT acquisition. For such regular-breathing patients the difference between phase and amplitude binning is not likely to be significant.

III. Segmentation

A. The Challenge

The CT datasets lack contrast between the blood-pool and wall. The heart has high contrast against the lungs; however, the densities of the heart and liver are similar. Thus, sometimes the slices are without clear heart-liver boundaries which make the segmentation of heart in non-contrast studies challenging. An example is given in Fig. 1(a). In this work a 3-D segmentation algorithm is presented. The approach was to assume an initial prior model for the heart shape as a first step in segmentation, with the assumption that the inherent spatial continuity of the model would help in segmenting the heart from the liver. To avoid building large training sets to find a prior model, our approach was to use a simple ellipsoid prior-shape that matched the major characteristics of the heart, as our starting point. This step was the only interactive part of the algorithm. Axial as well as coronal views were used to aid the placement of the prior. See Fig. 1(b) showing the cross section of the ellipsoid for the coronal slice and Fig. 1(g) showing that in a trans-axial slice for a second patient. Then a geometric active surface evolution was used to perform the segmentation [18]–[24]. This is done according to weighted minimal-surface distance map evolution as derived in [18]. The weights were inversely related to the gradient strengths, thereby attracting the surface to the nearest gradients. Fig. 1(c) and (h) shows the final segmentation for the examples. Considering Fig. 1(c), wherever the liver and heart have similar intensities and do not have too many strong gradients in between them, the surface would not evolve significantly there. Thus, the prior-shape boundary largely determined the indistinct parts of the boundary between the heart and the liver. The algorithm is described in more detail in what follows.

Fig. 1
Illustration of segmentation results for the first stage of respiration for three different segmentation methodologies for two different patients. The first row shows coronal slices for the first patient, and the second row shows transverse slices for ...

B. 3D Segmentation Algorithm

The 3-D Segmentation algorithm had 3 distinct steps.

1) First Step: 3-D Ellipsoid Prior Shape

First a 3-D ellipsoid is manually fit to the heart using coronal and axial views of the heart. At this very first step, we visually place the ellipsoid by manipulating the 12 degrees-of-freedom (DOF) of the ellipsoid to achieve the desired fit. The 12 DOF altered are three scale factors, three shears, three rotation angles, and finally three translations. In reality, just changing the three scales, the three translations and two angles (about y axis and z-axis) was sufficient to place the 3-D ellipsoid satisfactorily for all datasets considered. Our two main goals with this step were to clearly separate the heart and the liver inferior to the heart in 3D, and also to delineate the heart from other organs such as esophagus, ribs and so on, in the anterior-posterior (AP) direction. Some anatomical knowledge is necessary for placing the prior-shape about the heart. This is the only non-automated step in the 3-D algorithm. The next two steps are fully automatic.

2) Second Step: Level-Set Segmentation Bounded by 3-D Ellipsoid

Our initial surface is an ellipsoid. However at this stage, parts of this ellipsoid may be protruding into lungs; see Fig. 1(b). We would like to collapse these parts of the ellipsoid to the surface of the heart. We adopted a geometric active contour (surface) evolution method [18] for achieving this. This involves an implicit method of describing a surface as a zero level-set of a distance-map function, [var phi](x, y, z) of the surface from all the voxels. Creating a distance map for a general 3-D surface is computationally expensive. Ideally one would have to compute the nearest distance to the surface for all the voxels, which implemented in a brute force manner would take time proportional to Ns × Nx × Ny × Nz, where Ns is the number of points on the surface and (Nx, Ny, Nz) are the dimensions of the object. There are approximations that can be used to compute the distance map to a generalized surface [19], [24]. However in this work, we devised a simple, fast, intuitive method to collapse the ellipsoid to the lung-heart boundary as well as to compute the approximate distance map to this new surface in two fast steps, taking only roughly ~2 to 3 times (Nx × Ny × Nz).

We start by interpolating and “closing” the ellipsoid surface in all three directions. We next compute the mean m and standard-deviation σ inside a sphere of 5 pixel radius which is centered at the center of the ellipsoid and compute two thresholds, m ± 4σ. Then we start at the center of the spherfe and run a 3-D seed-growing algorithm to grow a region within these two thresholds, bounded by the closed surface of the ellipsoid. This region is called R henceforward. We used 6-connectivity in 3-D to obtain the region R. The resulting mask is a connected region of “0” and “1” inside the original ellipsoid, with “1” at voxel locations which were valued within the two thresholds. Now we proceed to find the approximate distance map for the outer-boundary of this mask.

In [18], weighted minimal surface evolution has been derived in 3-D in general to give the evolution equation

tφ(x,y,z)=g(I(x,y,z))div(φφ)φ+g(I)φ+vg(I)φ
(1)

where g(I(x, y, z)) is a general segmentation function, which will be (locally) minimized at the object boundaries, ν is the velocity term. The first two terms are in effect minimizing functional s g(I(xs, ys, zs))dS, where g(I(x, y, z)) is any appropriate function of the image I minimizing which would achieve the required segmentation. A popular [18], [24] “segmentation function” is g(I(x, y, z)) = 1/(1 + |[nabla]I(x, y, z)|p) where p is an appropriate power, usually p ≥ 1. It is also typical to use a smoothing operation on the image I before taking the gradient. The function then would be inversely related to the (gradient-based) edges of the image and the (1) then tends to evolve the surface such that this functional is minimized along voxels of the final surface. The velocity term helps the evolution when the surface is initialized (via an initial distance map) far from any edges. The final surface is the zero level-set of the (evolved) distance map function [var phi](x, y, z), that is the set of voxels for which [var phi](x, y, z) = 0. We used (1) with some minor modifications for our purpose (as explained in detail below).

First the mask m(x, y, z) defining the region R is smoothed with a Gaussian kernel to get another mask mg(x, y, z). The gradient of this smoothed mask [nabla]mg(x, y, z) will have high values at the boundary of the region R and zero or low values elsewhere. We compute a “segmentation function”

h(x,y,z)=11+mg(x,y,z)
(2)

This function h(x, y, z) has nearly a constant value of 1 where gradients are low or zero (non-boundary regions) and has low values at high gradients (boundaries). We then applied a modified version of the weighted minimal-distance surface evolution technique ((1)) to approximately find the distance map for the outer surface of the region (as explained below). We start with our small sphere (of radius 5 voxels) at the center of ellipsoid and analytically compute the distance map of the sphere to all the voxels. Then we evolve the sphere so that it minimizes h(x, y, z) along on its surface, or minimize the criterion s h(xs, ys, zs)dS [18], [24]. In other words, we want to find the surface such that the sum-total of the function h(xs, ys, zs) evaluated at all the surface points (xs, ys, zs) reaches a minimum.

For the purpose of finding the distance map of the boundary of the region R, instead of the outward velocity term being dependent on g(I) as in (1), we have made it proportional to m(x, y, z). That is, the velocity is “1” or “0” according to the unsmoothed mask of region R. We also found it useful to add an additional curvature flow term to make the surface evolution smoother. The iterative evolution used for Step 2 is therefore

φ(x,y,z)=φ(x,y,z)λf(h(I(x,y,z))div(φφ)φ+h(I)φ)λkdiv(φφ)φλvm(x,y,z)φ
(3)

The relative weights λf, λk and λv determine the relative importance of the “function” minimizing term, the curvature term, and the velocity terms, respectively. We used λf = 0.05, λk = 0.05 and λv = 0.1, found by trial and error. We found that 500 iterations were adequate to obtain an approximation of the distance map of the boundary of R for all the patients.

3) Third Step: Fine-Tuning Surface to Nearest Edge

We obtained a fairly good segmentation of the heart by just the region R. The protrusions into low-CT-number regions of the lungs were collapsed. However we may still have cases where the prior-shape boundary fell inside the heart. As a last step here we fine tune the surface obtained in Step 2 to the nearest edges in the original dataset. Here we use the weighted minimal-surface evolution given in (3), by minimizing the function

h(x,y,z)=11+I(x,y,z)
(4)

Note that in Step 3, compared to Step 2, we use a different, h(I(x, y, z)), composed of gradients of the original intensities and not those of the mask. We also do not use the ballooning force (the velocity term) since that would defeat the purpose of the delineation of the organs by the prior-shape and make the segmentation leak into the liver. Thus, in our evolution equation for Step 3, (3), was used with the values λf = 0.001, λk = 0 and λv = 0 and the function given in (4). This is also equivalent to the original (1) with the velocity term ν = 0 and g given by the function in (4). The number of iterations used was 100 for all patients.

At the regions where the intensities were similar across the boundaries of organs, such as the heart-liver interface, there would not be strong gradients. Hence, the segmentation at those regions remains largely determined by the prior-shape chosen in Step 1. A further step could possibly be to introduce another “shape-prior” term in (1) similar to the idea presented in [22], such as, β([var phi]; − [var phi]*), where [var phi]* is the distance map of the prior-shape, obtained by the Step 2. This term would make the evolving distance map stay close to the prior-shape with the parameter β determining the degree of importance to give the prior relative to the other terms. However, for our purposes we found it was not necessary to use this prior-shape term.

IV. Registration and Motion Estimation

Once the first stage of respiration was segmented for each patient, the 12-DOF affine transformation registering the other stages to the first was determined as the mechanism for estimating respiratory motion of the heart in that patient. The usage of the 12-DOF affine transformation model to represent respiratory motion of the heart is based on the studies of cardiac respiratory motion in breath-hold MRI by McLeish et al. [5], and the determination by Manke et al. [2]–[4] of its utility in correcting for respiratory motion in MR coronary angiography and also by A.P. King et al. [8] for X-ray and MRI-roadmap guided cardiac catheterization. Its usage is to reduce the motion to 12 parameters for reporting and modeling purposes. The 12 DOF affine transformation matrix T is given by

T(x,y,z)=(s1sh1sh2s2sh3s3)(R)(xyz)+(txtytz)
(5)

The first term on the right is an upper-triangular matrix with the three scale (s1, s2, s3) and the three shear (sh1, sh2, sh3) parameters. The second term is the rotation matrix R, composed of three rotation matrices about the x, y, and z axes. The final term is the translations. The rotations and translations were obtained about and along the x, y, and z axes which correspond to the Right-Left (RL), anterior-Posterior (AP), and Superior-Inferior (SI) axes, respectively. The center of rotation was taken as the center of segmented heart volume of the first stage of respiration for the patient. The initial three scaling values and three shear values were set at 1.0 and 0.0 respectively.

We use intensity based image registration [25]. The criterion used to determine the transformation matrices was the minimization of the sum squared difference (SSD) of the intensities which was calculated as

SSD=n[I1(x,y,z)I2(T(x,y,z))]2
(6)

where I1(x, y, z) are the intensities in the dataset for the first stage of respiration (end inspiration) and I2(x, y, z) are the intensities for the other datasets (the Stages 2–10 in turn) which are registered to the first. A gradient descent optimizer was used iteratively. Analytically derived gradients were used for minimization and parameters were updated with iteration i according to Ti+1 = TiL(d(SSD)/dTi) Here, the L are the learning rates, the order of magnitude of which depend on the image intensities. The relative sensitivity of the different parameters and algorithm speed also determines these learning rates. For all transformation estimations, only the voxels that correspond to non-zero voxels in the first (segmented) stage were considered. Tri-linear interpolation was used to interpolate all the voxels of the second dataset into the transformed co-ordinates. The applied transformation was actually the inverse of the estimated transformation. In other words, as our algorithm progressed, the second dataset was transformed and gradually registered to Stage-1, while the estimated transformation (inverse of that applied to second dataset, for every iteration) was actually the transformation necessary to register Stage-1 to the second dataset. The registration algorithm was executed in batches on a Linux cluster of more than 170 nodes. Each node had dual AMD processors and at least 2G of RAM.

Before registration all the volumes were re-sampled to obtain a uniform voxel grid of dimension 0.98 mm in all directions. The 9 other respiratory stages were registered with the first stage to first obtain the 12-DOF motion parameters. For the purpose of validation of the registration (as described in next section) the volumes were re-sampled back to original dimensions (0.98 mm × 0.98 mm × 2.5 mm).

V. Validation

A. Validation of Segmentation Algorithm

For each patient, a comparison was made between the semiautomatic 3-D segmentation and that of two physicians (a cardiologist and a radiologist) involved in this investigation. The physicians used a graphical-user-interface (GUI) to manually trace out the heart region on the CT coronal slices of the first stage of respiration using a mouse. They segmented all the slices they judged to be clearly belonging to the heart. The physicians included some of the blood-vessels such as the aorta and the inferior vena cava during interactive segmentation, since it was difficult to consistently exclude them.

The evaluation method was based on the method by Bland and Altman [26], where to compare across two methods the differences are plotted versus the means, for all the samples. It is typical also to draw the average of the difference across the samples, and ±1.96σ lines, which would indicate a 95% acceptance limit if the distribution of the difference was normal. We normalize the difference and mean by the average segmented volume of all the patients by all the three methods.

When evaluating the algorithm, any slices that the physicians left out un-intentionally or were unable to segment were automatically set to zero and not calculated as errors. The top few rows from each coronal slice were systematically cropped for all datasets, so as to be consistent in the amount of aorta included by the three different segmentations.

B. Validation of Overall Method: Registration and Segmentation

One way to validate our registration method would be to determine if selected features aligned at the same location in the reference dataset (Stage-1 or end-inspiration) and at the other stages of respiration once they had been registered to Stage-1. However, our clinical co-authors advised that it would be difficult if not impossible to reliably determine the location of features between different stages of respiration in these non-contrast-enhanced datasets. Thus, instead the entire contour boundaries of the heart for (the already segmented) Stage-1 and that segmented in Stage-7 (approximately end-expiration) after registration to Stage-1 were compared. The contour boundaries of the registered version of Stage-7 were obtained using our semi-automated segmentation algorithm, after the registered slices of Stage-7 were re-sampled back to their original dimensions (0.98 × 0.98 × 2.5 mm). Then the segmentation errors were computed with respect to the stage 1 segmented dataset by Bland and Altman method [26] s discussed in previous section.

VI. Results

A. Segmentation and Registration

The results of our semi-automated segmentation method were compared with those of the two physicians denoted as Physician 1, and Physician 2. By qualitative visual inspection of the coronal slices, the similarity of the segmentation appeared good overall between the three methods. The separation of the liver and heart was very good. The inclusion of major blood-vessels (mainly the aorta and inferior vena cava) varied between and across the two Physicians and the algorithm. Considering trans-axial slices, the 3-D algorithm (Algo) was observed to be similar to the segmentation of the two Physicians. In the anterior direction, Algo and Physician 1 had comparable results. Physician 2 did err on the side of caution for some patients and excluded some voxels in the anterior direction.

Results for two patients representative of the different sizes and different noise-levels in our datasets are shown in Fig. 1. Shown in each Fig. 1(c)–(e) is a coronal Stage-1 CT slice from a patient superposed by the contour estimated by 3-D segmentation algorithm, and the contours drawn by Physician 1 and Physician 2, respectively. Fig. 1(h)–(j) shows axial slices for a different patient. These results show the Algo as the smoothest of the three in the trans-axial direction and that all three algorithms include most of the same pixels.

Fig. 2 shows the plots of the normalized “difference” versus the normalized “mean” as described in Section IV-A, for quantitative comparison of the degree of agreement between the pair-wise comparison of the Algorithm and the two Physicians. Since the ambiguity in the blood-vessel segmentation superior to the heart was inconsistent between the physicians themselves and between the algorithm and the physicians and not very relevant to the segmentation of the heart, some of the voxels in this region were cropped for all three segmented datasets before comparing them quantitatively. Fig. 2 top-row shows the normalized difference versus normalized mean plotted for the pair-wise comparisons between the Algo, Physician 1, and Physician 2. We observed that most of the patients have differences below the mean. There are a couple of outliers above the upper-limit of the 95% agreement line, which is driving up the value of the mean and the standard-deviation. The bottom-left shows the Physician 1 and Physician 2 comparison. This has a tighter 95% bound and one outlier. The mean and the standard-deviation (across patients) of the normalized difference between Algo and Physician 1 is given by (0.11, 0.056) and that between Algo and Physician 2 is given by (0.095, 0.049). The mean and variance between Physician 1 and Physician 2 is given by (0.08, 0.034). Note that the means indicate on the average the percentage differences are about 11%, 10% and 8% of the average volume across all patients for the pair-wise comparison of Algo & Phys1, Algo & Phys2 and Phys 1 & Phys 2 respectively. However the variation across patients about the mean is smaller for the pair-wise comparison of Phys 1 & Phys 2.

Fig. 2
Two-way comparisons of segmentation normalized difference versus normalized mean between (a) Physician 1 (Phys1) and 3-D Algorithm (Algo); (b) Physician 2 (Phys2) and Algo; (c) Phys1 and Phys2. In (d) we show the segmentation results for comparison between ...

To obtain a quantitative measure of registration accuracy across all patients, Stage 7 after registration to Stage 1 was segmented for all the patients using the 3-D segmentation algorithm. A comparison was then made between voxels that were segmented in the registered Stage 7 to those segmented in Stage 1. For this case we normalized by the voxels in the original Stage 1. Fig. 2(d) shows the plots of the normalized difference versus the normalized mean of the two segmentations by the Algo. The mean and the standard variation for the normalized difference is (0.07, 0.036). The match is therefore within 7% of the average volume across all patients, with a couple of outliers.

B. Cardiac Respiratory Motion Estimation Results

For respiratory motion estimation via registration a fixed number of iterations of 3000 was used. For the three scaling and rotation parameters, the learning rates used for most patients were 5 × 10−6 while those for translations were 5 × 10−4. For the shear parameters they were 2 × 10−7.

The parameters of the 12-DOF registration constitute the respiratory motion estimates for the hearts in the 4D-CT datasets. For all patients, the magnitude of motion relative to end-inspiration (Stage 1) steadily increased till it maximized around Stage 6, 7, or 8, which was thus end-expiration for that patient. It then decreased back towards end-inspiration. The majority peaked around Stage 7. The largest displacement was superiorly along the SI axis. These trends are seen in Figs. 35 which show the mean 12 DOF parameters as a function of respiratory stage for the 3 different methods of segmentation. For display purposes, the 12-parameters were grouped into plots of three scales, three shears, three rotation angles and three translations. Stage 1 was the end-inspiration stage and all motions were relative to it. The first 6 to 7 stages belong to the expiration cycle, which was typically longer than inspiration. Typically Stage 6, and 7 (and on rare cases, Stage 8) was end-expiration. Since the SI translation was the most significant parameter, its one-standard-deviation bar was shown in these figures. Notice that the variation of the mean values for the other 11 parameters is small in comparison.

Fig. 3
The mean 12-DOF respiratory motion parameters resulting from use of 3-D Algorithm for segmentation, over the 32 patients plotted as a function of respiration stage. Note the motion in each case is relative to Stage 1 which is end inspiration. The 12 parameters ...
Fig. 5
The mean 12-DOF respiratory motion parameters resulting from use of Physician 2’s segmentation, over the 32 patients plotted as a function of respiration stage. Note the motion in each case is relative to Stage 1 which is end inspiration. The ...

Table I shows the mean, SD, and maximum values for the 12 parameters for Stage 7 for the three registration methodologies (using the 3-D algorithm segmentation, and the two Physicians). The SI translation was the largest with a mean value of 6–7 mm (depending on segmentation methodology), and maximum of 10–11 mm. The AP translation was the next largest with a mean of 1.5–2 mm, and a maximum of up to 5 mm. The lateral translation was the least with an average value of less than 1 mm, but maximums of over 5 mm. Of the rotations, the roll about the RL axis was the most significant, followed by a flip around the AP-axis; however, these are very small. The product of the scale-terms (D) was also calculated to determine if a volume change was detected. This parameter is shown in the last column of Table I. Overall D was found to be close to unity (0.96–0.97).

TABLE I
Mean, SD, Max of Parameters of Respiratory Motion Estimated Between Stage 1(End Inspiration) and Stage 7(Approximately End Expiration) for 3-D Algorithm (Algo) and That of the Two Physicians (Phys1 and Phys2). Scl is Scaling in the x (Right-Left), y (Posterior-Anterior), ...

VII. Discussion

While it is inexact to compare across studies because of differences in ways of controlling respiratory for imaging and motion estimation methodologies, the mean and maximum values (6.4 mm and 10.7 mm) of the SI motion parameters obtained herein were similar to what has been reported in the literature. Our estimates lie in between the free tidal breathing estimates from 8 coronary angiograms for patients by Shechter et al. [6] (mean 4.9 mm, max 8 mm) and the values for breath-holding at extents of normal respiration for 9 patients reported by McLeish et al. [5] (mean 8.3 mm, max 14.3 mm).

For respiratory motion estimation the segmentation was necessary to eliminate the liver and other structures which move differently with respiration and may lead to errors in motion estimation. To achieve this we developed a 3-D semi-automatic cardiac segmentation algorithm for non-contrast-enhanced CT studies, where we place an initial ellipsoid and then run a minimal-surface-based segmentation algorithm.

In the higher gradient areas such as against the lungs the ellipsoid placement is not very critical since the subsequent steps would achieve a good job of accurate segmentation. To aid the placement of the ellipsoid to delineate the lower-gradient regions such as the liver-heart boundary, we had the option of full 12-DOF affine parameters (scales, shears, rotations, and translations) to place the ellipsoid. However we found that we needed only a subset of those parameters to place the ellipsoid satisfactorily to delineate the heart and the liver for all the patients. Note the resultant ellipsoid may protrude into the lung-area, which is handled in the subsequent steps of segmentation. The task for which the semi-automatic segmentation was developed was to prune-down the voxels being compared in the registration step to those of the heart for purpose of motion estimation. While a perfect segmentation is not critical for this task, our success in achieving segmentations is seen in the good agreement with those of 2 experts.

The motion estimates were similar using the 3-D algorithm for segmentation vs the segmentations of the Physicians. Thus, it appears our segmentation methodology was sufficient for its purpose. The advantage of using the 3-D algorithm is that it is an order of magnitude faster than the manual segmentation by the physicians. The only non-automated part of the algorithm is placing the prior-ellipsoid and it is possible to achieve than within minutes. The Physicians took about an hour for each patient.

We estimated the rotation, scaling and shear about the centroid of the first segmented volume. Therefore, the translation estimates are effectively only for the centroid. Other points could be moving more. The apex of the pericardial sack of the heart for example would be flapping and rolling due to the rotation and having more movement for some patients. Note that in this study the patients were CT imaged at rest. Larger respiration motion would logically be expected following physical exercise as per patients undergoing exercise-stress imaging. Similarly the period of time taken for CT imaging is short compared to that of SPECT imaging. Thus, in SPECT patient have more time to alter their respiratory pattern during imaging and this will increase the impact of respiration in studies when it occurs.

One potential limitation of this study is that cardiac motion was ignored in this work. The cardiac cycle is typically less than 1 sec and the respiratory cycle is approximately 4–6 sec. We employed a time between successive CT image reconstructions of 0.45 sec or less. Thus, unless the cardiac cycle is an exact multiple of the time between successive reconstructions, the binning process would to some extent average over the cardiac cycle. Certainly the impact of sampling would be greater on cardiac motion than on respiratory motion. Also there is no blood-pool contrast in the datasets. Thus, we are essentially seeing the respiratory motion of the pericardial sack of the heart, and not the contractions and expansions of the myocardium. It has been demonstrated that the contents of the pericardial sack are relatively constant during contraction in the intact human thorax [27]. This is the result of the atria filling as the ventricles move towards the apex causing the myocardial walls to thicken, with the combined effect pumping out blood. Also it has been determined that that cardiac gating of transmission imaging for attenuation map formation is not needed when using both respiratory and cardiac gating of PET emission studies [28]. Thus, it seems likely that our lack of having available combined cardiac and respiratory datasets should not be a major limitation of the motion estimated in this study, but may have had some impact.

Given the magnitude and variability of cardiac motion due to respiration motion determined herein one can expect the impact will vary with patients and imaging modality. For cardiac SPECT, a mild impact on image quality is expected for patients imaged at rest, since the 1 cm or greater system spatial resolution will dominate image blurring. A greater impact would be expected for cardiac PET due to its superior spatial resolution compared to SPECT. One caveat is that cardiac emission acquisitions are much longer than the CT acquisitions used here. Hence, there could be potentially more variation in respiration during emission imaging resulting in a greater impact than might be expected based on the results of this study. For example, “upward creep” of the heart during SPECT due to a change in patient respiration primarily following exercise-stress has been know as a cause of artifacts [29].

VIII. Conclusion

In this work, motion of the heart due to respiration in 32 patients was investigated. To perform this investigation we developed a semi-automated 3-D segmentation algorithm which overcomes some of the challenges of segmenting non-contrast enhanced images of the heart. Our semi-automated segmentation results were compared against interactive hand segmentations of coronal slices by a cardiologist and a radiologist. On pair-wise comparison of the difference of segmentation between the Algo & Physician1, Algo & Physician2 and Physician1 & Physician2 are 11%, 10% and 8% respectively of the average segmented volume across the patients. The 3-D segmentation was an order of magnitude faster than the Physicians’ manual segmentation and represents significant reduction of Physicians’ time. The segmented first stages of respiration were registered to the other respiratory stages for 32 patients to obtain 12-DOF respiratory motion parameters. Our results indicated superior-inferior (SI) axis translation to be the largest component of motion with a maximum of 10.7 mm (moving superiorly on expiration) and a mean SI translation of 6.4 mm. Translation in the anterior-posterior direction was the next largest in magnitude with an average of 1.7 mm and maximum of 4.0 mm (moving posteriorly with expiration). Translation along the right-left axis was the least on average. A rotation of 1.6 degrees on average about the right-left axis was noted with a maximum rotation of 4 degrees. The other rotations, scaling, and shear parameters were all close to zero on average indicting the motion could be reasonably well approximated by rigid-body motion. However, the product of the three scale factors averaged about 0.97 indicating the possibility of a small decrease in heart volume with expiration. The motion parameters were similar whether we used the Physician’s segmentation or the 3-D algorithm.

Fig. 4
The mean 12-DOF respiratory motion parameters resulting from use of Physician 1’s segmentation, over the 32 patients plotted as a function of respiration stage. Note the motion in each case is relative to Stage 1 which is end inspiration. The ...

Acknowledgments

This work was supported by the National Heart, Lung, and Blood Institute under Grant HL50349, and by the National Institute of Biomedical Imaging and Bioengineering under Grant EB001457. The contents are solely the responsibility of the authors and do not necessarily represent the official views of the National Heart, Lung, and Blood Institute or the National Institute of Biomedical Imaging and Bioengineering.

References

1. Wang Y, Riederer SJ, Ehman RL. Respiratory motion of the heart: Kinematics and the implications for the spatial resolution in coronary imaging. Magn Reson Med. 1995;33(3):713–719. [PubMed]
2. Manke D, Rösch P, Nehrke K, Börnert P, Dössel O. Model evaluation and calibration for prospective respiratory motion correction in coronary MR angiography based on 3-D image registration. IEEE Trans Med Imag. 2002 Sep;21(9):1132–1141. [PubMed]
3. Manke D, Nehrke K, Börnert P, Rösch P, Dössel O. Respiratory motion in coronary magnetic resonance angiography: A comparison of different motion models. J Mag Res Imag. 2002;15(6):661–671. [PubMed]
4. Manke D, Nehrke K, Börnert P. Novel prospective respiratory motion correction approach for free-breathing coronary MR angiography using a patient-specific adaptive affine motion model. Mag Res Med. 2003;50:122–131. [PubMed]
5. McLeish K, Hill DLG, Atkinson D, Blackall JM, Razavi R. A Study of the motion and deformation of the heart due to respiration. IEEE Trans Med Imag. 2002;21(9):1142–1150. [PubMed]
6. Shechter G, Ozturk C, Resar JR, McVeigh ER. Respiratory motion of the heart from free breathing coronary angiograms. IEEE Trans Med Imag. 2004;23:1046–1056. [PMC free article] [PubMed]
7. Jahnke C, Nehrke K, Paetsch I, Schnackenburg B, Gebker R, Fleck E, Nagel E. Improved bulk myocardial motion suppression for navigator-gated coronary magnetic resonance imaging. J Mag Res Imag. 2007;26:780–786. [PubMed]
8. King AP, Boubertakh R, Rhode KS, Ma YL, Chinchapatnam P, Gao G, Tangcharoen T, Ginks M, Cooklin M, Gill JS, Hawkes DJ, Razavi RS, Schaeffter T. A subject-specific technique for respiratory motion correction in image-guided cardiac catheterisation procedures. Med Imag Anal. 2009;13:419–431. [PubMed]
9. Boucher L, Rodrigue S, Lecomte R, Benard F. Respiratory gating for 3-dimensional PET of the thorax: Feasibility and initial results. J Nucl Med. 2004;45:214–219. [PubMed]
10. Livieratos L, Rajappan K, Stegger L, Schafers KK, Bailey DL, Camici PG. Respiratory gating of cardiac PET data in list-mode acquisition. Eur J Nucl Med Mol Imag. 2006;33(5):584–588. [PubMed]
11. Klein GJ, Reutter BW, Huesman RH. Four-dimensional affine registration models for respiratory-gated PET. IEEE Trans Nucl Sci. 2001;48(3):756–760.
12. Dey J, Pan T, Smyczynski M, Pretorius PH, Choi D, King MA. Investigation of respiration motion of the heart based on semi-automated segmentation and modeling of cardiac-CT data. IEEE Nucl. Sci. Symp. Med. Imag. Conf. Rec; Puerto Rico. Oct. 2005.
13. Dey J, Pan T, Choi DJ, Smyczynski M, Pretorius PH, King MA. Knowledge-based segmentation of the heart from respiratory-gated CT datasets acquired without cardiac contrast-enhancement. Proc. 2006 SPIE 6144, Med. Imag; San Diego, CA. Feb. 11–16, 2006.
14. Park H, Bland PH, Meyer CR. Construction of an abdominal probabilistic atlas and its application in segmentation. IEEE Trans Med Imag. 2003 Apr;22(4):483–492. [PubMed]
15. Shimizu A, Ohno R, Ikegami T, Kobatake H, Nawano S, Smutek D. Segmentation of multiple organs in non-contrast 3-D abdominal CT images. Int J Comput Assist Radiol Surg. 2007 Dec;2:135–142.
16. Pan TS, Lee TY, Rietzel E, Chen TY. 4D-CT imaging of a volume influenced by respiratory motion on multi-slice CT. Med Phys. 2004;31:333–340. [PubMed]
17. Lu W, Parikh PJ, Hubenschmidt JP, Bradley JD, Low DA. A comparison between amplitude sorting and phase-angle sorting using external respiratory measurement for 4D CT. Med Phys. 2006(33):2964–2974. [PubMed]
18. Caselles V, Kimmel R, Sapiro G, Sbert C. Minimal surfaces based object segmentation. IEEE Trans Pattern Anal Mach Intell. 1997;19(4):394–398.
19. Sethian JA. Level Set Methods and Fast Marching Methods. Cambridge, U.K: Cambridge Univ. Press; 1999.
20. Malladi R, Sethian JA, Vemuri BC. Shape modeling with front propagation: A level set approach. IEEE Trans Pattern Anal Mach Intell. 1995;17(2):158–175.
21. Xu C, Yezzi A, Prince JL. A summary of geometric level-set analogues for a general class of parameteric active contour and surface models. Proc. 1st IEEE Workshop on Variational and Level Set Methods in Comput. Vis; Jul. 2001; pp. 104–111.
22. Leventon M, Grimson E, Faugeras O. Statistical shape influence in geodesic active contours. Proc. IEEE Conf. Comput. Vis. Pattern Recognit; 2000. pp. 316–323.
23. Chan TF, Vese LA. Active contours without edges. IEEE Trans Imag Process. 2001;10(2):266–277. [PubMed]
24. Kimmel R. Numerical Geometry of Images, Theory, Algorithms and Applications. New York: Springer; 2003.
25. Thevenaz P, Ruttimann UE, Unser M. A pyramid approach to subpixel registration based in intensity. IEEE Trans Imag Process. 1998 Jan;7(1) [PubMed]
26. Bland JM, Altman DG. Measuring agreement in method comparison studies. Stat Meth Med Res. 1999;8(2):135–160. [PubMed]
27. Hoffman EA, Ritman EL. Invariant total heart volume in the intact thorax. Am J Phys. 1985;249:H883–H890. [PubMed]
28. Klein GJ, Reutter BW, Budinger TF, Huesman RH. Cardiac gating of transmission data is unnecessary for attenuation compensation of double-gated emission scans. Proc. 1998 IEEE Med. Imag. Conf; 1999. pp. 2019–2022.
29. Friedman J, Van Train K, Maddahi J, Rozanski A, Prigent F, Bietendorf J, et al. Upward creep of the heart: A frequent source of false-positive reversible defects during thallium-201 stress-redistribution SPECT. J Nucl Med. 1989 Oct;30(10):1718–1722. [PubMed]