FACULTY OF ENGINEERING TECHNOLOGY CAMPUS DE NAYER

Automated Segmentation of the Mandibular Nerve Canal in CBCT images

Arnout KONINGSVELD

Supervisor: Ing. K. Van Beeck Mentor: J. Keustermans

Master Thesis submitted to obtain the degree of Master of Science in Engineering Technology: Electronics-ICT, ICT

Academic Year 2013-2014

© Copyright KU Leuven Without written permission of the supervisor(s) and the author(s) it is forbidden to reproduce or adapt in any form or by any means any part of this publication. Requests for obtaining the right to reproduce or utilize parts of this publication should be addressed to KU Leuven campus De Nayer, Jan De Nayerlaan 5, B-2860 Sint-Katelijne-Waver, +32-15-316944 or via e-mail [email protected]. A written permission of the supervisor(s) is also required to use the methods, products, schematics and programs described in this work for industrial or commercial use, and for submitting this publication in scientific contests.

FACULTY OF ENGINEERING TECHNOLOGY CAMPUS DE NAYER

Automated Segmentation of the Mandibular Nerve Canal in CBCT images

Arnout KONINGSVELD

Supervisor: Ing. K. Van Beeck Mentor: J. Keustermans

Master Thesis submitted to obtain the degree of Master of Science in Engineering Technology: Electronics-ICT, ICT

Academic Year 2013-2014

© Copyright KU Leuven Without written permission of the supervisor(s) and the author(s) it is forbidden to reproduce or adapt in any form or by any means any part of this publication. Requests for obtaining the right to reproduce or utilize parts of this publication should be addressed to KU Leuven campus De Nayer, Jan De Nayerlaan 5, B-2860 Sint-Katelijne-Waver, +32-15-316944 or via e-mail [email protected]. A written permission of the supervisor(s) is also required to use the methods, products, schematics and programs described in this work for industrial or commercial use, and for submitting this publication in scientific contests.

Acknowledgements I would like to thank my supervisor Kristof Van Beeck and my mentor Johannes Keustermans for their guidance and suggestions during my thesis. I would also like to thank Medicim for allowing me to work with them and providing the data required to complete this thesis. My special thanks are extended to my lecturers at De Nayer for sharing their knowledge over the past few years. I wish to thank my parents for supporting me and allowing me to pursue my interests. And finally, I would like to thank my brothers and sister, my friends, my family and my girlfriend for their continued support. Thank you.

Arnout Koningsveld Sint-Katelijne-Waver, May 2014

ii

Abstract In this thesis we try to develop and implement a method enabling the segmentation of the mandibular nerve canal in CBCT images in an automated manner. This can help planning software for dental implant surgery respect safety margins around the mandibular nerve canal and thus reduce the risk of damaging it. This canal runs through the lower-jaw bone and contains nerves that run towards the lower lip and chin. Our method starts with the combination of a smoothing and gradient filter to reduce noise and enhance the edges of the canal, preparing the image for a fuzzy-connectedness method. After applying fuzzy-connectedness, we interpolate the results to fill in gaps and correct any errors. Using this method, we are able to segment large parts of the mandibular nerve canal with a computation time of about 12 seconds.

iii

iv

Abstract In deze thesis proberen we een methode te ontwikkelen en te implementeren die het mandibulair zenuwkanaal op een automatische manier segmenteert uit CBCT beelden. Dit kan planning software voor tandheelkundige implantaten helpen veiligheidsmarges in het oog te houden zodat er geen schade kan worden toegebracht aan het mandibulair zenuwkanaal. Dit kanaal loopt door het onderkaakbeen en bevat zenuwen die naar de onderlip en kin lopen. Onze methode begint met de combinatie van een smoothing en gradi¨ent filter om ruis te verminderen en de randen van het kanaal beter zichtbaar te maken. Op deze manier trachten we het volume voor te bereiden op een fuzzy-connectedness methode. Na het toepassen van fuzzy-connectedness interpoleren we de resultaten om gaten op te vullen en eventuele fouten te corrigeren. Met behulp van deze methode zijn we in staat om grote delen van het mandibulair zenuwkanaal te segmenteren met een rekentijd van ongeveer 12 seconden.

v

vi

Korte samenvatting Tegenwoordig wordt de plaatsing van tandheelkundige implantaten meestal gepland met behulp van drie-dimensionele medische beelden en diverse softwarepakketten. Door de lagere kost en lagere stralingsdosis in vergelijking met conventionele Computed Tomography (CT), verkiest men momenteel Cone-Beam CT (CBCT) als modaliteit om deze medische beelden te verwerven. Tijdens de planning moeten veiligheidsmarges rond bijvoorbeeld zenuwen worden gerespecteerd om deze niet te beschadigen. Het mandibulair zenuwkanaal loopt door het onderkaakbeen en bevat zenuwen die naar de onderlip en kin lopen. Figuur 1 geeft enkele mogelijke locaties van het mandibulair zenuwkanaal weer. Er is een aanzienlijke individuele variatie in de exacte vormgeving van dit kanaal in de mandibula (Johnson and Moore, 1997). Daarom kan de segmentatie van het mandibulair zenuwkanaal de planning software helpen rekening te houden met de veiligheidsmarges voor elk individueel geval. Het doel van deze thesis is dan ook een methode te ontwikkelen die op een automatische manier het mandibulair zenuwkanaal segmenteert.

Literatuur Tabel 1 geeft een overzicht van de bestaande methodes. De laatste jaren worden de methodes meer en meer geautomatiseerd en tevens wordt de nauwkeurigheid ervan beter. Het is ook duidelijk dat methodes toegepast op CT beelden meestal een hogere nauwkeurigheid hebben in vergelijking met CBCT beelden. Dit is te verwachten aangezien de kwaliteit van CT beelden hoger is dan CBCT beelden. De methode van Llor´ens et al. (2012) lijkt significant betere resultaten te bekomen op CT beelden. We merkten ook op dat de meeste moderne methodes beginnen met het lokaliseren van de foramen mentale en foramen mandibulare, gevolgd door het zoeken naar het pad dat deze verbindt. vii

Figure 1: Illustratie van het mandibulair zenuwkanaal. Deze figuur toont enkele mogelijke locaties van het mandibulair zenuwkanaal. Afbeelding genomen uit Johnson and Moore (1997).

viii

Table 1: Overzicht van bestaande methodes. De interactie kolom geeft weer of interactie met de gebruiker nodig is. E avg is de gemiddelde fout en E std de standaard afwijking. ”/” is gebruikt waneer informatie niet is gevonden of te veel afhangt van de nauwkeurigheid van de input van de gebruiker. Methode Stein et al.(1998) Lobregt et al.(2001) Hanssen et al.(2004) Kondo et al.(2004) Rueda et al.(2006) Sotthivirat et al.(2006) Yau et al.(2008) Kainmueller et al.(2009) Kim et al.(2011) Kroon (2011) Moris et al.(2012) Llor´ens et al.(2012)

Modaliteit CT CT CBCT CT CT CBCT CT CBCT CT CBCT CBCT CT

Interactie ja ja ja ja nee ja ja nee nee nee nee nee

E avg(mm) / / / 0.69 3.4 / / 1.00 0.73 1.88 0.99 0.14

E std(mm) / / / 0.76 1.11 / / 0.60 0.69 1.05 1.13 0.02

#datasets 5 / / 1 62 1 / 106 10 13 15 20

Methode Uit de literatuurstudie blijkt dat de methode van Llor´ens et al. (2012), die fuzzy-connectedness gebruikt, zeer goede resultaten bekomt op CT beelden. Wij hebben een methode ontwikkeld die ook fuzzy-connectedness gebruikt, en op CBCT beelden moet kunnen toegepast worden. We hebben deze methode ge¨ımplementeerd met behulp van ITK1 in C++. Figuur 2 geeft een overzicht van onze methode.

Figure 2: Overzicht van de methode. Een eerste belangrijke stap is de combinatie van een smoothing filter gevolgd door een gradi¨ent filter en thresholding. Deze combinatie wordt vaak gebruikt voor edge-detectie en wij trachten hiermee dan ook de randen van het zenuwkanaal beter zichtbaar te maken. Deze eerste stappen hebben als doel enkele minder gewenste effecten van CBCT te compenseren en het volume voor te bereiden op de fuzzy-connectedness methode. Vervolgens passen we twee keer een opeenvolging van fuzzy-connectedness en interpolatie toe. De eerste keer om het mandibulair zenuwkanaal in 1

http://www.itk.org

ix

de linkerkant van de mandibula te segmenteren en de tweede keer voor die in het rechtse deel van de mandibula. Dit werkt door een thresholding operatie uit te voeren op een fuzzy-connectedness scene (of fuzzy scene). De fuzzy scene is een beeld waarin de gegevens van elke pixel de fuzzy-connectedness van die pixel naar het startpunt (seed point) voorstelt. Deze fuzzy-connectedness wordt gedefinieerd als het sterkste pad tussen de pixel en het seed point. Hierbij is een pad een lijst van pixels die de pixel met het seed point verbindt. De sterkte van een pad wordt gedefinieerd als de zwakste fuzzy affinity tussen de pixels die het pad vormen. Fuzzy affinity is gedefinieerd als een functie die de kans dat twee aanliggende pixels tot hetzelfde object horen berekent. In de ITK implementatie wordt dit gedefinieerd als een Gauss functie van het verschil tussen de pixelwaardes en het verschil van de geschatte gemiddelde pixelwaarde van het object. Voor meer informatie over fuzzy-connectedness verwijzen we naar Udupa and Samarasekera (1996). In de implementatie van ITK kan men ook maar 1 seed point gebruiken. Wij hebben dit uitgebreid zodat we meerdere seed points kunnen toevoegen en zo betere resultaten kunnen bekomen. Na deze stap interpoleren we de resultaten van de fuzzy-connectedness filter. Deze interpolatie heeft als doel de resultaten te verbeteren en ze betrouwbaarder te maken. Tenslotte voegen we nog de resultaten van de twee kanten samen.

Resultaten We zijn van nul af aan begonnen en hebben een werkende methode ontwikkeld en ge¨ımplementeerd. Voor de ontwikkeling van de methode baseerden we ons vooral op de resultaten van de literatuurstudie en onze eigen kennis en ervaring. De methode kan op zowel CT als CBCT volumes toegepast worden. De methode is volledig ge¨ımplementeerd in C++ met behulp van ITK. We hebben dankzij een zorgvuldige selectie van filter-implementaties de rekentijd kunnen reduceren tot ongeveer 20 seconden. Dit wordt ongeveer 12 seconden als men de lees- en schrijfoperaties niet meerekent. We hebben de fuzzy-connectedness implementatie van ITK ook aangevuld om meerdere seed points te ondersteunen. Figuur 3 geeft aan dat elke extra seed de rekentijd verhoogt met iets meer dan 1 seconde. Helaas waren we niet in staat om automatisch seed points te genereren, dus de gebruiker moet zelf zoveel seed points toevoegen als hij/zij nodig acht. Hierdoor kan de actuele versie van onze methode hoogstens semiautomatisch genoemd worden.

x

Figure 3: Rekentijd in functie van het aantal seed points. Het effect dat het toevoegen van seeds heeft op de rekentijd wordt weergegeven. Het lezen en schrijven van het volume is inbegrepen in deze metingen. Over het algemeen kunnen we wel concluderen dat de methode met succes grote delen van het mandibulair zenuwkanaal kan segmenteren. Maar om de resultaten te verbeteren moet men af en toe nog wat wijzigen aan enkele parameters. Hierdoor is er uiteraard nog wel ruimte voor verbetering. Zo zou er nog wat kunnen veranderen aan de implementatie van de interpolatie, en zou de rekentijd ook nog verder gereduceerd kunnen worden. Maar het toevoegen van automatische seed point generatie zou volgens ons de methode het meest verbeteren. Dit zou de methode gebruiksvriendelijker maken en de resultaten verbeteren. Er zou minder interpolatie nodig zijn en de kwaliteit van de interpolatie verbeteren daar waar het nog nodig is.

xi

xii

Contents 1 Introduction 1.1 Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Goal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 2 2

2 Literature Review 2.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . 2.1.1 Anatomy . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Computed Tomography and Cone beam computed tomography . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Segmentation . . . . . . . . . . . . . . . . . . . . . . . 2.2 Literature Overview . . . . . . . . . . . . . . . . . . . . . . . 2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 3 3 6 7 9 13

3 Method 3.1 Framework and tools . . . . . . . 3.2 Method Overview . . . . . . . . . 3.3 Method Details . . . . . . . . . . 3.3.1 Input data and reading . 3.3.2 Smoothing . . . . . . . . 3.3.3 Gradient . . . . . . . . . . 3.3.4 Fuzzy-connectedness . . . 3.3.5 Interpolation . . . . . . . 3.3.6 Recombine . . . . . . . . 3.3.7 Output data and writing

15 15 16 17 17 18 19 21 24 27 27

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

4 Results 29 4.1 General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.2 Computation time . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5 Conclusion

35

xiii

xiv

Acronyms 2D two-dimensional. 7, 11 3D three-dimensional. 1, 7, 9, 10, 12, 17–21, 25, 27, 28 AAM Active Appearance Model. 12 ASM Active Shape Model. 12 CBCT Cone-Beam CT. vii, ix, x, 1, 3, 7, 9, 11, 13, 16, 17, 21, 22, 29, 31–33, 35 CT Computed Tomography. vii, ix, x, xv, 1, 3, 6, 7, 10–13, 16, 22, 29, 31, 35 DICOM Digital Imaging and Communications in Medicine. 15–18, 27, 28, 30, 35 FMM Fast Marching Method. 10, 12 MEDLINE Medical Literature Analysis and Retrieval System Online. 9 MRI Magnetic Resonance Imaging. 6 SNR Signal-to-Noise Ratio. 7, 13, 21 US Ultrasound. 6

xv

xvi

List of Figures 2.1 2.2 2.3 2.4

Illustration of the mandible . . . . . . . . . . . . Illustration of the mandibular nerve canal . . . . Illustration of the filtered backprojection method CT slice compared to CBCT slice . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

4 5 6 8

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10

Method Overview . . . . . . . Method: Reading input . . . Method: Smoothing image . . Smoothing: Effect of Sigma . Method: Gradient . . . . . . Gradient: Effect of Sigma . . Method: fuzzy-connectedness Method: Interpolation . . . . Method: Recombine . . . . . Method: Writing output . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

16 17 18 19 19 20 21 24 27 28

4.1 4.2 4.3

Computation time in function of the number of seed points . Results: CBCT volume . . . . . . . . . . . . . . . . . . . . . Results: CBCT volume (continued). . . . . . . . . . . . . . .

30 32 33

xvii

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

xviii

List of Tables 2.1

Overview of existing methods . . . . . . . . . . . . . . . . . .

14

4.1

Specifications of our test system

30

xix

. . . . . . . . . . . . . . . .

xx

Chapter 1

Introduction Nowadays dental implant surgery is usually planned to ensure its succes. With the help of three-dimensional (3D) medical images and various software systems, the clinician can make a virtual planning of the procedure. Due to its lower cost and lower radiation dose compared to conventional Computed Tomography (CT), Cone-Beam CT (CBCT) is currently the modality of choice for the acquisition of these medical images. During the planning of dental implants safety margins around, for example, nerves should be respected to reduce the risk of damaging them. The mandibular nerve canal runs through the lower-jaw bone and contains nerves that run towards the lower lip and chin. There is significant individual variation in the exact arrangement of this canal within the mandible (Johnson and Moore, 1997). This is why extraction of the mandibular nerve canal can help the planning software in taking the safety margins for each individual case into account.

1.1

Setting

This thesis is made as a master’s thesis for the Katholieke Universiteit Leuven (KU Leuven) in collaboration with Medicim1 . Medicim originally started as a KU Leuven spin-off in 2002. They develop innovative medical imaging solutions for treatment planning and guided surgery. In 2008 it was acquired by Nobel Biocare2 , which focusses on dental implants. One of their products is the software NobelClinician3 , which can be used by dentists for diagnosis, treatment planning and transfer of this planning to the patient and operating room. They also have a software product called Maxilim which can be used for orthognathic surgery planning. 1

http://www.medicim.com http://www.nobelbiocare.com 3 http://www.nobelclinician.com/ 2

1

2

1.2

CHAPTER 1. INTRODUCTION

Goal

The main goal of this thesis is to develop an algorithm enabling the visualisation of the mandibular nerve canal in an automated manner while maintaining an acceptable level of accuracy. In order to achieve this a suitable method for its segmentation must be implemented. We will search and compare existing methods in literature to validate their suitability. For this evaluation we will take their accuracy and level of automation into account. Afterwards we will develop and implement our own method.

1.3

Outline

In chapter 2 we review the existing methods found in literature after briefly explaining some theoretical background on the subject. Next we present our own method in chapter 3, starting with a brief explanation of the tools we have used to implement it in section 3.1, followed by an overview in section 3.2 and ending with a more detailed explanation in section 3.3. Chapter 4 displays and discusses our results. Finally, we give our conclusions and recommendations for future work in chapter 5.

Chapter 2

Literature Review In this chapter we will review various methods for segmenting the mandibular nerve canal found in literature. We will start with an explanation on the theoretical background needed to understand the subject-matter of these methods. Next, we give an overview of the proposed methods and in the last section of this chapter we will compare and discuss them.

2.1

Theoretical Background

In this section we will introduce and explain several key concepts to support the contents of the literature review and this thesis in general. First we briefly dive into the anatomy of the lower jaw to familiarize ourselves with the location and surroundings of the mandibular nerve canal. Then we will give some information about CT and CBCT, which is the technique used to obtain the anatomical data that make up the medical images. Finally the concept of image segmentation will be explained and an overview of some of its basic methods will be given.

2.1.1

Anatomy

The lower jaw, or mandible, illustrated in figure 2.1 consists of 3 main parts. A horseshoe-shaped body and two rami on each end of the body. The junction between the body and ramus nearly forms a right angle and is aptly named angle. Initially two halves of the mandible are united by a joint called the symphysis menti, but during the first year of postnatal life the halves fuse together obliterating this joint. The line of fusion creates a vertical ridge below the space between the central incisor teeth. This ridge runs into a raised area called the mental protuberance. From a lateral point of view this raised area may be more prominent on each side of the mental protuberance. This is called the mental tubercle. A small ridge 3

4

CHAPTER 2. LITERATURE REVIEW

Figure 2.1: Illustration of the mandible. Image taken from http://oralmaxillo-facialsurgery.blogspot.be/2010/05/facial-boneanatomy.html.

runs from the mental tubercle backwards and upwards towards the ramus. This ridge becomes more prominent below the last molar tooth and ends up running continuous with the edge of the ramus. This ridge is called the oblique line. The mental foramen lies just above this oblique line in the region of the premolar teeth. It is a hole through which the mental nerve, which branches of the inferior alveolar nerve, and blood vessels pass. From here the mandibular canal, which contains this nerve, runs through the body to the mandibular foramen in the centre of the ramus. Figure 2.2 shows some possible locations of the mandibular nerve canal.

2.1. THEORETICAL BACKGROUND

5

Figure 2.2: This image shows some possible locations of the mandibular nerve canal. Image taken from Johnson and Moore (1997).

6

CHAPTER 2. LITERATURE REVIEW

Figure 2.3: Illustration of the filtered backprojection method applied on the Shepp Logan phantom. From left to right and top to bottom: Original phantom, Resulting sinogram, Sinogram filtered with a rampfilter, Reconstructed phantom.

2.1.2

Computed Tomography and Cone beam computed tomography

Computed Tomography (CT) is a technique often used for medical imaging purposes. Although other imaging techniques like Ultrasound (US) and Magnetic Resonance Imaging (MRI) exist, CT is most commonly used for this application(Suetens, 2009). It is also most relevant for this thesis. In 1972 Godfrey N. Hounsfield developed the first CT scanner. His work was based on methods developed by A. M. Cormack. In 1979 they received the Nobel Prize in Physiology or Medicine for their work. CT produces cross-sectional images or slices that represent the x-ray attenuation properties of the scanned object. To produce an image an x-ray source transmits beams through the body which attenuates them. On the opposite side x-ray detectors measure the attenuated beams. The source and detectors rotate around the body and the process is repeated for a large number of angles. The output of these measurements is the Radon Transform of the body. A sinogram is a visual representation of this raw data. Each point in the body will give rise to a sinusoid in the sinogram. These measurements can then be used to reconstruct an image using for example the filtered back projection method (Gonzalez and Woods, 2002), illustrated in figure 2.3.

2.1. THEORETICAL BACKGROUND

7

The Hounsfield values, H, or sometimes called CT-values represent the x-ray attenuation and are displayed in the resulting image. They are defined as: µ − µH 2 O H= µH2 O where µ is the linear attenuation coefficient. CT-values are expressed in Hounsfield Units (HU). From this definition we can see that water has a CT-value of 0 HU. Air has a CT-value of around -1000 HU and the CT-value of bone varies from 250 HU up to 1000 HU and higher. Unfortunately we will soon see that the Hounsfield scale is not usable for CBCT because of distortion (De Vos et al., 2009). CT is commonly used to look at bony structures. A scan is fast (less than 5 minutes) and gives a relatively high radiation dose (several mSv) (Suetens, 2009). Cone-Beam CT (CBCT) is a CT image acquisition technique where the x-ray beam is cone-shaped and a flat panel detector is used. One rotation around the body produces a series of two-dimensional (2D) images which can be combined to produce a 3D volume. CBCT produces images with higher resolution using a lower radiation dose and is cheaper compared to conventional CT. However, a CBCT scan typically takes longer and thus is more prone to motion artifacts. The size of the detector is also limited resulting in a limited field of view. Moreover, due to photon scattering and the use of multiple rows of detectors, the Signal-to-Noise Ratio (SNR) is lower. The lower radiation dose further enhances this effect. Finally, CBCT images suffer from inhomogeneity of the gray values. Thus, as mentioned before, the Hounsfield scale is not usable. Figure 2.4 shows a CT slice (2.4(a)) next to a CBCT slice (2.4(b)). We can clearly see that the CBCT slice has a much lower SNR.

2.1.3

Segmentation

In this context segmentation is the partitioning of an image or volume into a number of connected regions. Every pixel or voxel should belong to one (and only one) of these regions. In medical applications it is typically used to segment for example organs, different types of tissue, tumors, or other anatomical structures. Images or volumes are generally segmented based on similarity/difference, uniformity or homogenity of the pixels or voxels. These criteria can be evaluated in position, intensity value, color, graininess, texture, orientation, variance or other features. The most basic methods of segmentation are pixel (or voxel) based, which means that they only take the features of each individual pixel into account. Thresholding is an example of such a simple and yet frequently used method. A certain value (e.g. gray value) is selected and every pixel with a value

8

CHAPTER 2. LITERATURE REVIEW

(a) Slice of CT volume.

(b) Slice of CBCT volume.

Figure 2.4: CT slice compared to CBCT slice. The CBCT slice clearly has a worse SNR. higher than the chosen value will be part of one class, while the pixels with a value lower than the chosen value will be part of another class. The drawback of this simple method is the difficulty to introduce spatial knowledge. One way to improve this is using adaptive threshold by changing the threshold proportional to the local average, which means also taking the features of other pixels in the region into account. This leads us to the region based segmentation methods. An obvious example for this type is the region growing method. This method starts from a seed (the starting point) and includes neighbours if they are sufficiently similar to the already included neighbours. A problem here is that it can be difficult to decide when to include and when to exclude a neighbour. Leakage occurs when the pixels of the edge between two regions aren’t considered sufficiently different resulting in one big region instead of two smaller ones. On the other hand oversegmentation occurs when pixels are too often considered sufficiently different resulting in many small regions instead of a few larger ones. A third type of segmentation methods are so called edge based and can for example be used in combination with region growing to prevent leakage by detecting the edges. A simple example of this type is the live wire method (Mortensen and Barrett, 1995) which computes the cheapest path between two points with the cost based on certain criteria like the ones mentioned earlier. Many other more complex methods are hybrids of these types of methods. After these methods some post-processing steps are usually applied. Sim-

2.2. LITERATURE OVERVIEW

9

ple examples of this are splitting regions into smaller ones or merging regions into larger ones.

2.2

Literature Overview

Now that we have a basic understanding of the theoretical background we can start explaining the existing methods. The purpose of this section is to give an overview of the methods found in literature. We will then compare and discuss these methods in the next section. We will briefly explain where and how we found the literature in a first paragraph. Afterwards the methods will be presented in a chronological order which should give you an idea of how these methods evolved. We have searched for literature regarding methods for segmentation of the mandibular nerve canal in CBCT data in several libraries and databases. First of all we used Limo1 , the search platform for the collections and libraries of KU Leuven. We also made use of our access to LiU’s UniSearch2 to scan through their collection as well. A third source for literature was PubMed3 , which is a well-known free database comprising more than 23 million citations for biomedical literature from Medical Literature Analysis and Retrieval System Online (MEDLINE), life science journals, and online books. Google Scholar4 , a web search engine that indexes academic literature, was also used. Finally, we followed the cross-references in the resulting literature to expand our search even further. The search resulted in the following articles (ordered chronologically): 1998 ”Tracing of Thin Tubular Structures in Computer Tomographic Data” by Stein et al. (1998). In this paper the authors present a detection algorithm that works as follows: a user selects a start- and endpoint. Afterwards they treat the voxel set as a graph and apply a modified Dijkstra’s algorithm to construct the path between the start- and endpoint. To reduce computation time they apply a series of erosions and dilations in 3D in order to reduce the search space. The resulting path is then used as a guidewire to initialize a balloon snake which will inflate and fill the canal.

1

http://bib.kuleuven.be/ http://www.bibl.liu.se/?l=en 3 http://www.ncbi.nlm.nih.gov/pubmed/ 4 http://scholar.google.com/ 2

10

CHAPTER 2. LITERATURE REVIEW

2001 ”Dental implant surgery: planning and guidance” by Lobregt et al. (2001). Here the authors have chosen to model the canal with a 3D polybezier path. To construct it a user should select a few points in several views. The constructed path is then displayed in all views. 2004 a) ”Nerves - Level sets for interactive 3D segmentation of nerve channels” by Hanssen et al. (2004). The method described in this paper is similar to the method presented by Stein et al. (1998). It also requires a user to provide a start- and endpoint. Afterwards a Fast Marching Method (FMM) is used to compute a path between these points that passes through the nerve canal. This path is used to initialize an active surface which is further computed using geodesics and level sets. b) ”Computer-based extraction of the inferior alveolar nerve canal in 3-D space” by Kondo et al. (2004). In this paper the described method works by computing panoramic images from CT data, followed by detection of hollow canals and finally a 3D line tracking algorithm. The computation of panoramic CT images works as follows: a slice of CT image that cointains the mandible is selected. Then thresholding is applied to extract objects. Holes in the objects are filled and if multiple objects remain a user should select the target object (the mandible). Afterwards the midline of the mandible is determined and curve fitting to this midline is applied to set up some sample points to be used for interpolation of intensity values. This interpolation scheme is repeated for other CT images that contain the mandible. To detect hollow canals the panoramic CT images are first enhanced by a gray-level transformation. Then gradient orientation is used for the detection of hollow canals. Finally, the 3D line tracking algorithm uses a guide mask and is performed bi-directionally between two terminal points provided by a user. The guide mask is used to constrain the direction of the path and thus steering it towards the terminal points and performing it bi-directionally should reduce the likelihood of errors introduced by spurious lines.

2.2. LITERATURE OVERVIEW

11

2006 a) ”Automatic Segmentation of Jaw Tissues in CT Using Active Appearance Models and Semi-automatic Landmarking” by Rueda et al. (2006) . For this method an Active Appearance Model approach is used. The training set used to train the model requires some manual input of landmarks in 2D images. Based on these landmarks and contours found by thresholding a large number of pseudolandmarks are defined. Once the model is trained it can be used to segment new images. b) ”Automatic Detection of Inferior Alveolar Nerve Canals on CT Images” by Sotthivirat and Narkbuakaew (2006). The method of this paper is applied on panoramic CT images. These are computed with the help of user-defined points on several slices. Then they apply linear contrast stretching to enhance the contrast and a Gaussian lowpass filter for noise reduction. After this pre-processing Canny edge detection is applied, followed by the distance transform of the result which is then converted to a binary image using thresholding. Next, morphological opening and closing operations are applied to remove unwanted connections. Finally they use the characteristics and properties of the nerve canal to extract it from the remaining groups of pixels. 2008 ”An adaptive Region Growing Method to Segment Inferior Alveolar Nerve Canal from 3D Medical Images for Dental Implant Surgery” by Yau et al. (2008). This method also starts by computing some sort of panoramic volume with the help of user-defined points. In a next step they use an adaptive region growing method on cross sectional images from this volume. A user should provide one seed point and select a terminal image. The user-provided seed point is used to initialize the iterative generation of new seed points. After the terminal image has been processed the method will come to an end. 2009 ”Automatic Extraction of Mandibular Nerve and Bone from ConeBeam CT Data” by Kainmueller et al. (2009). This paper presents one of the first methods that is applied on CBCT data. They used datasets labeled by experienced dentists to generate a statistical shape model. This statistical shape model is used to perform some sort of registration by adapting it to fit the targeted CBCT data. The result can already yield an approximation of the

12

CHAPTER 2. LITERATURE REVIEW nerve but it can also be used as an initialization for a next step, which applies a modified Dijkstra algorithm to optimize the result.

2011 a) ”Automatic Extraction of Inferior Alveolar Nerve Canal Using Feature-Enhancing Panoramic Volume Rendering” by Kim et al. (2011). The method described in this paper works in three main steps. First they automatically segment the mandible from the CT images using high-value thresholding and connected component labeling followed by region growing and morphological closing. In the next step they search for the mental and mandibular foramina which will serve as the start- and endpoints for a line-tracking algorithm in the final step. They generate two 3D panoramic volume rendering images, one for the detection of the mental foramen in a front view and one for the mandibular foramen in a back view, using a ray-casting algorithm. Texture analysis is performed on these images to automatically detect the mental and mandibular foramina. Finally a modified Dijkstra’s algorithm is used to track the canal path from mental to mandibular foramen. The full canal is then extracted using a Fast Marching Method (FMM) on the resulting path. b) ”Segmentation of the Mandibular Canal in Cone-Beam CT Data” by Kroon (2011). Kroon discusses several methods in this thesis and tries to improve an Active Shape Model (ASM) method as used by Kainmueller et al. (2009) and an Active Appearance Model (AAM) method as used by Rueda et al. (2006). The ASM method is altered by adding several post-processing methods like a filter to remove noise, using 3D snake optimization and using fast-marching. These modifications however did not succeed in improving the accuracy. The improvement of the AAM method was a succes though. By using weighted intensities the mean error decreased from 2.0mm to 1.88mm. 2012 a) ”Automated tracking of the mandibular canal in CBCT images using matching and multiple hypotheses methods” by Moris et al. (2012). This method also starts with detection of the mental and mandibular foramina by using the recurring discontinuities in the border

2.3. DISCUSSION

13

of the mandible. Next, they track the path bi-directionally and combine the results to produce a single path. b) ”Jaw tissues segmentation in dental 3D CT images using fuzzyconnectedness and morphological processing” by Llor´ens et al. (2012). The approach used in this method is similar to the one by Kim et al. (2011). They split the jaw into five regions and for each slice try to determine which region it belongs to. The regions are defined as follows: everything before the mandibular foramen, the mandibular foramen, everything in between the two foramina, the mental foramen and everything after the mental foramen. This way they know the nerve enters the canal in the second region, runs through the third and exits it in the fourth region. In a next step they use a fuzzy connectedness method. To reduce the chances of failure for this method they use pseudoorthopantomographic projections to determine where it is hard to detect the nerve and to estimate a curve that fits the nerve. The curve provides the seed points for the fuzzy-connectedness method. For the cross-sections where it is hard to detect the nerve they use interpolation instead of the fuzzy connectedness method to fill in the gaps of the segmentation. Finally a marching cubes algorithm is used to reconstruct the 3D volume.

2.3

Discussion

Now that we know the basic functionality of these existing methods we will discuss and compare them in this section. The purpose of this section is to argument our selection of suitable methods from all of the methods mentioned and give another overview of the existing methods. To facilitate the comparison we made Table 2.1 which gives an overview of the reviewed methods with respect to the criteria previously mentioned in the goal of this thesis. Over the years the methods have evolved from manual to semi-automated and automated operation. Moreover, the accuracy of the automated methods has improved significantly. We can see that the modality has a clear impact. Methods applied on CT images generally acquire more accurate results. This is expected as CBCT images tend to have a lower SNR (see section 2.1.2 on page 6). The method proposed by Llor´ens et al. (2012) also seems to perform significantly better than others, albeit on CT images. Most modern methods start by locating the mental and mandibular foramina, followed by finding the path of the nerve in between those points. It seems clear that we should also start this way.

14

CHAPTER 2. LITERATURE REVIEW

Table 2.1: Overview of existing methods. The interaction column displays if user interaction is required. E avg represents the average error and E std the standard deviation. ”/” is used when this information was not found or too dependent on the accuracy of the users input. Method Stein et al.(1998) Lobregt et al.(2001) Hanssen et al.(2004) Kondo et al.(2004) Rueda et al.(2006) Sotthivirat et al.(2006) Yau et al.(2008) Kainmueller et al.(2009) Kim et al.(2011) Kroon (2011) Moris et al.(2012) Llor´ens et al.(2012)

Modality CT CT CBCT CT CT CBCT CT CBCT CT CBCT CBCT CT

Interaction yes yes yes yes no yes yes no no no no no

E avg (mm) / / / 0.69 3.4 / / 1.00 0.73 1.88 0.99 0.14

E std (mm) / / / 0.76 1.11 / / 0.60 0.69 1.05 1.13 0.02

#datasets 5 / / 1 62 1 / 106 10 13 15 20

Chapter 3

Method In this chapter we will present our method. It comprises a brief explanation of the tools used for the implementation, an overview of the method itself, and finally, a section in which we examine each part of the method in more detail.

3.1

Framework and tools

We used a number of tools along the way to implement our method. This section will describe some of those which may be less known and explains why we used them. CMake CMake1 is a cross-platform, open-source build environment. It is written in C++ and can generate native build files for many platforms and IDEs. It was developed because ITK, which will be explained soon, needed a cross-platform build environment. DICOM Digital Imaging and Communications in Medicine (DICOM) is a standard which defines how medical images should handle, store, print and transmit their information. It defines a file format and a network communications protocol that is built on TCP/IP. The file format consists of a header and the image data. The header contains meta-information about the image like parameters and details about the scanner, modality, information about the patient, and many more. This ensures that the image is always accompanied by its meta-information.

1

http://www.cmake.org

15

16

CHAPTER 3. METHOD

ITK ITK2 stands for Insight Segmentation and Registration Toolkit. It is an open-source application development framework for performing registration and segmentation on multidimensional data. ITK uses CMake to manage its build process, thus making it cross-platform. It is written in C++ and makes frequent use of the C++ template feature and generic programming. MeVisLab MeVisLab3 is a cross-platform application framework. It is designed for medical image processing and visualisation, and features an IDE for graphical programming. It is written in C++ and Python in combination with Qt4 for the graphical user interfaces. We used ITK to implement every filter in our method while MeVisLab was used to visualise our results. Our input datasets are stored in DICOM format and we also write our results back in DICOM format.

3.2

Method Overview

After our literature review it was clear that the method proposed by Llor´ens et al. (2012), which uses fuzzy-connectedness, performed better than others. However, their method is applied on CT images. We have developed our own method, also using fuzzy-connectedness, and apply it on CBCT images. An overview of our method is given in figure 3.1.

Figure 3.1: Overview of the method. After reading an image we perform a smoothing operation to reduce noise. Then we calculate the gradients of the smoothed image. After applying a threshold the pipeline splits in two, one for the left nerve and one for the right nerve. We use a fuzzy-connectedness method followed by interpolation to fill in gaps and correct errors on both sides. Finally we recombine the results into one image and write this result as output. 2

http://www.itk.org http://www.mevislab.de 4 http://qt-project.org/ 3

3.3. METHOD DETAILS

17

Thresholding an image created by a smoothing filter followed by a gradient filter is a common way to detect edges. This way we enhance the edges of the nerve canal and try to reduce the negative aspects of CBCT images before applying fuzzy-connectedness. For the interpolation we use a B-spline interpolation filter to fit our data. All input and output volumes are written in DICOM format.

3.3

Method Details

In this section we will take a closer look at the method by explaining the input, output, parameters, inner workings and choices we have made for each individual step in the method. Each step will be accompanied by a figure depicting the input, output and parameters. Horizontal lines represent input and output, while vertical lines represent parameters.

3.3.1

Input data and reading

Figure 3.2: Reading input data. The location of the DICOM directory is a parameter, the output is the 3D image. This step loads the image into memory. Figure 3.2 shows the part of the block diagram that is relevant for this step. Input DICOM files Output 3D image Parameters Input DICOM directory: The location of the directory that contains the input image in DICOM-format. Implementation • itk::GDCMImageIO • itk::GDCMSeriesFileNames

18

CHAPTER 3. METHOD • itk::ImageSeriesReader

The first series of DICOM files in the given directory are read and the 3D image is loaded into memory. Some meta-data is also collected, this is to be used when we write our final result in DICOM-format in the final step.

3.3.2

Smoothing

Figure 3.3: Smoothing image. The parameter Sigma (σ) determines how much smoothing is applied on the input image. In this step a smoothing operation is applied to the image in order to reduce noise. Figure 3.3 shows the part of the block diagram that is relevant for this step. Input 3D image Output Smoothed 3D image Parameters Sigma (σ): This determines how much smoothing is applied on the input image. Implementation itk::SmoothingRecursiveGaussianImageFilter The image is smoothed by approximate convolution with Gaussian kernels. The effects of the parameter sigma are illustrated on a single slice in figure 3.4. It is important not to apply too much smoothing to the image, because then the details needed to segment the nerve might get lost (see figure 3.4(d)). On the other hand, if not enough smoothing is applied to the image, the gradient filter in the next step might pick up gradients resulting from noise, propagating problems to the next steps. These problems will be illustrated in the next section.

3.3. METHOD DETAILS

19

(a) Original image

(b) σ = 0.2

(c) σ = 0.5

(d) σ = 0.7

Figure 3.4: Effect of Sigma when smoothing the image. Higher σ means more smoothing. If σ is too high the image loses too many details.

3.3.3

Gradient

Figure 3.5: Calculate the gradient of the image. The gradient of the smoothed image is calculated. Figure 3.5 shows the part of the block diagram that is relevant for this step. Input Smoothed 3D image

20

CHAPTER 3. METHOD

Output Thresholded gradient 3D image Parameters Threshold: The value for the thresholding operation. Implementation • itk::GradientMagnitudeImageFilter • itk::ImageDuplicator • itk::ThresholdImageFilter The combination of a smoothing and gradient filter followed by thresholding is frequently used for edge detection. Our goal here is to make the edges of the canal more visible. The effects of the smoothing step on our gradient image are illustrated in figure 3.6.

(a) Gradient of original image

(b) Gradient of image after smoothing operation with σ = 0.2

(c) Gradient of image after smoothing oper- (d) Gradient of image after smoothing opation with σ = 0.5 eration with σ = 0.7

Figure 3.6: Gradient: Effect of Sigma.

3.3. METHOD DETAILS

21

In figure 3.6(a) we can see that the filter picks up a lot of gradients in the noise. After performing a smoothing operation on the image, this effect is reduced but still visible in figure 3.6(b). In figure 3.6(c) this effect is almost invisible, but the image has started to lose some of its details. Even more details are lost and even the edges start to become less pronounced when we increase σ further. This can be seen in figure 3.6(d). The main reason why we use the gradient filter is to make it possible to use our method on CBCT images. As previously mentioned in section 2.1.2, a big problem that we encounter with CBCT images is that they suffer from inhomogeneity of the gray values and a low SNR. The gradient only takes the difference in intensity between pixels/voxels into account. It doesn’t evaluate absolute intensity levels, it evaluates them relative to each other. This means that if an offset is added to a volume, it will still produce the same results. And even if an offset is added only locally, it will still produce the same results there. After this step we apply a quick thresholding operation to filter out any remaining small gradients picked up in the noise. The threshold value shouldn’t be too high in order to keep the edges of the nerve canal in the image and only filter out the gradients caused by remaining noise. Afterwards we split the pipeline in two. One branch to segment the left nerve canal, and one branch for the right nerve canal. The gradient image is duplicated and each instance passes through the next set of filters. This is done because the fuzzy-connectedness filter can only segment one object at a time. We create two instances of the same set of filters and feed them different seed points.

3.3.4

Fuzzy-connectedness

Figure 3.7: Calculate the fuzzy-connectedness. In this step we apply a fuzzy-connectedness method to extract the nerve. Figure 3.7 shows the part of the block diagram that is relevant for this step. Input Thresholded gradient 3D image Output Binary 3D image

22

CHAPTER 3. METHOD

Parameters • Mean: The estimated object mean • MeanDif: The estimated mean difference between neighboring pixels for the object • Seed: Seed point, believed to be inside the object • Threshold: Threshold value to perform the thresholding operation on the calculated fuzzy scene • Variance: The estimated object variance • VarianceDif: The estimated variance of the difference between pixels for the object • Weight: The weight of the first term in the affinity computation Implementation • itk::MaximumImageFilter • itk::SimpleFuzzyConnectednessScalarImageFilter In our literature review we noticed that a method using fuzzy-connectedness performed best on CT images. We decided to apply this on CBCT images. The previous two steps were used to prepare our CBCT image for segmentation using fuzzy-connectedness. To perform the segmentation thresholding is applied on a fuzzy-connectedness scene, or fuzzy scene. The fuzzy scene is an image where the data of each pixel represents the fuzzy-connectedness of that pixel to the seed point. The fuzzy-connectedness of a pixel is defined as the strongest path between the pixel and the seed point. A path is a list of pixels that connect a pixel with the seed point. The strength of a path is defined as the weakest fuzzy affinity between the pixels that form the path. Fuzzy affinity is defined as a function that calculates the probability that two neighboring pixels belong to the same object. In ITK’s implementation this is defined as a Gaussian function of the pixel difference and the difference of the estimated object mean. For more information about fuzzy-connectedness, we refer to Udupa and Samarasekera (1996). The Threshold parameter determines the value for the thresholding operation on the calculated fuzzy-connectedness scene. A high value means that the fuzzy-connectedness must be high, or in other words the strongest path between a pixel and the seed point must be strong enough. So increasing the threshold usually decreases the size of the segmented object. A lower value will result in pixels with a weaker path, thus less certain to be part of the object, to be added to the result. The seed point is arguably the most important parameter of all, because in ITK’s implementation only one seed point can be specified. This means

3.3. METHOD DETAILS

23

the whole fuzzy scene is calculated using the fuzzy-connectedness of each pixel to that one seed point. Because of this, finding one suitable seed proved to be extremely difficult. To solve this problem we added multiple seed support by calculating multiple fuzzy scenes for different seed points and adding them together using itk::MaximumImageFilter. Listing 3.1 shows how we implemented this for the left side. We should note that by adding more seeds, the computation time increases because a larger number of fuzzy scenes have to be calculated and combined. Each seed point added will result in an extra iteration of the loop in listing 3.1. Listing 3.1: Code snippet to combine the fuzzy scenes from the left fuzzyconnectedness filter. for(int i=0; iGetNumberOfPoints(); i++){ ImageType::IndexType seed; leftSeedPointSet->GetPointData(i, &seed ); FuzzySegmentationFilterType::Pointer leftFuzzysegmenter = FuzzySegmentationFilterType::New(); leftFuzzysegmenter->SetInput( leftImage ); leftFuzzysegmenter->SetObjectSeed( seed ); leftFuzzysegmenter->SetThreshold( threshold ); leftFuzzysegmenter->SetParameters(mean,variance,meandif,vardif,weight); leftFuzzysegmenter->Update(); if(i==0){ maxImageFilter->SetInput(0,leftFuzzysegmenter->GetFuzzyScene()); } else{ maxImageFilter->SetInput(0,maxImageFilter->GetOutput()); } maxImageFilter->SetInput(1,leftFuzzysegmenter->GetFuzzyScene()); maxImageFilter->Update(); }

The result of the code is a fuzzy scene that should have more high certainty values than any other fuzzy scene calculated using just one seed. We could therefore increase the threshold parameter and still get a good, or possibly even better result. Normally ITK’s implementation of the fuzzyconnectedness filter returns a binary image that is the result of thresholding the input according to the data contained in the fuzzy scene. But since we combine multiple results by only extracting the fuzzy scenes, we still have to manually threshold the volume according to the final fuzzy scene. Listing 3.2 shows how we implemented this.

24

CHAPTER 3. METHOD

Listing 3.2: Code snippet to threshold the left image according to the results of the final fuzzy scene. typedef itk::ImageRegionIteratorWithIndex FuzzyIteratorType; FuzzyIteratorType fuzzyLeftIt( maxImageFilter->GetOutput(), maxImageFilter->GetOutput()->GetLargestPossibleRegion() ); typedef itk::ImageRegionIteratorWithIndex ImageIteratorType; ImageIteratorType fuzzyLeftOt( leftImage, leftImage->GetLargestPossibleRegion()); while( !fuzzyLeftIt.IsAtEnd()) { if( fuzzyLeftIt.Get() > threshold ) { fuzzyLeftOt.Set( itk::NumericTraits::max() ); } else { fuzzyLeftOt.Set( itk::NumericTraits::min() ); } ++fuzzyLeftIt; ++fuzzyLeftOt; }

Naturally, the same operations are also applied on the right side. All other parameters are used in the Gaussian function for the affinity computation and should reflect the statistic properties of the desired object.

3.3.5

Interpolation

Figure 3.8: Interpolate results. The results of the fuzzy-connectedness filter are interpolated. Figure 3.8 shows the part of the block diagram that is relevant for this step.

3.3. METHOD DETAILS

25

Input Binary 3D image Output Binary 3D image Parameters • Spline order: order of the spline • Number of control points: defines the number of control points in each dimension at the coarsest level. • Number of levels: defines the number of levels used in the multilevel approach • Close dimension: defines if the spline should be periodic in certain dimensions • Size: the size of the parametric domain • Spacing: the spacing of the parametric domain • Origin: the origin of the parametric domain Implementation itk::BSplineScatteredDataPointSetToImageFilter We perform an interpolation of the results to make them more reliable. For each slice in the Anterior to Posterior direction, the mean of the xyz coordinates of each point (that is thought to be part of the nerve) is calculated and the resulting point is added to a pointset. This is a data structure to represent a set of points in n-dimensional space. Listing 3.3 shows our implementation of this part for the left side. Listing 3.3: Code snippet to calculate the mean of the xyz coordinates of each point (that is thought to be part of the nerve) and add the resulting point to a pointset. const unsigned int ParametricDimension = 1; const unsigned int DataDimension = 3; typedef float RealType; typedef itk ::Vector VectorType; typedef itk ::Image VectorImageType; typedef itk::PointSet PointSetType ; PointSetType::Pointer leftPointSet = PointSetType::New (); leftImage->SetPixel(leftMental, itk::NumericTraits::max());

26

CHAPTER 3. METHOD

leftImage->SetPixel(leftMandibular, itk::NumericTraits::max()); typedef itk::ImageSliceConstIteratorWithIndex< ImageType > SliceIteratorType; SliceIteratorType leftIt( leftImage, leftImage->GetRequestedRegion() ); leftIt.SetFirstDirection(0); leftIt.SetSecondDirection(2); double avgx=0; double avgy=0; double avgz=0; double count=0; leftIt.GoToBegin(); while( !leftIt.IsAtEnd() ) { while( !leftIt.IsAtEndOfSlice() ) { while( !leftIt.IsAtEndOfLine() ) { if( leftIt.Get() == itk::NumericTraits::max() ) { avgx=(leftIt.GetIndex()[0]+avgx*count)/(count+1); avgy=(leftIt.GetIndex()[1]+avgy*count)/(count+1); avgz=(leftIt.GetIndex()[2]+avgz*count)/(count+1); count++; } ++leftIt; } leftIt.NextLine(); } if(avgx!=0){ PointSetType::PointType point; unsigned long i = leftPointSet -> GetNumberOfPoints (); point[0]=i; leftPointSet ->SetPoint( i, point ); PointSetType::PixelType V( DataDimension ); V[0] = static_cast( avgx); V[1] = static_cast( avgy); V[2] = static_cast( avgz); leftPointSet -> SetPointData ( i, V ); } avgx=0; avgy=0; avgz=0; count=0; leftIt.NextSlice(); }

3.3. METHOD DETAILS

27

The idea here is to eliminate possible outlying points and only add the ”centre” of the nerve to the pointset. We feed the resulting pointset to the interpolation filter to fit a spline to these points. Afterwards we sample the result of the interpolation filter and add the samples to our image. This will fill in possible gaps in our data and should result in a smooth nerve canal. For more information about the specific interpolation method used inside the filter we refer to Tustison and Gee (2006).

3.3.6

Recombine

Figure 3.9: Recombine left and right image. The split pipeline is joined again by recombining the left and right image. Figure 3.9 shows the part of the block diagram that is relevant for this step. Input • Left binary 3D image • Right binary 3D image Output Binary 3D image Parameters none Implementation ikt::AddImageFilterType We merge the pipeline by simply adding the two binary 3D images together. This is implemented as a simple pixel-wise summation of the two images, as the left and right nerve shouldn’t overlap.

3.3.7

Output data and writing

Finally we write our final image in a new DICOM directory. Figure 3.10 shows the part of the block diagram that is relevant for this step. Input Binary 3D image Output DICOM files

28

CHAPTER 3. METHOD

Figure 3.10: Writing output data. Parameters Output DICOM directory: The location of the directory that should contain the output image in DICOM format. Implementation • itk::GDCMImageIO • itk::GDCMSeriesFileNames • itk::ImageSeriesWriter The final 3D image is written in DICOM format in the given output directory. The meta-data collected in the first step is used for the details of the DICOM format. The result is a dataset with the same meta-data, but different image data.

Chapter 4

Results In this chapter we will review our results. First our own accomplishments will be presented. Next, our testing environment will be described briefly along with the computation time. And finally, we will discuss the results of our method.

4.1

General

We have developed and implemented a functioning method for the segmentation of the mandibular nerve canal from scratch. There was no previous work available to us to start from. For the development of the method we relied heavily on the results of the literature review and our own knowledge and experience. Our resulting method can be applied on both CT and CBCT images. This will be illustrated and explained in section 4.3 The implementation is written entirely in C++ using ITK. We have spent a significant amount of time reducing the computation time from 30+ minutes to within one minute by careful selection of filter implementations and other implementation details. We also added multiple seed support to complement ITK’s fuzzy-connectedness implementation. However, we were unable to automatically generate suitable seed points, so a user should still supply as many seed points as he/she considers necessary. Because of this, the current version of our method can at best be called semi-automated.

4.2

Computation time

The specifications of our test system are displayed in table 4.1. On this system our method, when applied on a CBCT volume, has a computation time of about 20 seconds when we include reading and writing the volume 29

30

CHAPTER 4. RESULTS Table 4.1: Specifications of our test system. Operating system IDE ITK version CMake version MeVisLab version Processor Memory

Windows 8.1 Pro 64-bit Visual Studio 2012 3.20.1 2.8.12.2 2.5.1 TM i7 920 @ 2.67 GHz R Intel Core 8 GB

in DICOM format. The computation time drops to about 12 seconds when reading input and writing output are excluded. As previously mentioned in section 3.3.4, adding more seeds increases the computation time because more fuzzy scenes have to be calculated and combined. This effect is illustrated in figure 4.1. The minimum amount of seeds required is two (one for each side). Each seed that is added afterwards increases the computation time by a little more than 1 second.

Figure 4.1: Computation time in function of the number of seed points. The effect that adding more seeds has on the computation time is displayed. Reading and writing the volume is included in these measurements. We believe the computation time can be reduced even further by decreasing the size of the images that are used as input for the filters. We currently use the entire volume size for all filters while we are only interested in the regions around the mandibular nerve canal. Setting a smaller region of interest for each filter should further reduce the computation time.

4.3. METHOD

4.3

31

Method

We had 3 CT and 5 CBCT datasets on which we could test our method. Unfortunately, none of our datasets were annotated so it is difficult to measure the accuracy of our method. The seed points are supplied by the user. Their amount and location also both have an impact on the accuracy. This means we currently have no hard data regarding the accuracy of our method. Figures 4.2 and 4.3 display some slices of the resulting volume after applying our method on a CBCT volume. In figure 4.2(a) a slice of the original volume is displayed. There is quite some noise present. In the other figures the original volume is displayed in red and the pixels that are thought to be part of the mandibular nerve canal are displayed in black. In figure 4.2(b) we can see part of the right nerve is succesfully found. However, in this slice our method has failed to find a part of the left nerve. To solve this we could add a seed point near this slice or decrease the threshold of the fuzzy-connectedness for this side. Figure 4.3(a) clearly shows a tubular structure on both sides. On the left side our interpolation has failed to fill in all the gaps. This means we should best add another seed point near that gap. On the right side it seems we have succesfully added a nice path, but some pixels that shouldn’t have been added are part of the solution. To solve this we should probably increase the threshold of the fuzzy-connectedness step for this side. Finally in figure 4.3(b) the tubular structure is still visible on the right side, but too many pixels are added near the front. This is because the interpolation is pulling it towards the mental foramen. On the left side only a small part of the canal is found. This is another indication that we should decrease the threshold of the fuzzy-connectedness for the left side. Overall, we can conclude that our method can succesfully segment large parts of the mandibular nerve canal. However, to improve the results we usually have to tinker with some parameters or add some seeds. An automated generation of seed points would be a welcome addition to our method. When suitable seed points are automatically generated, the thresholds for both sides could be increased. This would result in less errors and thus better results. Our interpolation implementation also doesn’t always yield the best results. Adding more automatically generated seeds should also increase the quality of the interpolation, and at the same time reduce the need for interpolation. Regardless of this, we believe an improvement of the implementation would benefit the results. It is hard to sample the results of the interpolation at the right points.

32

CHAPTER 4. RESULTS

(a) Slice of original volume.

(b)

Figure 4.2: Results: CBCT volume. Some slices of the resulting volume after applying our method on a CBCT volume.

4.3. METHOD

33

(a)

(b)

Figure 4.3: Results: CBCT volume (continued). Some slices of the resulting volume after applying our method on a CBCT volume.

34

CHAPTER 4. RESULTS

Chapter 5

Conclusion The goal of this thesis was to develop an algorithm enabling the visualisation of the mandibular nerve canal. This meant developing and implementing a method for its segmentation in CBCT images. We reviewed literature about this subject and found a very effective method for CT images that uses fuzzy-connectedness. Then we developed our own method that uses this fuzzy-connectedness and can be applied on CBCT images. To achieve this we try to counter the negative aspects of CBCT images via a combination of a smoothing filter, a gradient filter and a thresholding operation to reduce noise, enhance edges and become less dependent on intensity levels. After these steps the volume should be ready to be processed by the fuzzyconnectedness operation. To make the results of the fuzzy-connectedness operation more reliable we apply interpolation on them. We fully implemented this method from scratch using ITK and visualised the results with MeVisLab. Our method, which includes reading and writing the volume in DICOM format, has a computation time of about 20 seconds on our test system. When reading input and writing output are excluded, the computation time drops to about 12 seconds on our test system. It can be applied on both CT and CBCT images. The user should supply as many seed points as he/she considers necessary. The amount and location of these seed points have an impact on the accuracy. This combined with the absence of annotated datasets means that we can’t support any claims regarding the accuracy of our method with hard data. However, by looking at the results, we can see that our method can succesfully segment large parts of the mandibular nerve canal. But there is some room left for improvement.

Future Work We believe our method would benefit most from automated seed generation. This would make the method more user-friendly and improve the results. It would reduce the need for interpolation, and improve the quality of the 35

36

CHAPTER 5. CONCLUSION

interpolation where it is still needed. Regardless of this, we believe an improvement of the interpolation method would still benefit the final results. Generating hard data regarding the accuracy of our method is of course also very important. When automated seed generation is added, this would also be easier and more consistent. And finally, we believe the computation time could be reduced by setting a smaller region of interest as input for each filter. We currently use the entire volume size for all filters while we are only interested in the regions around the mandibular nerve canal. Implementing this would improve the computation time even further.

Bibliography De Vos, W., Casselman, J., and Swennen, G. (2009). Cone-beam computerized tomography (cbct) imaging of the oral and maxillofacial region: a systematic review of the literature. International Journal of Oral and Maxillofacial Surgery, 38(6):609–625. Gonzalez, R. C. and Woods, R. E. (2002). Digital image processing. Prentice Hall. Hanssen, N., Burgielski, Z., Jansen, T., Lievin, M., Ritter, L., von RymonLipinski, B., and Keeve, E. (2004). Nerves-level sets for interactive 3d segmentation of nerve channels. In Biomedical Imaging: Nano to Macro, 2004. IEEE International Symposium on, pages 201–204. IEEE. Johnson, D. and Moore, W. (1997). Anatomy for Dental Students. Oxford University Press, Incorporated. Kainmueller, D., Lamecker, H., Seim, H., Zinser, M., and Zachow, S. (2009). Automatic extraction of mandibular nerve and bone from cone-beam ct data. In Medical Image Computing and Computer-Assisted Intervention– MICCAI 2009, pages 76–83. Springer. Kim, G., Lee, J., Lee, H., Seo, J., Koo, Y.-M., Shin, Y.-G., and Kim, B. (2011). Automatic extraction of inferior alveolar nerve canal using featureenhancing panoramic volume rendering. Biomedical Engineering, IEEE Transactions on, 58(2):253–264. Kondo, T., Ong, S., and Foong, K. W. (2004). Computer-based extraction of the inferior alveolar nerve canal in 3-d space. Computer methods and programs in biomedicine, 76(3):181–191. Kroon, D.-J. (2011). Segmentation of the mandibular canal in cone-beam CT data. University of Twente. Llor´ens, R., Naranjo, V., L´ opez, F., and Alca˜ niz, M. (2012). Jaw tissues segmentation in dental 3d ct images using fuzzy-connectedness and morphological processing. Computer methods and programs in biomedicine. 37

38

BIBLIOGRAPHY

Lobregt, S., Schillings, J., and Vuurberg, E. (2001). Dental implant surgery: planning and guidance. Medicamundi, 45(4):30–35. Moris, B., Claesen, L., Sun, Y., and Politis, C. (2012). Automated tracking of the mandibular canal in cbct images using matching and multiple hypotheses methods. In Communications and Electronics (ICCE), 2012 Fourth International Conference on, pages 327–332. IEEE. Mortensen, E. N. and Barrett, W. A. (1995). Intelligent scissors for image composition. In Proceedings of the 22nd annual conference on Computer graphics and interactive techniques, pages 191–198. ACM. Rueda, S., Gil, J. A., Pichery, R., and Alca˜ niz, M. (2006). Automatic segmentation of jaw tissues in ct using active appearance models and semiautomatic landmarking. In Medical Image Computing and ComputerAssisted Intervention–MICCAI 2006, pages 167–174. Springer. Sotthivirat, S. and Narkbuakaew, W. (2006). Automatic detection of inferior alveolar nerve canals on ct images. In Biomedical Circuits and Systems Conference, 2006. BioCAS 2006. IEEE, pages 142–145. IEEE. Stein, W., Hassfeld, S., and Muhling, J. (1998). Tracing of thin tubular structures in computer tomographic data. Computer Aided Surgery, 3(2):83–88. Suetens, P. (2009). Fundamentals of Medical Imaging. Cambridge medicine. Cambridge University Press. Tustison, N. J. and Gee, J. C. (2006). Generalized nd c k b-spline scattered data approximation with confidence values. In Medical Imaging and Augmented Reality, pages 76–83. Springer. Udupa, J. K. and Samarasekera, S. (1996). Fuzzy connectedness and object definition: theory, algorithms, and applications in image segmentation. Graphical models and image processing, 58(3):246–261. Yau, H., Lin, Y., Tsou, L., and Lee, C. (2008). An adaptive region growing method to segment inferior alveolar nerve canal from 3d medical images for dental implant surgery. Computer-Aided Design and Applications, 5(5):743–752.

FACULTY OF ENGINEERING TECHNOLOGY CAMPUS DE NAYER (@Thomas More) Jan De Nayerlaan 5 2860 SINT-KATELIJNE-WAVER, België tel. + 32 15 31 69 44 [email protected] www.iiw.kuleuven.be