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 defi ned 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 intensitybased algorithm, which uses local gradients on the static reference image to derive a displacement vector fi eld as a transformation of the registered image. Another non-rigid method with semi-elastic properties is the Affi ne transformation. The Affi ne 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 identifi able 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 identifi cation could benefi t 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 defi ned thoracic landmarks, defi ned 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 fi rst 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 defi ned as the reference scan (fi xed CT). The two subsequent CT scans were co-registered to the reference CT volume using three different non-rigid methods: Demons, B-spline and Affi ne. Succeeding image registration the bronchial airways were segmented and bronchial branch point identifi ed.

Landmark defi nition
A seed point in trachea, defi ned as the lowest CT density value inside the body volume of the fi rst CT slice, initialized a wavefront propagation described by Stephansen et al. [12] and incorporated as MAT-LAB 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, defi ned 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 defi nition. A bronchial branch point was defi ned as the fi nal 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 identifi ed 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.

Volume similarity
Image similarities using Dice Similarity Coeffi cient (DSC) [13] were defi ned for the overall patient volume (body) and a sub volume of the body covering the airways (lungs). The body volume was defi ned by creating a binary image, covering voxels above -200 HU, and post-processing to fi ll 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 defi ned as the sub volume of voxels within an upper threshold of -700 HU and connected to the trachea. Finally, a post-processing procedure fi lling 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 Affi ne 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 Affi ne image registrations, used a rigid CT-CT procedure, prior to the fi nal 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 Affi ne 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 MaximumStep-Length 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% confi dence 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. Signifi cant differences were considered for p-values less than 0.05.

Results
A total number N ϭ 738 corresponding set of bronchial branch points were identifi ed (ranging from 43 to 103 in individual patients). The airway segmentation process defi ned 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, defi ned 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 signifi cantly different between patients (p Ͻ 0.01). The tested algorithms differed signifi cantly (p Ͻ 0.01). Median Bronchial branch point deviations (total n ϭ 738) for all fi ve patients divided into the three algorithms and superior/inferior position inside the lungs.  (Table II), grouped by algorithms and superior-inferior positions distinguish the two locally deformable algorithms (Demons and B-spline) medians from the Affi ne. For the locally deformable algorithms, 95% of the bronchial branch points were registered within 5 mm. For the Affi ne 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 Affi ne 13 mm (upper and middle). Figure 1B demonstrates a local deformation with consequently misaligned airway structures.
The image registrations produced a higher DSC for the body, compared to the lung volume (p Ͻ 0.01). Additionally, signifi cant different DSCs were found between the algorithms (p Ͻ 0.01) ( Table III). The B-spline image registrations yield the maximum DSC, whereas the Affi ne transformation produced the lowest accuracy according to the DSC. No signifi cant differences of the DSCs among the fi ve 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 Affi ne 3.1 mm (2.5; 3.7) mm.

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, identifi cation of corresponding branch points revealed a time consuming task. Primary due to classifi cation 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 infl uenced 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. Identifi cation 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 Affi ne 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 defi ned 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 Affi ne 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 specifi ed increased accuracy for the B-spline method compared to a Demons algorithm. The work from Varadhan et al. dedicated a smoother deformations fi eld 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 signifi cant higher image similarity for the body volume relative to the lungs. However, this might be related to the mathematical defi nition of the DSC, consequently the DSC may not be appropriate for large volume associations. The DSC confi rmed signifi cant differences among all three tested image registration algorithms, in the favor of the B-spline algorithm. In contrast to the DSC quantity landmark identifi cations provided additional information related to local misaligned image registrations.
A main limitation, using the bronchial tree as landmark identifi cations, is the ability to segment the airway tree and thus defi ne 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. spec-ifi ed 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 Affi ne algorithm. However, local registrations errors of approximately 15 mm were identifi ed for both the Demons and B-spline algorithms using the bronchial branch points for validation of image registration.