Development of image registration methods for aided diagnosis and monitoring in computed tomography

Abstract

Extended Abstract The term interstitial lung disease (ILD) refers to more than 200 chronic lung disorders that are classified to the same group because of their similar clinical, radiological, physiologic and pathologic features [1]. Accurate chest CT quantification of ILD extent is crucial for patient management, since no robust biomarkers for monitoring disease progression exist [2]. Chest computed tomography, (CT) imaging is a powerful modality for detection, diagnosis, and follow up of interstitial lung Disease. CT bases assessment of ILD extent estimation is performed by several semi-quantitative scoring systems, estimating extent, utilizing high-resolution CT (HRCT) protocols [3], which are deprived of information regarding localization of disease, demonstrating moderate inter- and intraobserver agreement [4]. Thus, ILD extent estimation (quantification) and progression by means of advanced image analysis methods, is emerging offering increased accuracy and reproducibility. [5][6][7]. In the frame of quantitative image-based follow-up ILD monitoring and response to therapy in CT, image registration methods have the important role to ensure that any measured volume change between follow-up scans is caused by ILD patterns change and not by patient’s breathing or positioning during CT scanning. Multiresolution non-rigid registration, capable for capturing local lung tissue deformations, accounts for a commonly used approach for lung CT registration utilized 2 in all resolution levels (fully non rigid schemes), or in combination with rigid transforms applied in low resolution levels (hybrid schemes). Up-to-now only one study has reported on ILD follow-up monitoring, utilizing deformable image registration (DIR) [8]. However, registration scheme parameters were adopted rather than obtained from a systematic analysis taking in to account their effect to registration accuracy. Evaluation methods for DIR accuracy assessment in CT of the thorax has receive considerable attention in the frame of adaptive radiotherapy [9], however remains an open issue [10][11]. The current thesis addresses, selection of optimal registration schemes in case of ILD follow-up in CT. The current thesis addresses selection of suitable multiresolution registration schemes for CT based follow up analysis of ILD with the following specific objectives: Development of an evaluation methodology to obtain high accuracy schemes Analyzing the effect of scheme components on both registration accuracy and time efficiency Investigation of the effect of mass preserving cost function to registration accuracy and time efficiency. Investigation of the effect of regularization term in artificially generated follow up data. The evaluation methodology, considers two stages: the first stage utilizes artificially warped ILD follow-up data to identify candidate registration schemes, while in the second stage the performance of candidate schemes, identified in the first stage, are verified with clinical follow-up data . At the first stage, the determinant of the Jacobian matrix of the each voxel of deformation field is the main tool to identify schemes including folding areas to be excluded from subsequent analysis. The performance of the remaining schemes is assessed and ranked in terms of their displacement error in two anatomical regions, i.e. Normal Lung Parenchyma (NLP) and ILD affect regions, by means of Euclidean distance of homologus points between baseline and register follow up pair. Statistical analysis was performed to select near optimal schemes, considering both NLP and ILD ranking lists of candidate schemes. In the second stage, evaluation of registration accuracy of candidate registration schemes is verified on their clinical follow up volumetric scans using the displacement error as well. Finally selected registration 3 schemes, was also performed by statistical analysis, considering both NLP and ILD ranking lists of near optimal schemes. A clinical dataset consisting of 10 pairs of CT scans corresponding to 10 patients diagnosed with ILD secondary to connective tissue diseases, radiologically manifested with ground glass and reticular patterns, at two different times, abstaining in time approximately two years was acquired. The total extent of ILD presence in these 10 scans ranges from 5% to 80% (mean value: 32%). Artificial warped data [12] is introduced to ensure that registration error is caused by registration algorithm alone and not by intrinsic data variability. Artificially warped data were generated using a single level non rigid thin plate kernel spline model and applied to baseline lung segments, resulting in artificial follow up lung segment simulating realistically lung deformations, preserving size and slice thickness of the original data set. The basic components of a typical registration schemes are the transformation model, the cost function the optimizer and the type of pyramid used. A total of 128 registration schemes was generated, by considering the following combinations of components: Four (4) transformations: Euler Transform - ET, Similarity Transform - SM, Affine Transform - AT, and 3rd order B-spline Transform - BST, applied to the first 3 resolution levels of the pyramid in order to obtain a coarse initial alignment, while the 3rd order B-Spline transform was utilized for the 4th resolution level, corresponding to the highest image resolution, in all schemes generated. Although fully non-rigid transformation models seem the natural choice for lung field registration due to the elastic nature of the lung tissue, hybrid schemes (including rigid and non rigid transformation models in different resolution levels), are also been recently proposed. As the optimizer accounts as for a critical component of the registration process, 4 gradient decent optimizers were utilized for the optimization step: SGD, regular step gradient decent - RSGD, adaptive stochastic gradient decent - ASGD and finite difference gradient decent - FDGD. Additionally, two 2 types of pyramids were considered: Gaussian Pyramid - GP, that applies smoothing and downsampling by a factor of 2 in all three dimensions, and Recursive Pyramid – RP, that applies no smoothing but only downsampling by a factor of 2 in all three dimensions. Finally, 4 different cost functions were considered (Sum of Square Difference – SSD, Normalized Correlation Coefficient – NCC, Mutual Information – MI, and Normalized Mutual Information – NMI). Taking 4 in to account, the susceptibility of SSD cost function to air quantity inside the lungs [13][14], due to varying breathing phase, a mass preserving variant of SSD was also considered in a separate experiment . All registration schemes evaluated in this study, are intensity based and utilize the multiresolution approach to avoid local minimum traps and speed-up calculations, using Elastix version 4.5, based on the open source software Insight Toolkit (ITK) version 4.0. Additi

Εκτεταμένη Περίληψη Ο όρος διάμεση πνευμονοπάθεια (ΔΠΝ) αναφέρεται σε περισσότερες από 200 χρόνιες παθήσεις των πνευμόνων που κατατάσσονται στην ίδια ομάδα, λόγω των παρόμοιων κλινικών, ακτινολογικών, και παθοφυσιολογικών χαρακτηριστικών τους [1]. Η ακριβής ποσοτικοποίηση της ΔΠΝ στην υπολογιστική τομογραφία- ΥΤ (Computed Tomography) θώρακος είναι ζωτικής σημασίας για τη διαχείριση των ασθενών, αφού δεν υπάρχουν ευρέως αποδεκτοί βιοδείκτες για την παρακολούθηση της εξέλιξης της νόσου [2]. Η ΥΤ θώρακος, αποτελεί ισχυρό εργαλείο ανίχνευσης, διάγνωσης και παρακολούθηση διάμεσης πνευμονοπάθειας. Η ακτινολογική εκτίμηση της έκτασης της ΔΠΝ με βάση την ΥΤ πραγματοποιείται με χρήση διάφορων ημι-ποσοτικών κλιμάκων βαθμολόγησης καθώς και χρήση δισ-διαστατου πρωτοκόλλου υψηλής ανάλυσης (High Resolution CT) [3]. Χρήση ημι-ποσοτικών κλιμάκων, ως μεθόδου εκτίμησης έκτασης της νόσου όμως υστερείται πληροφοριών σχετικά με τον χωροταξικό εντοπισμό της νόσου, καταδεικνύοντας μέτρια δια- και ενδο- συμφωνία παρατηρητή [4].Αποτέλεσμα των ανωτέρω αποτελεί η τρέχουσα τάση αξιοποίησης προηγμένων μεθόδων ανάλυσης εικόνας για ποσοτικό προσδιορισμό της έκτασης της ΔΠΝ, οι οποίες χαρακτηρίζονται από αυξημένη ακρίβεια και επαναληψιμότητα. [5] [6] [7]. Η ποσοτικοποίηση και παρακολούθηση της έκτασης της ΔΠΝ βάσει CT, στο πλαίσιο διαχείρισης ασθενούς σχετικά με την ανταπόκριση στη θεραπεία, αναδεικνύει το ρόλο των μεθόδων Αντιστοίχισης εικόνας. Ο ρόλος αυτός συνίσταται στη διασφάλιση ότι κάθε μετρούμενη μεταβολή της έκτασης της νόσου μεταξύ διαδοχικών χρονικά σαρώσεων οφείλεται στη μεταβολή της έκταση των απεικονιστικών προτύπων της ΔΠΝ και όχι σε παράγοντες όπως η φάση αναπνοής του ασθενούς ή τοποθέτηση κατά τη διάρκεια της τομογραφικής λήψης. 8 Η μέθοδος της εύκαμπτης Αντιστοίχισης εικόνων (deformable image registration-dir) πολλαπλών κλιμάκων ανάλυσης, αντιπροσωπεύει μια ευρέως χρησιμοποιούμενη προσέγγιση για την Αντιστοίχιση των πνευμονικών πεδίων σε ΥΤ, για τον προσδιορισμό μετασχηματισμών εικόνων, ικανών για την προσομοίωση των τοπικών παραμορφώσεων των ιστών του πνεύμονα (πνευμονικού παρεγχύματος). H μέθοδος DIR μπορεί να χρησιμοποιηθεί σε όλα τα επίπεδα ανάλυσης συστήματος πολλαπλής κλίμακας, οδηγώντας σε πλήρως εύκαμπτα συστήματα πολλαπλών κλιμάκων ανάλυσης (fully deformable multiresolution registration schemes), ή σε συνδυασμό με άκαμπτους (rigid) μετασχηματισμούς εικόνας που αξιοποιούνται σε επίπεδα ανάλυσης χαμηλής κλίμακας, οδηγώντας σε υβριδικά συστήματα εύκαμπτης Αντιστοίχισης εικόνων πολλαπλών κλιμάκων ανάλυσης (hybrid deformable multiresolution registration schemes). Μέχρι σήμερα Αντιστοίχισης εικόνων μόνο μία μελέτη έχει αναφερθεί στη βιβλιογραφία σχετικά με την χρήση συστήματος εύκαμπτης Αντιστοίχισης εικόνων πολλαπλών κλιμάκων ανάλυσης στα πλαίσια συστήματος ποσοτικοποίησης έκτασης και παρακολούθησης της ΔΠΝ, [8]. Ωστόσο, οι παράμετροι του συστήματος Αντιστοίχισης που χρησιμοποιήθηκαν υιοθετήθηκαν και δεν αναλύθηκαν συστηματικά. Οι μέθοδοι αξιολόγησης απόδοσης των DIR συστημάτων για την περίπτωση της υπολογιστικής τομογραφίας θώρακος αν και έχουν συγκεντρώσει αρκετό ενδιαφέρον ειδικά στην περίπτωση της ακτινοθεραπείας [9], εξακολουθούν να παραμένουν ανοιχτό πρόβλημα [10][11]. H παρούσα διατριβή αποτελεί μία πρώτη προσπάθεια εστίασης στην επίδραση των επιμέρους δομικών στοιχείων (παραμέτρων) συστημάτων εύκαμπτης Αντιστοίχισης εικόνων) πολλαπλών κλιμάκων ανάλυσης, στο πλαίσιο επιλογής των βέλτιστων συστημάτων Αντιστοίχισης στην περίπτωση παρακολούθησης ΔΠΝ σε ΥΤ, με τους ακόλουθους ειδικούς στόχους: Ανάπτυξη μεθοδολογίας αξιολόγησης ακρίβειας συστημάτων εύκαμπτης Αντιστοίχισης εικόνων πολλαπλών κλιμάκων ανάλυσης Ανάλυση της επίδρασης των δομικών στοιχείων του συστήματος Αντιστοίχισης τόσο ως προς την ακρίβεια Αντιστοίχισης όσο και ως προς τον υπολογιστικό χρόνο (κόστος) 9 Διερεύνηση της επίδρασης του όρου διατήρησης μάζας σε συγκεκριμένη συνάρτηση κόστους ως προς την ακρίβεια Αντιστοίχισης όσο και ως προς τον υπολογιστικό χρόνο (κόστος) Διερεύνηση της επίδρασης του όρου ομαλοποίησης (regularization ) ως προς την ακρίβεια Αντιστοίχισης όσο και ως προς τον υπολογιστικό χρόνο (κόστος), με χρήση σε τεχνητών δεδομένων παρακολούθησης, που προσομοιάζουν την ακτινολογική εμφάνιση πνευμονικών πεδίων. Η μεθοδολογία αξιολόγησης, θεωρεί δύο στάδια: το πρώτο στάδιο χρησιμοποιεί τεχνητά δεδομένα παρακολούθησης ILD μέσω στρέβλωσης των αρχικών εικόνων [12] για τον εντοπισμό υποψήφιων συστημάτων αντιστοίχισης, ενώ στο δεύτερο στάδιο οι επιδόσεις των υποψηφίων συστημάτων, τα οποία προσδιορίζονται στο πρώτο στάδιο, επαληθεύονται με κλινικά δεδομένα παρακολούθησης ΔΠΝ. Στο πρώτο στάδιο, η ορίζουσα της Ιακωβιανής μήτρας (Determinant of the Jacobian Matrix-DETJ) του κάθε ογκοστοιχείου (voxel) του πεδίου παραμόρφωσης (deformation field) είναι το βασικό εργαλείο για τον εντοπισμό συστημάτων που περιλαμβάνουν περιοχές αναδίπλωσης (folding areas or singularities). Τα συστήματα αυτά που συμπεριλαμβάνουν τέτοιες περιοχές αποκλείονται από την περαιτέρω ανάλυση. Η απόδοση των υπόλοιπων υποψήφιων συστημάτων αξιολογείται και κατατάσσεται με όρους σφάλματος μετατόπισης τους (displacement error) ως προς σε δύο ανατομικές περιοχές, δηλαδή στην περιοχή φυσιολογικού πνευμονικού παρεγχύματος (Normal Lung Parenchyma - NLP) και στην παθολογική περιοχή (ILD), μέσω της Ευκλείδειας απόστασης των ομόλογων σημείων μεταξύ της αρχικής εικόνας και τεχνητής εικόνας παρακολούθησης. Η επιλογή των υποψήφιων συστημάτων έγινε επί τη βάσει κριτηρίου στατιστικής σημαντικότητας, λαμβάνοντας υπόψη και τις δύο λίστες κατάταξης των υποψηφίων συστημάτων ως προς το σφάλμα μετατόπισης ομόλογων διακριτικών ανατομικών σημείων (anatomical landmarks) στις περιοχές NLP και ILD. Στο δεύτερο στάδιο, η αξιολόγηση της ακρίβειας των υποψήφιων συστημάτων Αντιστοίχισης επαληθεύεται σε κλινικά δεδομένα παρακολούθησης χρησιμοποιώντας επίσης το σφάλμα μετατόπισης ομόλογων διακριτικών ανατομικών σημείων, στις ανωτέρω δύο ανατομικές περιοχές. Τέλος, η επιλογή των βέλτιστων συστημάτων έγινε επίσης επί τη 10 βάσει κριτηρίου στατιστικής σημαντικότητας, λαμβάνοντας υπόψη και τις δύο λίστες κατάταξης. Στην παρούσα μελέτη, αξιοποιήθηκε κλινικό δείγμα 10 ζεύγών εικόνων τρισ-διάστατης σάρωσης, που αντιστοιχούσαν σε 2 χρονικά στιγμιότυπα, από σύστημα υπολογιστικής τομογραφίας πολλαπλών ανιχνευτών, Οι 10 ασθενείς έφεραν διάγνωση ILD με δευτερεύουσες παθήσεις συνδετικού ιστού που ακτινολογικά παρουσίαζαν απεικονιστικά πρότυπα θαμβής υάλου (ground glass opacities, ggo) και δικτυωτά (reticular) πρότυπα, συμπεριλαμβανομένου του προτύπου μελικηρύθρας (honeycomb). Οι σαρώσεις έγιναν σε χρονικές στιγμές που απείχαν μεταξύ τους περίπου 2 χρόνια. Η συνολική έκταση της ILD σε αυτές τις 10 σαρώσεις κυμαίνεται από 5% έως 80% (μέση τιμή: 32%). Η προσέγγιση χρήσης τεχνητών δεδομένων παρακολούθησης υιοθετείται ώστε η διαδικασία αξιολόγησης του πρώτου σταδίου να εστιάσει στο σφάλμα Αντιστοίχισης που προκαλείται από τα διάφορα υπό μελέτη σχήματα (αλγόριθμους) αντιστοίχισης και