A reference dataset for deformable image registration spatial accuracy evaluation using the COPDgene study archive

The COPDgene study protocol indicates the acquisition of dual inspiratory and expiratory chest CT scans for all study subjects to identify and quantify airway abnormalities, emphysema, and expiratory air trapping (Regan et al 2010). The inspiratory BH-CT is to provide thorough assessment of small airway wall thickness and emphysema, whereas the expiratory BH-CT is performed to assess for air trapping. To facilitate image-based measurements to assess small airway wall thickness and abnormalities, images are reconstructed using a high-resolution recon kernel. The algorithm used is vendor-specific (i.e., GE: BONE, SIEMENS: B46f, PHILIPS: Detail (D)), and represents a balance between enhancement of image resolution via higher spatial frequency cutoff, and increase in image noise. Since each algorithm previously has been selected as meeting the COPDgene study criteria, it is not expected that there will be CT vendor specific issues with the quantitative image analysis.

In this work, ten BH-CT image pairs were randomly selected from the COPDgene study cases. Each patient had received CT imaging of the entire thorax in the supine position at normal expiration and maximum effort full inspiration. The CT imaging was performed with a GE VCT 64-slice scanner (GE Healthcare Technologies, Waukesha, WI) with a pitch of 1.375 mm, speed of 13.75 mm per rotation, 120 kV_p, 0.5 s per rotation, 400 mA per rotation for inhale BH, and 100 mA per rotation for exhale BH. The images used in this study were reconstructed with high-resolution (BONE) recon kernel, with lung diameter setting the field of view, and with contiguous 2.5 mm slice spacing. Each image set was reconstructed to image dimensions 512 × 512 × N, (N ∈ [102135]), with in-plane voxel dimensions ranging from (0.590 × 0.590)–(0.652 × 0.652) mm². The image characteristics of the ten COPD cases utilized in this study are given in table 1.

Additional quantitative analyses were performed characterizing lung volume change and relative image noise between BH-CT pairs. To approximate lung volume from each BH-CT image, the set of voxels with CT number in the range [−1000 −200] was initially selected. Two-dimensional region growing was applied slice-by-slice in the transverse plain to remove background objects defined by connectivity to the image border (Castillo et al 2009). A voxel erosion and subsequent growing technique was applied to remove extraneous voxels. The final segmented regions were visually assessed and manually modified as necessary. The volume of each segmented lung region was calculated and the difference between corresponding iBH- and eBH-CT determined for each case. Volume changes were determined as absolute magnitude change (in ml) and percent change relative to the segmented eBH-CT lung volume.

Estimates of relative noise between BH-CT pairs were determined by characterization of image intensity distributions within manually placed regions of interest (ROIs). For each case, an expert in thoracic imaging placed a fixed 1.5 cm diameter sphere on corresponding iBH- and eBH-CT within the aorta, approximately mid-position relative to the superior—inferior (SI) lung volume extent (figure 3). The choice of aorta for ROI placement was based on the availability of centrally located, uniform measurements within the greater lung volume of interest, and the expected constant blood image intensity between inhale and exhale breathing states. Care was taken to avoid regions with prominent shadowing effects by calcifications in the immediate vicinity of the major vessel. For each ROI, the coefficient of variation (COV) was determined as standard deviation (SD) divided by mean CT number within the sphere. The ratio of eBH-to-iBH COVs was determined for each case as a surrogate measure of the relative noise content within corresponding image pairs.

Zoom In Zoom Out Reset image size

A MATLAB-based software interface named APRIL (Assisted Point Registration of Internal Landmarks), which has been previously described (Castillo et al 2009), was utilized to facilitate manual selection of landmark feature pairs between the volumetric eBH- and iBH-CT images. The software incorporates basic image processing and visualization tools that are common to most commercial treatment planning systems, including independent window and level settings for all displays, digital magnification and pan, visualization of equivalent voxel positions in the principal orthogonal plains, transparent volume rendering, and interactive tools for segmentation of lung voxels from the image data. To determine corresponding features, the user initially must manually delineate a candidate voxel in the designated 'source' image via mouse click over the desired position. The target image voxel corresponding to the same underlying anatomy is then similarly selected to complete the feature pairing. The process is repeated until the desired sample size and uniformity of distribution have been achieved. No implanted fiducials or added contrast agents were used to aid in the selection of landmark features, which typically included vessel and bronchial bifurcations.

To assist in localization of the anatomic features within the target image space, the software provides an optional feature localization tool based on normalized cross-correlation of the local feature neighborhood depicted in the source image with a larger, user-adjustable target window in the eBH-CT (Castillo et al 2009, Lewis 1995, Gonzalez and Woods 2008). The image content within the feature neighborhoods is converted to a binary mask according to a user-defined intensity threshold, such that the cross correlation is computed only over the local structural content. The intensity threshold controls the level of structural detail included in the feature masks. The target voxel representing the maximum of the three-dimensional correlation function is highlighted with a crosshair within the target display, and represents an estimate of the feature correspondence. Previously we found that localization estimates based on the cross-correlation coefficient were suitable for expediting the manual feature pair registration process using treatment planning 4D CT. However, preliminary investigation as part of this work demonstrated little to no utility for use with the COPDgene imaging data. Most notably, this is due to the drastic deformation of prominent vasculature structures commonly identified as anatomic landmark features. In addition, the relatively large displacements, which require excessively large search windows for both images, further increase the probability that the correlation function maximum does not correspond to accurate approximation of the underlying physical displacement.

As an alternative, to facilitate maximally efficient workflow, we have upgraded the APRIL software to include an additional feature localization capability, which incorporates the current set of registered feature pairs into moving least squares (MLS) interpolation applied to the current source feature position to determine the corresponding target estimate. MLS is a highly generic and versatile tool for approximating an unknown function by fitting polynomials to function samples given at uniform or non-uniform locations (Bos and Salkauskas 1989, Lancaster and Salkauskas 1981). Though most commonly employed to reconstruct surfaces from noisy, unstructured point cloud data, MLS has recently found utility in DIR (Castillo et al 2010a, 2012a, Schaefer et al 2006). For our purposes, the MLS function operates on the current source voxel position $\vec x$, for which we would like to obtain an estimate of the corresponding target position. An affine function $A_{\vec x}$ is determined by minimizing the expression:

$$\begin{equation} \min \sum\limits_i {\omega _i \left\| {A_{\vec x} \left( {\vec s_i } \right) - \vec t_i } \right\|^2 } , \end{equation} \tag{ 1.1 }$$

where $\vec s_i$ and $\vec t_i$ are the ith source and target feature pair positions, respectively. The scalar weights ω_i are of the form:

$$\begin{equation} \omega _i = \frac{1}{{\left\| {\vec s_i - \vec x} \,\right\|^2 + \varepsilon }}\,\varepsilon > 0. \end{equation} \tag{ 1.2 }$$

In this way, each newly selected feature-pair provides additional contribution to estimates of subsequent feature locations. While an initial set of registered feature pairs is required before the MLS estimation procedure provides non-trivial support, in practice we have found that approximately four input pairs are necessary before the estimates contribute to improved target feature localization. This significantly reduces the time spent navigating the target volume space, and is independent of the magnitude displacements or image content.

In this work, both cross-correlation and MLS-based estimation procedures were available to all users. Cross correlation and MLS-based estimates of target features are presented as a crosshair in the target image. In practice, multiple cross correlations varying both the feature neighborhood dimension and/or intensity threshold may be performed for a given feature. The corresponding target estimate will vary according to the particular search criteria. Similarly, the MLS-based estimates will vary according to the current set of previously matched feature pairs. As more features are matched, the estimates become more robust. However, for both estimation procedures, the user ultimately must designate the feature correspondence via mouse click within the target image. The reader visually determines every point; there is not a mechanism to automatically accept a feature position based on either automated assistance operation.

An expert in thoracic imaging, beginning at the apex of the lung, systematically selected the source feature points on the ten image pairs. The manual registration process consists of initial feature localization via mouse click within the source image. The feature position is highlighted with a crosshair, and the user is prompted for confirmation. Following confirmation, the optional estimation tools become enabled to provide computer assistance for target feature localization based on either cross correlation or MLS operations. The reader is free to navigate the target volume along each of the three principal orthogonal axes, with the option to display the source feature coordinate via crosshair within the target image space for reference. The matched target feature position is identified via mouse click and indicated with a crosshair. Optionally, the user may overlay the discretized image voxel grid for assistance in precise localization of the intended anatomic feature. Following confirmation, the coordinate position is stored. At any time, the reader is free to enable a review mode to scroll through the current set of matched feature pairs. The reader may specify a feature pair index, or sequentially review the current set via scroll bar. Previously matched features are presented in both source and target displays, with crosshair centered according to the current magnification setting in each window. Previously accepted source and target positions can be modified as necessary.

In this study, the iBH-CT was designated as the 'source' image, which served as the primary dataset within which anatomic features were initially selected. This was due to the better image quality (table 3) and the higher contrast of vessels relative to the inflated lung background. Points were selected with an initial goal of >5 for each lung per axial image slice. This approach ensured the collection of ⩾447 validation point pairs for each case distributed throughout the lungs. Following feature selection for a given case, all landmark pairs were visually reviewed by the primary reader a second time and the locations adjusted on either image as necessary. The verification step was a required part of the initial registration process performed by the primary reader.

A subset of landmark points was re-registered by the primary reader, to estimate intra-observer variance in target localization, and by secondary readers, to estimate the inter-observer variance. For each case, the APRIL software was used to extract a random sample of 150 feature pairs from the corresponding complete set generated by the primary reader. For each of the 150 sampled pairs, the reader was presented with the primary iBH-CT feature position via crosshair within the source image space. The re-registration process consists of locating the corresponding point within the eBH-CT. All feature estimation and visualization tools were available in both primary and secondary registration processes. For all secondary readings, observers were blinded to the primary target selections; all secondary readings were performed independently and without prior knowledge of the primary point matches. Observer variation in target point selections for a given image pair were estimated based on the three-dimensional Euclidean distance in units of millimeters between the primary and secondary target positions. The mean distance, with SD, between corresponding targets was determined over the range of subsampled data. Inter- and intra-observer variance in feature localization was determined on a case-by-case basis.

Proper and efficient utilization of the assisted landmark registration software is a learning process that requires a certain degree of practice and experience to realize. In general, care should be taken to ensure that variations in target point selection both within and between users are true reflections of variability in individual judgment and selection criteria, rather than reflections of variation in time-specific aptitude or experience with the registration software. To address this concern, all observers were required to complete a minimum training set of 1000 independent feature-pair registrations, using image pairs not included in the study, prior to repeat measurements over the reference data. Additionally, potential degradation in CT image quality due to noise or image acquisition/reconstruction artifacts is likely to have some impact on the quantity of prominent features an observer is able to match between images. Rather than filter the image data, we will accept all clinically acquired cases, such that the reference library will accurately reflect those images that are likely to be encountered in routine clinical application.

Estimated lung volumes, with corresponding magnitude and relative changes are shown in table 2 for the set of BH-CT pairs included in this study. Magnitude volume change between image pairs varied considerably among cases, with range from 974–2819 ml. To provide context for these measurements, we previously reported estimates for resting tidal volumes derived from seven treatment planning 4D CTs using an analogous image segmentation-based approach, which yielded measured volume changes in the range 462–926 ml (Castillo et al 2010b). The extreme discrepancy is expected, and reflective of the markedly different clinical scenarios investigated in each study (i.e., tidal breathing versus maximum effort BH). Relative changes between BH-CT pairs ranged from 25%–106%, indicating in some instances greater than two-fold increase in lung volume between corresponding inhale and exhale BHs. Figure 4 provides visual context for these measurements, illustrating an example image pair (case 5) for which the relative change was >100%. Maximum distance between the inferior contoured lung boundaries is on the order of 40 image voxels.

Zoom In Zoom Out Reset image size

Estimates of the relative noise between corresponding iBH- and eBH-CT image pairs are summarized in table 3, in terms of the COV ratio within manually placed spherical ROIs (figure 3). The measurements indicate relative noise greater than the theoretical (2:1) given by consideration of the relative acquisition tube currents alone. However, the high-resolution BONE reconstruction filter, which will increase the noise in greater relative proportion in the exhale image, also contributed to the increase in the overall comparison ratio (Goldman 2007). Figure 3 shows representative transverse and coronal views through an example iBH- and eBH-CT image pair (case 7), with cross-section of the spherical ROI also indicated. The example views illustrate the obvious increase in mottled appearance of soft-tissue structures depicted in the eBH-CT.

The reference landmark characteristics for all cases included in this study are summarized in table 1. A total of 7298 anatomic landmark features were manually paired between the 10 sets of images, while the registered feature pairs per case ranged from 447 to 1172. Approximately 10–12 h were required to complete the primary feature pair dataset for each case. For secondary readers, approximately 6–8 h were required to re-register the 150 sampled source features onto the eBH-CT. Average three-dimensional Euclidean landmark displacements varied substantially among cases, ranging from 12.29 (SD: 6.39) to 30.90 (SD: 14.05) mm. Average magnitude displacements in component right–left, anterior–posterior (AP), and SI directions ranged from 2.11 (SD: 1.67) to 6.06 (SD: 3.97), 6.11 (SD: 3.50) to 24.57 (SD: 6.25), and 6.20 (SD: 5.45) to 21.20 (SD: 14.75) mm, respectively. The component displacement magnitudes clearly indicate that the largest displacements occur along the AP axis, which is unlike the previously reported 4D CT reference images, for which the most significant tidal motion occurred along the SI axis (Castillo et al 2009). Repeat registration of uniformly sampled subsets of 150 landmarks for each case yielded estimates of observer localization error, which ranged in average from 0.58 (SD: 0.87) to 1.06 (SD: 2.38) mm for each case. Figure 5 graphically summarizes each individual observer's localization errors per case. Figure 6 shows the set of reference displacement vector field projections for each case. Transparent isosurface renderings of the maximum effort iBH-CT lung fields are shown overlain in oblique projection. The color scale shown at right is fixed for each projection, illustrating the substantial variability in motion field characteristics among the ten cases, which is due in part to the variable extent of their COPD disease.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Of notable concern for feature-based DIR evaluation is the density of anatomical landmarks with respect to the organ surface. Analyses based on samples that are biased toward high contrast, centrally located structures are likely to misrepresent algorithm performance at the periphery where the landmark density is low. To characterize the overall feature pair distribution with respect to the lung surface, the minimum distance to the segmented lung contour was determined for each landmark in the set of eBH-CTs. The cumulative distribution of distance measurements was compiled for all cases and shown in figure 7. For the dataset generated in this study, approximately 11% of the reference features are within 5 mm of the contoured exhale lung surface, while approximately 46% are within 1 cm. The measurement process is illustrated in figure 7 for an example transverse section in which 20 features are identified ranging in distance to surface from 3–37 mm.

Zoom In Zoom Out Reset image size