A new method to validate thoracic CT-CT deformable image registration using auto-segmented 3D anatomical landmarks

ABSTRACT

Background. Deformable image registrations are prone to errors in aligning reliable anatomically features. Consequently, identification of registration inaccuracies is important. Particularly thoracic three-dimensional (3D) computed tomography (CT)-CT image registration is challenging due to lack of contrast in lung tissue. This study aims for validation of thoracic CT-CT image registration using auto-segmented anatomically landmarks.

Material and methods. Five lymphoma patients were CT scanned three times within a period of 18 months, with the initial CT defined as the reference scan. For each patient the two successive CT scans were registered to the reference CT using three different image registration algorithms (Demons, B-spline and Affine). The image registrations were evaluated using auto-segmented anatomical landmarks (bronchial branch points) and Dice Similarity Coefficients (DSC). Deviation of corresponding bronchial landmarks were used to quantify inaccuracies in respect of both misalignment and geometric location within lungs.

Results. The median bronchial branch point deviations were 1.6, 1.1 and 4.2 (mm) for the three tested algorithms (Demons, B-spline and Affine). The maximum deviations (> 15 mm) were found within both Demons and B-spline image registrations. In the upper part of the lungs the median deviation of 1.7 (mm) was significantly different (p < 0.02) relative to the median deviations of 2.0 (mm), found in the middle and lower parts of the lungs. The DSC revealed similar registration discrepancies among the three tested algorithms, with DSC values of 0.96, 0.97 and 0.91, for respectively Demons, B-spline and the Affine algorithms.

Conclusion. Bronchial branch points were found useful to validate thoracic CT-CT image registration. Bronchial branch points identified local registration errors > 15 mm in both Demons and B-spline deformable algorithms.

The number of different medical imaging modalities has been increasing steadily since the early introduction of computed tomography (CT). In recent years positron emission tomography (PET) and magnetic resonance imaging (MRI) scanners are becoming clinical standard in oncology [1]. Consequently, multiple modalities of different three-dimensional (3D) medical image examinations are used. In order to fuse information from multiple modality images or to correlate information between single modality images, image registration are required. As an example, single modality intra-subject image registration is becoming essential in radiotherapy to detect pulmonary changes [2,3].

In general, the concept of image registration is a way to create a spatial transformation from one image into another image based on voxels properties or image features [4,5]. For registration of medical images non-rigid approaches defined as deformable image registration (DIR) [6,7] are increasingly becoming the preferable registration procedure. The DIR algorithms can be categorized as parametric and non-parametric models [5]. A parametric model, such as the B-spline, interpolates the transformation with spline models using a grid of control points [8]. A widely used non-parametric model is the Demons method [4,9]. The Demons method is an intensity-based algorithm, which uses local gradients on the static reference image to derive a displacement vector field as a transformation of the registered image. Another non-rigid method with semi-elastic properties is the Affine transformation. The Affine image registration use the rigid linear transformation combined with isotropic scaling and shearing transformations [8].

Potentially, the DIR algorithms may produce unrealistic deformations and consequently validation of registration accuracies is important. Previous studies designed for evaluation of DIR methods and algorithms have provided different methods for validation. Typically, grayscale differences, volume coincident or identifiable anatomical landmarks have been used to quantify image similarities. A validation method based on anatomical structure (landmark) appears to be a better choice for detection of unrealistic deformations [10,11]. Manual correlation of landmarks may be a time consuming process with potential interobserver variance. Consequently, an automatically process for landmark identification could benefit validation in image registration.

The aim of the present study was to validate anatomical landmarks for 3D image co-registration of thoracic CT-CT scans. The validation method was based on semi-automatically defined thoracic landmarks, defined as bronchial branch points, extracted from a segmented airway bronchial tree.

Material and methods

Five lymphoma patients were CT scanned three times within a period of 18 months (Supplementary Table I available online at http://informahealthcare.com/doi/abs/10.3109/0284186X.2015.1061215). All patients followed a protocol with a low dose full body CT scan in combination with a PET scan, as a follow-up after primary treatment (chemotherapy). The patients were scanned in supine position with the arms above the head. The CT scanner settings were 120 kVp, 35–90 mA, reconstructed with a slice thickness of 0.625 mm. The CT slices selected for image registration were from the first thoracic vertebra and 400 CT slices below covering 25 cm of the thorax including the total lung volume. For each patient the initial CT scan was defined as the reference scan (fixed CT). The two subsequent CT scans were co-registered to the reference CT volume using three different non-rigid methods: Demons, B-spline and Affine. Succeeding image registration the bronchial airways were segmented and bronchial branch point identified.

Table I. The median bronchial branch point deviations following image registration.

Variable	Median (mm)	95% CI (mm)
Patient
1	1.4	(1.2; 1.5)
2	2.2	(1.9; 2.5)
3	1.8	(1.6; 2.0)
4	2.5	(2.2; 2.8)
5	2.0	(1.7; 2.3)
Algorithm
Demons	1.6	(1.4; 1.7)
B-spline	1.1	(1.0; 1.2)
Affine	4.2	(3.8; 4.7)
Branch position
Upper lung	1.7	(1.6; 1.9)
Middle lung	2.1	(1.9; 2.3)
Lower lung	2.0	(1.8; 2.2)

Landmark definition

A seed point in trachea, defined as the lowest CT density value inside the body volume of the first CT slice, initialized a wavefront propagation described by Stephansen et al. [12] and incorporated as MATLAB functions. The wavefront propagation expanded the airway segmentation through multiple bronchial branch points using an initial adaptive grayscale threshold of -750 HU. A series of wavefronts propagated until reaching the HU threshold or a bronchial branch point. By increasing the threshold iteratively additional voxels were included within the wavefront. The wavefront expanded until leakage detection occurred, either by excessing a wavefront expansion factor or a segment expansion factor. A wavefront expansion factor of 2.8 and a segment expansion factor of 1.5 similar to the proposed values by Stephansen et al. [12] were used. A skeleton, defined as the centroid of the airway segment within each CT slice, was extracted with information of both parental segments and possible children segments. The skeleton information was used for bronchial branch point definition. A bronchial branch point was defined as the final skeleton coordinate (x,y,z), in the presence of at least two children segments. Coordinates of the bronchial branch points were automatically assigned relative to the origin of each of the CT volumes. 3D Euclidean distance deviations between equivalent bronchial branch points formed a metric of residuals following image registration. Corresponding bronchial branch points were identified manually by following the skeleton through equivalent branches. Only corresponding branch points, with identical number of child segments, were selected for evaluation of branch point deviations. Coordinates of the bronchial branch points were categorized into three superior-inferior positions within the lungs: upper, middle and lower positions, as demonstrated in Figure 1A.

Figure 1. (A) Airway segmentation for a reference image (blue) with overlay of a registered image (red). Three superior-inferior positions for the bronchial bifurcations were defined as: upper, middle and lower lung positions. The bronchial branch point in the upper right and left lung segment were classified to the upper position. The right middle lung- and the ligulae segment in the left lung were classified as middle position and the right and left lower lung segment were classified to the lower position. (B) Bronchial segment following image registration where local deformations result in misaligned airway segment.

Volume similarity

Image similarities using Dice Similarity Coefficient (DSC) [13] were defined for the overall patient volume (body) and a sub volume of the body covering the airways (lungs). The body volume was defined by creating a binary image, covering voxels above -200 HU, and post-processing to fill low-density cavities. Through indexing of the 3D binary image, into objects according to enclosed number of voxels, non-patient support equipment was removed by keeping the largest object (patient body). The lungs were defined as the sub volume of voxels within an upper threshold of -700 HU and connected to the trachea. Finally, a post-processing procedure filling detached cavities (lung parenchyma) within the lung volume, associated airway intensity values above -700 HU to the lung volume.

Registration methods

The Demons algorithm was a commercial available algorithm, Varian Medical Systems Smart Adapt Ver. 11.0, described by Thirion et al. [9]. The B-spline and Affine algorithms consisted of MATLAB functions, incorporated in procedures developed for this study. The Demons image registration features the options of both rigid and deformable registration, which allowed a rigid registration prior to the deformable registration. The volume of interest for the rigid image registrations was primary the CT scanner couch. The following volume of interest, for the deformable registration, contained the entire thoracic volume including the patient surface but excluding the scanner couch. The two MATLAB established methods, B-spline and Affine image registrations, used a rigid CT-CT procedure, prior to the final registration. The rigid registration synchronized the CT scans according to couch position, origin and reconstruction diameter. The B-spline image registration algorithm consisted of MATLAB optimization functions from Mathworks File Exchange [14]. The B-spline image registration optimization used an initial uniform grid of control points with an isotropic internal spacing corresponding to 30 mm (x,y,z). The optimizer used a Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) with a maximum number of 50 iterations and a squared sum difference for image similarity metric. The Affine image registrations used a 3D image registration procedure with the ‘Image Processing Toolbox’ for MATLAB 2013b. This algorithm used a regular step gradient descent for optimization with 100 iterations and default settings except for the initial step in the optimizer, which was set at 0.15 (named MaximumStepLength in MATLAB Image Processing Toolbox). As an image similarity metric a mean squared difference of voxel intensities were used.

Statistical analysis

Corresponding bronchial branch point were analyzed using a multi-way ANOVA comparison test adjusting for variations within algorithms, patients and superior-inferior position inside the lungs. Deviations are expressed as median (mm) with 95% confidence interval (95% CI). The DSC, expressed with 95% CI, were analyzed using a multi-way ANOVA comparison test adjusting for variations within algorithms, volumes (body and lungs) and patients. Significant differences were considered for p-values less than 0.05.

Results

A total number N = 738 corresponding set of bronchial branch points were identified (ranging from 43 to 103 in individual patients). The airway segmentation process defined the bronchial tree up to 15 cm from the main trachea bifurcation. The 3D reproducibility in branch point position were 0.6 mm (upper lung), 0.1 mm (middle lung) and 0.2 mm (lower lung). The reproducibility, defined by one standard deviation of bronchial branch points, originated from a single test CT by forcing the algorithm to segment from an initial threshold values in the range -650 to -800 HU (similar to the range used by the iterations). Table I contains a summary of bronchial branch point deviation divided into groups of patients, algorithms and superior/inferior position inside lungs. For patients the median deviation ranged from 1.4 to 2.5 mm and were significantly different between patients (p < 0.01). The tested algorithms differed significantly (p < 0.01). Median deviations ranged from 1.1 mm (B-spline), 1.6 mm (Demons) to 4.2 mm (Affine). A significant lower deviation of bronchial branch point (p = 0.02) was observed in the superior part of the lung, compared to the middle/lower lung part. The median deviations were 1.7 mm, 2.1 mm and 2.0 mm in the upper, middle and lower part of the lungs, respectively. No significant difference was found between the middle and lower sited bronchial branches. A quantile distribution of the branch point deviations (Table II), grouped by algorithms and superior-inferior positions distinguish the two locally deformable algorithms (Demons and B-spline) medians from the Affine. For the locally deformable algorithms, 95% of the bronchial branch points were registered within 5 mm. For the Affine algorithm, 95% of the branch points were registered within 10 mm. The observed maximum bronchial deviations were for the Demons 16 mm (middle lung), B-spline 15 mm (upper lung) and Affine 13 mm (upper and middle). Figure 1B demonstrates a local deformation with consequently misaligned airway structures.

Table II. Centile distribution of the bronchial branch point deviation in (mm).

	Centile	5%	25%	50%	75%	95%	Max deviation
Demons	Combined	1	1	2	2	5	16
n = 108	Upper	1	1	1	2	4	7
n = 60	Middle	1	1	2	2	12	16
n = 71	Lower	1	1	2	2	3	6
B-spline	Combined	0	1	1	2	5	15
n = 105	Upper	0	1	1	2	3	15
n = 80	Middle	0	1	1	2	6	10
n = 86	Lower	0	1	1	2	8	9
Affine	Combined	1	3	5	6	10	13
n = 102	Upper	2	3	4	5	10	13
n = 56	Middle	1	3	4	6	10	13
n = 70	Lower	1	4	5	6	9	12

Bronchial branch point deviations (total n = 738) for all five patients divided into the three algorithms and superior/inferior position inside the lungs.

The image registrations produced a higher DSC for the body, compared to the lung volume (p < 0.01). Additionally, significant different DSCs were found between the algorithms (p < 0.01) (Table III). The B-spline image registrations yield the maximum DSC, whereas the Affine transformation produced the lowest accuracy according to the DSC. No significant differences of the DSCs among the five patients were found. The mean lung surface differences, using Euclidian distance metric with 95% CI, were for the Demons algorithm 2.2 mm (2.1; 2.2) mm, B-spline 2.1 mm (2.0; 2.2) mm and the Affine 3.1 mm (2.5; 3.7) mm.

Table III. The body and lung DSC following image registration.

Dice coefficients	Body volume 95% CI	Lung volume 95% CI
Demons	0.99 (0.98; 0.99)	0.96 (0.96; 0.96)
B-spline	0.99 (0.99; 0.99)	0.97 (0.97; 0.98)
Affine	0.96 (0.96; 0.97)	0.91 (0.88; 0.95)

DSC calculated for all patients and scans separated into two volumes and the three algorithms.

Discussion

The present study aimed to evaluate different DIR methods, based on automatically segmented bronchial branch points using corresponding anatomical landmarks. Despite the automatically process of detecting bronchial branch points, identification of corresponding branch points revealed a time consuming task. Primary due to classification of superior/inferior position and the different number of sub-bronchial segmentations, as illustrated in Figure 1. The image material for this study consisted of low dose CT scans, which influenced on the image quality compared to a normal CT. The low dose CT may limit the bronchial tree segmentation in each lung volume. Normal dose or even better high resolution diagnostic CT would most likely give an enhanced segmentation of the bronchial tree. Excluding the most distant parts of the bronchial tree would most likely eliminate the tedious manual task of checking for corresponding landmarks.

The present study revealed discrepancies of image registrations between the three tested algorithms. Despite a good median accuracy for image registration with local deformations (Demons and B-spline), the maximum branch point deviations showed a potential risk of large local deformation errors of approximately 15 mm. These maximum bronchial branch deviations appeared in the middle lung for the Demons algorithm and in the upper lung for the B-spline algorithm. Particularly the bronchial branches in the upper lung demonstrated the lowest median deviation of 1.7 mm, but for the B-spline algorithm a misalignment of 15 mm as well. Identification of such a registration error would be crucial for high dose thoracic radiotherapy with the bronchi as an organ at risk [15]. The bronchial branch point deviations for the Affine method was expected to be higher relative to the deformable algorithms due the lack of ability local deformations within the image set. A study by Senthi et al. [16], comparing a rigid and deformable algorithm defined as contour differences, demonstrated improved registration accuracy by 3 mm when using the deformable registrations. In the present study a similar improvement of the median bronchial branch points were observed between the Affine and B-spline algorithm. Regarding the Demons and B-spline algorithms, both the DSC and the deviation between the bronchial branch points maintained a lower deviation for the B-spline. A general lower bronchial branch deviation for a B-spline algorithm was in agreement with a study by Varadhan et al. [5], which specified increased accuracy for the B-spline method compared to a Demons algorithm. The work from Varadhan et al. dedicated a smoother deformations field for the B-spline method relative to a Demons algorithm as a reason for a more reliable deformation.

Landmarks as the bronchial branch point in the present study provided a reasonable method to validate image registrations comparable to the DSC. Using the DSC as a metric for image registration accuracy provided a significant higher image similarity for the body volume relative to the lungs. However, this might be related to the mathematical definition of the DSC, consequently the DSC may not be appropriate for large volume associations. The DSC confirmed significant differences among all three tested image registration algorithms, in the favor of the B-spline algorithm. In contrast to the DSC quantity landmark identifications provided additional information related to local misaligned image registrations.

A main limitation, using the bronchial tree as landmark identifications, is the ability to segment the airway tree and thus define the bronchial branch points. This could be of concern for patients with partial airway obstruction or severe lung density changes due to radiotherapy [3,17]. Chow et al. specified the requirement for DIR for proper detection and follow-up of lung injury succeeding chemoradiotherapy. For images decomposed by density changes, the segmented airways as well as a map of corresponding bronchial branch point have prospect for usage within DIRs. A study by Vasquez Osorio et al. [18] revealed CT-MR image registration, established on automatic vessel segmentation of the liver in the two image modalities. Their method has the potential for adaption in thoracic CT-CT co-registration, where local deformation within the lung volume could be improved using segmented airways.

In conclusion, the present study demonstrated a lower median overall registration error when using the B-spline algorithm compared to the Demons and Affine algorithm. However, local registrations errors of approximately 15 mm were identified for both the Demons and B-spline algorithms using the bronchial branch points for validation of image registration.

Supplementary material available online

Supplementary Table I available online at http://informahealthcare.com/doi/abs/10.3109/0284186X.2015.1061215.

Declaration of interest: The work is supported by DAKFO – Danish Graduate School in Clinical Oncology and CIRRO – The Lundbeck Foundation Center for Interventional Research in Radiation Oncology and The Danish Council for Strategic Research. No conflicts of interest to declare. The authors alone are responsible for the content and writing of the paper.