Interactive segmentation of abdominal organs from 3D CT and MRI images

Department of Science and Technology

Institutionen för teknik och naturvetenskap

Linköping University Linköpings Universitet

SE-601 74 Norrköping, Sweden

601 74 Norrköping

LiU-ITN-TEK-A--10/025--SE

Interactive segmentation of

abdominal organs from 3D CT and

MRI images

Olof Hansson

LiU-ITN-TEK-A--10/025--SE

Interactive segmentation of

abdominal organs from 3D CT and

MRI images

Examensarbete utfört i medieteknik

vid Tekniska Högskolan vid

Linköpings universitet

Olof Hansson

Handledare Chunliang Wang

Examinator Örjan Smedby

Upphovsrätt

Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare –

under en längre tid från publiceringsdatum under förutsättning att inga

extra-ordinära omständigheter uppstår.

Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner,

skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för

ickekommersiell forskning och för undervisning. Överföring av upphovsrätten

vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av

dokumentet kräver upphovsmannens medgivande. För att garantera äktheten,

säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ

art.

Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i

den omfattning som god sed kräver vid användning av dokumentet på ovan

beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan

form eller i sådant sammanhang som är kränkande för upphovsmannens litterära

eller konstnärliga anseende eller egenart.

För ytterligare information om Linköping University Electronic Press se

förlagets hemsida http://www.ep.liu.se/

Copyright

The publishers will keep this document online on the Internet - or its possible

replacement - for a considerable time from the date of publication barring

exceptional circumstances.

The online availability of the document implies a permanent permission for

anyone to read, to download, to print out single copies for your own use and to

use it unchanged for any non-commercial research and educational purpose.

Subsequent transfers of copyright cannot revoke this permission. All other uses

of the document are conditional on the consent of the copyright owner. The

publisher has taken technical and administrative measures to assure authenticity,

security and accessibility.

According to intellectual property law the author has the right to be

mentioned when his/her work is accessed as described above and to be protected

against infringement.

For additional information about the Linköping University Electronic Press

and its procedures for publication and for assurance of document integrity,

please refer to its WWW home page: http://www.ep.liu.se/

Abstract

Within the medical field, image segmentation is an important tool which can be used by radiologists and surgeons who want quantitative measurements of a lesion or organ. To be clinically useful, the tool has to be fast and easy to use. This work comprises implementation of the image foresting transform for segmentation using the Dijkstra algorithm and compares computation time between the implemented algorithm and a previous implemented algorithm, Bellman-Ford. These algorithms solve the shortest path with minimum cost problem. For a given cost function, similar results both in computation time and visual results were obtained with the two algorithms. Changing the cost functions, on the other hand, yielded very different segmentation results. The volume of liver and kidney was compared with manually delineated organs regarding seed planting and execution time. A graphical user interface has also been implemented.

Acknowledgments

First of all I would like to thank my excellent supervisor, Chunliang Wang, Ph.D student at CMIV for his encouragement, patience and support to me during this thesis. He has been a major influence with his sense of humour, positiveness, joy and knowledge.

Thanks also to my examiner Professor Örjan Smedby, who gave me this opportu-nity and has come with straight suggestions whenever I felt vacillating. Thanks to all CMIV colleagues for a warm and open enviroment making my thesis time enjoyable.

I would also like to thank my family and friends for supporting me and believe in me through all these years.

There is a particular friend that I want to thank. Avraz Hirori, originially from Kurdistan, but now living in France. With a rare intellect he has influenced me most of everybody I have met during the last six years of my life in Norrköping and Linköping. The stories about his life were breathtaking and exciting, remembering me how easy life is in Sweden compared to other countries in the world. He gave me new persepectives of life by which one could see and distinguish which parts in life that are most important. He also learned me about choosing hard decisions in life and how to take care of them.

I miss our discussions about life and culture.

I miss our nights at school when studying math until late meanwhile having a great time and for the first time really understanding math.

I miss his joy, fantasy and energy. I miss my friend.

Finally, I thank my lovely Marie for supporting me through all this work. She has been an enormous energy source with her positiveness, joy and love. I will never forget your face expression when carefully riding your new bike.

For the first time in my life I really understand the change mathematics and physics can do. I think about the people that have contributed to this part of science before me. This work is built upon that knowledge.

Contents

2.1.1 CT . . . 6

2.1.2 MRI . . . 7

2.2 Segmentation . . . 8

2.3 Minimum-cost path search . . . 10

2.3.1 Simple case - Dijkstra . . . 11

2.3.2 Simple case - Bellman-Ford . . . 13

2.3.3 2D case . . . 14

2.3.4 Connectivity map . . . 14

2.3.5 Label map . . . 15

2.3.6 Moving on to 3D . . . 15

2.4 Cost of the path . . . 16

2.5 Cost of edges . . . 16

3 Implementation 19 3.1 Minimum-cost path search algorithms . . . 19

3.2 I/O of the Image Foresting Transform Module . . . 20

3.2.1 Version 1 . . . 21

3.2.2 Version 2 . . . 22

3.3 An application for abdominal organ segmentation . . . 23

4 Results 25 4.1 Computation performance . . . 25

4.1.1 Segmentation results - Dijkstra vs Bellman-Ford . . . 26

4.1.2 Cost functions evaluation - Kidney . . . 27

4.1.3 Cost functions evaluation - Liver . . . 28

6 Contents

5 Discussion 31

5.1 Computation results . . . 31

5.2 Visual results of Dijkstra and Bellman-Ford . . . 32

5.3 Cost functions and volume calculation . . . 32

5.4 Future work . . . 33

5.4.1 Improving Dijkstra . . . 33

5.4.2 Improving user interface . . . 33

5.4.3 Cost functions . . . 35

5.4.4 Conclusion . . . 35

Bibliography 37 A User interface manual 39 A.1 The interface . . . 39

Contents 7

List of Tables

2.1 Table of Dijkstra algorithm . . . 12

2.2 Table of Bellman-Ford algorithm . . . 13

4.1 Seeding and computation time . . . 26

4.2 Results when subtracting Dijkstra from Bellman-Ford segmented data 26 4.3 Results when calculating volume of kidney with different cost func-tions . . . 27

4.4 Volume of kidneys manual segmented . . . 28

4.5 Results when calculating volume of kidney with different cost func-tions . . . 29

4.6 Volume of liver manual segmented . . . 29

List of Figures

2.2 The label map - in theory . . . 15

3.1 The Connectivity Map . . . 20

3.2 The Label Map . . . 20

3.3 The 2D interface part . . . 24

3.4 The 3D interface part . . . 24

4.1 Comparison of a segmented slice using different algorithm . . . 26

4.2 Comparison of segmented slice using different cost functions . . . . 27

4.3 Comparison of segmented slice of liver using different cost functions 29 A.1 The 2D interface part . . . 39

Chapter 1

Introduction

1.1

Problem description

During the last ten years there has been a substantial and rapid development of medical image modalities like CT and MRI. The best scanners today produce thousands of 2D images from one abdominal scan with a resolution as small as 0.3 mm. This leads to huge amounts of data to handle. A big advantage of the high resolution images is that it is now possible for physicians to inspect organs in detail but also to visualize the lesion from different views. The disadvantage, on the other hand, is that radiologists need to interpret huge amounts of data. This is a challenge, both regarding time for inspection of the data and finding methods for presenting the results to clinicians without showing all data.

Traditional visualization of the data by displaying 2D transverse images in a grid on screen or scrolling through the 2D slices of the volume by paging up and down is yet a reliable and accurate method in most cases. Even though this method is still the first choice of many radiologists in daily practise, it becomes more and more time consuming as the amount of data increases from 20-30 images per pa-tient to several hundreds, and more recently up to 2000-3000 images (1-1.5GB) per patient. This forces medical image scientists to find new solutions to improve the efficiency of clinical reviews. Common visualization techniques like Maximum Intensity Projection (MIP) and Direct Volume Rendering (VRT) were introduced to partly overcome the limit of the 2D screen by construction a 3D view pro-jected to the screen instead and this has been succesful in serveral fields [1] [2]. However, it is still problematic to use these techniques in areas where anatomy is complicated, for example the concealed neural vessels in the skull and the circu-lation system enclosed by the ribs in the thoracic region. Due to the overlapping of the organs in many angles, this makes it hard to interpret the 3D visualized data in the relevant regions. A rational solution to this would be to isolate the organs of interests first before visualizing them. On the other hand, in certain cases, such as for surgery planning, it is often necessary to have a quantitative measurement of the dimension or volume of the organ or lesion. In both cases it

2 Introduction

requires a tool that segments the targeted organ or lesion [3]. Although this can be done manually by marking the organs’s contour slice by slice in the volume, it is usually time-consuming and labor-intensive. An automatic or semi-automatic segmentation is desired to speed up this procedure with much less user interaction.

There are many segmententation methods published [4] [5], ranging from easy thresholding techniques to more advanced ones like neural network, from graph based algorithms to PDE (Partial Differential Equation) based solutions such as Fast Marching and Level Set technique. All produce good results, but mostly only in special cases. Often the trivial solutions are considered to have better performance in term of speed, but these are often sensitive to image quality. More sophisticated methods may give more reliable results in noisy and blurry images, but they require longer computation time. In order to balance the requirements on speed and robustness of the algorithm, this project will concentrate on a class of methods called the Image Foresting Transforms [5] as the foundation to the abdominal parenchymal organ segmentation.

The Image Foresting Transform unifies a family of algorithms where the main idea is to create segmented areas as spanning trees from the seed points (which the user gives) in a graph (2D/3D image). The construction of these trees is based on the minimum cost path search which gives as a result that all the leaves (pixels) are connected to a seed point where the cheapest path was found. In turn, this method can work as a region-based method as well as a boundary tracking method given a sequence of pixels along the object. To be used in different applications, the algorithm can be changed by using different cost functions to construct the tree. Using a binary cost function result in an algorithm similar to region growing, a discretized cost function results in an algorithm like fuzzy connectedness, and using analog accuracy cost function as in Dijkstra results in an algorithm as Fast Marching. This thesis present results using the Image Foresting Transform for abdominal organ segmentation.

1.2 Aim 3

1.2

Aim

The goals of this thesis are:

• To develop a segmentation method based on the Image Foresting Transform with the shortest minimum-cost path search approach.

• To compare the Dijkstra algorithm with a modified Bellman-Ford algorithm, both in computation time and visual results.

• To compare different cost functions, both in visual results and in quantitative measurements.

• To compare the calculated quantitative results from the Dijkstra algorithm with corresponding results from manual delineation of abdominal organs. • To build a user interface which has the possibility to give quantitative

mea-surements and to give radiologists’ a program to train with.

1.3

Reader prerequisites

The reader is expected to have some kind of engineering background. Knowledge about image processing, basic visualization techniques and data structures needs to be fresh.

Chapter 2

Background

2.1

Medical images and Acquisition techniques

In this report, medical images mainly refer to CT and MR images. Like all digital images they consist of small units called pixel in a 2D image and voxel in a 3D volume. The numerical value in this unit describes different properties from the original scanned data depending wether it was scanned using CT or MRI. For example, in CT the pixel or voxel describes the density of the scanned object. While the pixel is often a mean within the area of the pixel, in image processing it is often calculated upon the centre point of the area. One can see the image as a 2D or 3D Cartesian grid where the centre points of the pixels or voxels are the grid points. The distance between points in this grid thus depends on the pixel or voxel size, this is called spacing. Mostly three spacing parameters (x,y,z directions) are needed for image processing. Due to limitations of the scanners, CT- and MR-scanner will produce many 2D slices of the volume. This implies that the spacing along the z-direction is usually not the same as in the x and y directions of spacing (z spacing is usually larger). In some image processing techniques one needs to have isotropic voxel size as input. This forces the user to do a resampling of the volume with an apprioriate interpolation method. An overview of the MRI and CT techniques is provided in the next sections.

6 Background

2.1.1

CT

Computed Tomography (CT) is a technique to scan the human body with X-rays which was introduced in the late 1960s. The primary data aquisition part is sim-ilar to ordinary X-ray technique. The detector receives an amount of radiation depending on how much of the radiation is absorbed in the body. The absorption depends on which organs the rays have gone through: the denser the organ is, the more X-ray absorption it is. But unlike traditional X-ray technique, which acquires an image in one direction, in CT X-ray source and detector rotate around the patient and produce images in many directions that cover 180 degrees of the scanning field. Due to limits in the technique, most CT scanners only capture data with an array of detector on the rotation plane. In each direction it produce one row data if the scanner has one array of detectors, and several rows of data if the scanner has more arrays of detectors (in current models up to 320 rows). The whole rotation then produce a 2D image called sinogram [6]. This projection pro-cedure is known as the Radon transform in mathematics. An inverse of the Radon transform gives the CT image of the scanned slice. There are many different types of reconstruction algorithms but the most common one is filtered back projection. CT is usually used for tissue examination such as in the abdomen, pelvics and brain. The CT value is calibrated by defining air as -1000 and water as 0 H. Sometimes iodine-rich solutions are injected into patient’s vascular system which enhance the blood CT value to 300-700 HU and increases the contrast between vessels and surrounding soft tissue (whose CT value remains within 50-150 HU). This type of exam is known as a Computer Tomography Angiography (CTA) [1], and is a good way of showing the vessels, to trace them, and to analyze if there is any stenosis in the arteries. Despite wide use of CT exams, radiation generated during CT scan can harm cells in the body and increase the risk of developing cancer. Physicians must balance between the benefits and this drawbacks.

2.1 Medical images and Acquisition techniques 7

2.1.2

MRI

The technique behind MRI, Magnetic Resonance Imaging, is quite different com-pared to CT. Magnetic resonance, also referred to as Nuclear Magnetic Resonance (NMR) is a phenomenon where the nuclei of certain atoms absorb energy from an applied electromagnetic pulse of a certain frequency when placed in a strong magnetic field. The energy absorbed is emitted back in form of an electromagnetic wave when the external electromagnetic pulse is turned off. As different atoms’ nuclei respond to different frequencies, an MRI scanner is normally designed to excite one type of nuclei, hydrogen is most common. The signal depends on the concentration of this type of nuclei as well as its physical properties that differ between water, fat and edema. A challenge in MRI is to distinguish which part of the body is sent with the signal. When the patient is in the magnetic field all the nuclei in the body contribute to the signal picked up by the reciever. This is solved by having a gradient magnetic field on top of the strong homogeneous magnetic field. Because nuclei at different magnetic field respond to different reso-nance frequencies, now only one 2D slice in the volume will respond to the external electromagnetic pulse, which is called slice selection. This gradient magnetic field is then turned off when the external pulse stops. Two variable gradient fields are then applied parallel to the 2D slice, resulting in different parts of the 2D plane generating radiofrequency responses at different frequency and phase, which is called frequency and phase coding, respectively. After this encoding, the final 2D image can be produced using the inverse Fourier Transform. Depending on the timing of reading out the signal (so called MR sequence), the image magnitude can reflect different properties of the resonance nuclei, such as T1 value and T2 value. An image that reflects mainly the difference in T1 value inside the body is called a T1-weighted image, and one that reflects the difference in T2 value is called a T2-weighted image. At the MR scanner, a large number of parame-ters can be changed for a specific examination and the images produced are often relatively noisy. Although the method MRI is a very multifaceted technique to scan the human body because of its superior soft tissue contrast compared to CT, it is widely used for diagnoses of tumours in the brain as well as other parts of the body, Multiple sclerosis (MS), stroke, infections, in particular in the spine and brain, myocardial ischemia, parenchymal diseases in the abdomen, musculoskeletal disorders and other conditions.

8 Background

2.2

Segmentation

Segmentation is a part of image processing which deals with partitioning an image into several regions of interests depending on different properties in the image for example intensity, gradient, spatial and temporal position and variance. The use of segmentation is multifarious and the technique is important in many different fields today, especially in the medical field. Segmentation is used in various appli-cations such as volume calculation of organs or lesions, preprocessing of data before visualizing it in 3D, locating objects in different kinds of images but it can also be used within industries to check whether a product has the right specification. It can be used to study anatomical structures like blood vessels, muscles, bones and brain areas and more. It is used within face and finger recognition systems.

There exist three different segmentation methods or classes - pixel based, region based, and border based methods. These methods differ regarding speed, ro-bustness and complexity. Pixel based methods are methods such as thresholding and classification. Region based methods include Region growing, Watershed and Level sets. Border based methods include Live wire and Active contours (snakes) and are used to find edges between objects. This report concentrates on region based methods.

To be able to segment images from the CT and MRI scanner one has to have a robust method because of the noise and blur in the images. In this report, we will build on earlier work dealing with segmentation of vessels in MRA and CTA [1] and [2]. In these both cases a method called Fuzzy Connectedness was used to segment the vessels. This method is suitable when segmenting images where the objects have high intensity compared to the background, as the vessels in CTA. However, the Fuzzy Connectedness approach used in [1] using the minimum inten-sity along the paths as cost will most likely fail to segment organs in the abdomen, where the organs are not surrounded by lower-intensity regions. Therefore, an-other approach is necessary.

In [5], Falcao, Stolfi, and Lotufo describe the Image Foresting Transform. This transform can be modified into several segmentation algorithms by a small change in the algorithm. The Image Foresting Transform defines a relation between graphs and images. By considering the vertices as pixels and arcs as adjacency relations between the pixels, one can create spanning forests from start positions in the algorithm. The path-cost along the arcs will from now on be called cost function. It usually depends on image intensity, gradient and position in the pixels along the calculated path or scale property. The starting points are in this project where the user adds the seed points. An algorithm is needed to solve the problem of calculating the path from the pixels to the user defined seed points. Dijkstra and Bellman-Ford are algorithms that solve this problem. With the help of these algo-rithms one can trace a pixel to its nearest and cheapest starting point (seed point) and label the pixels along the path depending on the seed point type. Doing this calculation for the whole image this will give rise to segmented areas which can

2.2 Segmentation 9

be seen as forest-like areas. Although by thresholding this spanning trees from a single seed one can get a segmentation result, a more natural way is to let the user define two or more different types of seed. Comparing the cost from one pixel to different seeds, the membership of this pixel can be determined using knowledge about the lowest cost of the surrounding pixels which implies more probability that they are related to the seed point. Using minimum-cost path search with Dijkstra and Bellman-Ford one can implement these algorithms to work in a general way by letting the user change the cost calculation along the path instead of using the minimum intensity along the path as in [1] and [2]. With this approach, different segmentation results can be obtained depending on the form of the cost function. However, one has to analyze how these algorithms work in order to segment with them. A brief description of the two algorithms minimizing the cost along paths are described in the next sections.

10 Background

2.3

Minimum-cost path search

Considering graphs, minimum-cost path search deals with finding the shortest path between nodes or vertices at minimal cost.

A graph G = (V, E) consists of a set of vertices V = v1, v2, v3, v4....vn 6= ø and a

set of edges E = e1, e2, e3, e4....en ⊆ V × V . Every edge e is assigned a weight or

in other words the cost from vertex u to vertex v. In this project we will consider directed graphs where the edges e between vertices are ordered pairs of vertices (u, v) such that u 6= v. When going from a pixel to another one also considers a cost between them, which will be defined in different ways while considering digital images later in the report. The definition of a weighted graph is that here, every edge e is a triple (u, v, w) where u and v is a pair of vertices and w is a weight or in other words the cost going from u to v. An important definition in this context is the definition of an sparse and dense graph. A graph is sparse if m = O(n) and a graph is dense if m = O(n2_{) [7] where m is the total number of edges in the graph.}

A path p in the graph G from vertex v1to vkis a sequence of vertices v1, v2, v3, v4....vk

such that (vi, vi+1) ∈ E for 1 ≤ i < k − 1. The weighted path length for p is

ac-cording to equation 2.1 the sum of all the costs along the path. Later we will consider an alternative equation for finding the minimum-cost of a path.

k−1

X

i=1

wi,i+1 (2.1)

This basically means summing the costs along the path. The minimum-cost path search in graph G amounts to solving a path problem between two given points with minimum cost, in other words, the shortest path with minimum weight (cost) path from a vertex u to a vertex v. The Dijkstra and Bellman-Ford algorithm both solve this problem. This general problem is widely addressed in many different fields today. The problem is this: having a start vertex s ∈ V , one wants to find the shortest minimum weighted (cost) path from s to every other vertex in V . These algorithms are very general in their nature and are used in many differ-ent fields today. Apart from image segmdiffer-entation, the algorithms are extremely important in computer networks where the router decides which way to send the package according to a table. Many algorithms have been proposed to solve this minimum-cost problem. Among them, the Dijkstra and Bellman-Ford algorithms are most used because of their efficiency. The next two sections describe these two algorithms in more detail.

2.3 Minimum-cost path search 11

2.3.1

Simple case - Dijkstra

Figure 2.1: A graph system - a simple case

In 1959, the Dutch computer scientist Edsger Dijkstra conceived an algorithm which calculates the shortest path from a single source to every other vertex in a graph by iteratively selecting the closest point from the neighboring vertices [8] [9] [10].

At initialization the algorithm will apply infinite distance or cost to every vertex in the graph. These distances are then compared with the accumulated distance or cost along the actual path. Consider a start vertex, V3 in figure 2.1, the algorithm first assigns zero distance at V3, then updating the distance to neighbor points (V5,V4) with the costs on the edges. The closest point (V5) can then be easily found, the shortest path from V3 to V5 is then secured, in other words this value will be permanent as there will be no path that can be shorter than the current path. The path V3-V5 will be the shortest path. After the distance to V5 is secured the distance to V7 will be updated (to 10) by adding the edge cost (5) on top of distance at V5(5). Also V2 is updated by adding the cost from V5(3) giving it the total cost of 8. V2 is then elected from the new neighboring vertices (V7, V2, V4) using the updated value. The distance to from V3 to V2 (8) is secured and then used to update its neighbors distance (V7: 10 to 9, V4: 10, V1:∞ to 12 and V8:∞ to 15). Doing this procedure again will results in updating the V8:15 to 13 by going through vertex V6 instead of V2. The algorithm will eventually go through all the vertices in a greedy way until all the points are secured or the desired point are reached. The worst performance of Dijkstra algorithm is

O(m + n ∗ log(n)) [11] where m is the total number of edges and n number of

12 Background

it also saves the previous vertex the shortest path came from and eventually this can be used to track the shortest path backwards from any secured points in the graph to the source point, like figure 2.2. Checking the minimum neighbor and saving the previous vertex is done by the code below. u is the neighbor vertex and v is the current vertex and w the weight or cost. This also shows how the comparison and update is done when neighbors to a vertex have infinite distance or cost at initialization.

if (dist(u) > dist(v) + w(v,u)) {

dist(u) = dist(v) + w(v,u) path(u) = v

}

When the algorithm has gone through all the vertices we get information with distances and paths like table 2.1. This information is actually all that one needs

Table 2.1: Table of Dijkstra algorithm

V 1 2 3 4 5 6 7 8

dist/cost 12 8 0 10 5 11 9 13

path 2 5 3 3 7 2 6

to trace a path from one vertex to another. Investigating the table we can see that every vertex has a path - the vertex we come from. With this information one can recursively go back until one comes to the starting point where distance = 0. That is described in table 2.1.

2.3 Minimum-cost path search 13

2.3.2

Simple case - Bellman-Ford

The Bellman-Ford algorithm [12] [13] is proposed to solve the same single source shortest path problem as the Dijkstra algorithm does [14]. Although it is also based on iteratively updating the distance value at each point, it does not distin-guish the secured vertices, neighbor vertices and far away vertices. All the edges will be checked in every iteration. Consider the same system as in figure 2.1. In the same way as in the Dijkstra algorithm the chosen start vertex V3 is assigned a distance of 0, and all the other points are assigned a distance of infinity. Then algorithm starts to relax all the edges in the graph. Bellman-Ford will now relax (or update) all vertices in the graph. In the first iteration, vertices V5 and V4 will be updated as the distance value at V3 plus the cost of the edges to those vertices are less than the initial value (infinite). The minimal cost to go to V5 and V4 is now 5 and 10. It should be pointed out that all the edges in the graph are checked. The distance value at V1 might also be changed if the edge cost 4 and 6 are checked after V2 or V4’s value are updated. This is also the case for V7 and V8. Even V6 might be updated in this iteration. But in the worst case only V5 and V4’s value will change while the other points remain infinite. In the second iteration the distance values at V5 and V4 will not change anymore. Then the algorithm goes on to check all the other edges and update their further points value if necessary.

The algorithm does not care if one vertex has reached the minimum-cost from the start node, it checks all edges and vertices in each iteration. One has to do this relaxing part n − 1 times, where n is the total amount of nodes in the system to secure the convergence of the shortest path at each point. This results in a performance of O(m ∗ n) [14]. The Bellman-Ford algorithm is sensitive to dense graphs which implies more scans due to the large number of edges in the graph. However in reality, the convergence of the algorithm can be reached far earlier than n − 1 times of relaxing the edges. The real convergence can be determined when no distance value is changed in one iteration. Still, the performance is not comparable to Dijkstra algorithm as every edge has to be examined in each itera-tion. In [1], Wang and Smedby presented a variation of Bellman Ford algorithm, in which only the changed points, neighbors will be checked in the next iteration. This results in a performance of O(m ∗ (1 − qn)/(1 − q) where q is the ratio of the number of changed point compared to the last iteration and n is the iteration times. The table of the minimum costs in Bellman-Ford algorithm is shown in table 2.2.

Table 2.2: Table of Bellman-Ford algorithm

V 1 2 3 4 5 6 7 8

dist/cost 12 8 0 10 5 11 9 13

14 Background

2.3.3

2D case

In the example in figure 2.1 some of the vertices has a number of vertices as neighbours. The vertex V6 only has one possibility to go to, V2 instead has four different vertex possibilities. Consider now a 2D image with width x and height y. Let each vertex now denote each pixel in the image. For every pixel in the image which is one pixel from the border it exists eight pixel neighbours. That makes the calculation part many times more computationally demanding. As described earlier, one can notice the same situation as vertices V1 and V6 which has arrows to each other. This will happen to every pixel and to its neighbouring pixels in the image.

2.3.4

Connectivity map

Consider a 2D image where the nodes is the pixels in the image and the arcs are a adjacency relationship between them. The user adds different types of seed points in the image. Every seed point is basically one voxel where the distance is set to zero in both the Dijkstra and the Bellman-Ford algorithm. In contrast to the simple example above, with one starting point, now there will be as many starting points as the seed points the user adds. The algorithms update all the distances to all pixels in the image in their special way. It is these distances that are showed in the Connectivity Map and, depending on the type of cost function applied, the result will be very different in each case. Although the Connectivity Map itself from a visual perspective is not that important, its information is and the Label Map is going to use this information to create the forest-like segmented parts of the image [5]. Once looking at a Connectivity Map one would recognize the shifting of distances depending on which cost function has been used.

2.3 Minimum-cost path search 15

2.3.5

Label map

From the Connectivity Map every pixel has an accumulated minimum distance value depending on how many and where the seeds points are in the image. The distance in every voxel is calculated as the shortest distance at minimum-cost from the specific seed point. All the voxels with the minimum distances to an endpoint (seed point) form a path. The Label Map is constructed by inspecting the paths backwards from every pixel to a seed point. Depending on the seed point type we assign those voxels a label. If this is done for all pixels in the image, one will get a forest-like segmentation [5] where the areas with the specific labels can be seen as the forests and the seed points can be treated as the roots to this forests.

Figure 2.2: The label map - in theory

2.3.6

Moving on to 3D

Now add depth to the image and let each vertex denote a voxel in the 3D volume. Let the arcs between a voxel neighbors be an adjacency relationship of the voxel and its 26 neighbors. A volume of size 512x512x90 contains ≈ 23,6 million voxels. In this domain, the structure of the algorithm and the way to solve the problem will show how sensitive the algorithm is to the dramatic increase of data calculations even when running on modern computers. The algorithm used calculates the Connectivity Map which stores all minimum-cost accumulated distance in every voxel in the volume. Later, one traces the path using the path array which stores the previous voxel index. This is done from every voxel in the image until one ends up in a seed point. Tests in this report show that the computation time of the algorithms not only depend on the size of the data but also on the placement of seed points.

16 Background

2.4

Cost of the path

Conventional minimum cost path searching algorithms calculate the cost of a path by accumulating the total cost of all the edges on the path. However this need not to be the only choice when implementing image segmentation. Often one is interested in the minimum or the maximum value on the path. This introduces a new subcategory of image foresting transform methods - fuzzy connectedness [4]. The advantage of the method is that the total cost of the path will not decay as the nodes on the path increases. This makes it best suitable for specific applications such as for vessel segmentation in CTA [1] [2]. Another variation of calculating the cost of a path exist using a binary condition. If all points on the path is within an upper and lower threshold it is true, otherwise it is false. This results in another method known as region growing. The binary condition dramatically reduces the computation load and allows generating results in almost real-time. However, this type of method is very sensitive to noise and can be only used in the condition where the noise level is relative low. In this project equation 2.2 is used for calculating the strength of a path p.

Sf(p) =

X

wi∈p

(f (g(wi), g(v))) (2.2)

The strength of a path p with aspect to f is the sum of the minimum costs along the path p with the cost function f (g(wi), g(v)) where v is the centre voxel in the

neighborhood and wi a neighboring pixel.

2.5

Cost of edges

The cost function is the function that decides the cost from going to a pixel v to a neighboring pixel wi. It is normally calculated from image properties such

as intensity, gradient, position, scale or an combination of those factors. A few common cost function is listed below:

1. Intensity - Intensity itself or its inverse can be used as a cost function. This will make the trees attempt to stay inside a dark area (or bright area if inverse of intensity is used). Different trees are then separated by the bright ridges (watershed) in the image. This type of cost function is implemented in a fuzzy connectedness framework mentioned above and was used by Wang and Smedby in [1] for CTA segmentation where they try to separate contrast enhanced vessels (bright in CTA) from bone structure (bright) and contrast enhanced heart chambers (also bright).

2. Intensity difference - Two types of intensity difference are often used to generate cost function for image foresting transfer - the intensity difference between two neighbor pixels and the intensity difference from any pixel to the seed point. Notice that the absolute value of the difference is used, since the Dijkstra algorithm can not handle negative costs. The function g denotes the gray scale function. Three types of such cost functions are used in this project.

2.5 Cost of edges 17

Absolute value of difference

f (g(wi), g(v)) = |(g(wi) − g(v))| (2.3)

Exponential of the absolute difference

f (g(wi), g(v)) = e|(g(wi)−g(v))| (2.4)

The difference squared

f (g(wi), g(v)) = (g(wi) − g(v))2 (2.5)

These cost functions will make the tree stay in an area where the pixel intensities are relatively similar. Compared to using the intensity of one pixel as mentioned in [4], this cost function suites more applications, but adds more computation time.

3. Gradient - Gradient is an important operator in image processing which can be used as cost function. The gradient is a vector which points in the direction where the rate of change is greatest. Considering a function f , the gradient of f is given by equation 2.6.

∇f (x, y, z) = (df dx, df dy, df dz) (2.6)

In medical images, a pixel with greater gradient usually means that it is on the boundary of an organ. A boundary tracking method like Live Wire uses the inverse of gradient as the cost function to attract the minimum path to the edges [15].

4. Relative position - Cost function based on the positions and distance to the specific object can be used. By adding the relative position from nodes to nodes along the path, one can calculate the geometric distance of the path. The minimum cost path searching then degenerates into the distance transform. This is useful in many cases for example to locate the centerline of a binary shape [5].

5. Scale - The Scale of a pixel is defined by Udupa [4] as the radius of the largest hyperball centered at that pixel which lies entirely in the same ob-ject region. In practice, the homogeneity measurement inside the ball area is used to judge whether it enters another object’s region. The radius of the ball starts with one and is increased by one in each iteration until the homogeneity measure drops below a certain threshold, then the ball is con-sidered to have entered a significant different region. It is predictable that the scale of pixels on the edges will be small, which can be converted to higher penalty on the cost, and in turn, limits where the trees grow only in relative homogeneity areas. The advantage of scale measure over intensity based cost function is that it avoids the leaking problem by giving a higher

18 Background

penalty on small bridge structures. Consider two objects that have similar intensity and are connected through a small line of noise with the same in-tensity. Using intensity based cost function the trees in one object can easily grow to the other object via the bridges, but using scale based cost function, those small bridges will be punished because of their small scale.

6. High dimensional vectors and tensors - More general structures of the cost function can be built upon theory of vectors and tensors. Although minimum cost path searching algorithms itself only consider scalar value as the cost of an edge, it is possible to describe the difference of vectorial data with a scalar value. An application example is found in [16]. Other high dimension property like tensor structures can also be used as a cost function with proper formula [17]. These advanced cost functions go beyond the content of this thesis, and interested readers are referred to [17]. These cost functions are not always isolated from each other. One can design a cost function from a combination of several different properties by giving them a proper weight factor. In [16], the vector properties, scale information and intensity difference are all combined into one formula.

Chapter 3

Implementation

In this project, an Image Foresting Transform module is implemented in MeVis-Lab as an image processing module. MeVisMeVis-Lab is a visual programming platform for image processing and visualization research and development with a focus on medical imaging [18]. It provides many advanced medical imaging methods for segmentation, registration, volume visualization in terms of modules by which the user is allowed to create her own pipeline by connecting different modules together. In this way, MeVisLab provides fast integration and testing of algo-rithms and development of application prototypes that have potential to be used in clinic enviroments. There are clinical prototypes which have been developed in MeVisLab such as assistance for neuro-imaging, dynamic image analysis, surgery planning and vessel analysis. Its interface was developed using Qt and it worksR

on several platforms such as Windows, Linux and MacOS X. MeVisLab also pro-vides interaction and building user interfaces with help of scripting languages like Python and Javascript to send variables between modules. OpenGL and theR visualization and interaction toolkit OpenInventor are big parts of how the visu-alization is done. MeVisLab is free for non-commercial use and supports various file formats like DICOM, JPEG, TIFF, 3DS (3D studio max), STL, VRML and many more.

3.1

Minimum-cost path search algorithms

Two versions of the Dijkstra algorithm are implemented. In the first version no neighbor queue is implemented. In each iteration, all n points are scanned to find the next closest point, which makes a performance in O(n2_{). The second}

imple-mentation of the Dijkstra algorithm uses a priority queue to store the neighboring zone in an ascending order depending on the cost. In each iteration the cheapest neighbor can be directly located at the head of the queue which means a read-ing time of O(1), and the insert time is O(log(k)) (k is the number of neighbor points in the queue). A modified Bellman-Ford algorithm is also implemented as described in [2]. The implementations of the Dijkstra algorithm are given in pseudocode in section 3.1.1 and 3.1.2.

20 Implementation

3.2

I/O of the Image Foresting Transform Module

The Image Foresting Transform module requires two inputs and two outputs. The first input is the original data such as CT or MRI data sets, the second is a seed map where all seeds represent different integers depending on the different type of seed point. There are different integers for different types of organs to segment. The first output consists of a Connectivity Map containing the accumulated dis-tance values from the Dijkstra algorithm when every point has obtained a final value. As mentioned in 3.1, the distances to all voxels in the volume are saved. The Connectivity Map consists of these distances at their corresponding place in the image. The map helps experienced users to find the optimal location to put the seed points. Another reason is to keep this output when new seeds are added, the distance only need to be partly updated instead of being reconstructed. An example of a Connectivity Map can be found in figure 3.1. The second output is

Figure 3.1: The Connectivity Map

the Label Map, in which the segmentation results is given by assigning the seeds index to every pixel of integer number of the seed in the seed map. The Label Map is calculated from the tree structure output from the image forresting transform module. Recall that the Dijkstra algorithm stores the path, by recording the index to the previous voxel from where the shortest path comes. This can be used to recursively go back along the path and end up at a seed point - the seed which the voxel is more connected to, in other words: this path is cheapest. Then one marks these voxels along the path with a label of the specific seed type. This is done for all voxels in the volume and it will eventually look like the example in the figure 3.2.

3.2 I/O of the Image Foresting Transform Module 21

3.2.1

Version 1

Algorithm 1 Dijkstra algorithm O(n2₎

1: _{procedure Dijkstra(seeds, distancemap, path, visited)} 2: for all v = 0 → sizeOf V olume do

3: initializing distancemap, path and visited arrays 4: end for

5: for all n = 0 → numberOf Seeds do

6: distancemap(seeds(n)) ← 0

7: visited(seeds(n)) ← true

8: U pdating the neighbors distance

9: end for

10: for all v = 0 → sizeOf V olume do

11: F ind out which vertex/voxel has the minimum distance

12: M ark that voxel as true, evolve the algorithm

13: end for

14: while true do

15: U pdating the neighbor0s distance

16: F ind out which vertex/voxel has the minimum distance

17: if nothing has been changed to initial values then

18: break

19: end if

20: M ark that voxel as true, evolve the algorithm

21: end while

22 Implementation

3.2.2

Version 2

Algorithm 2 Dijkstra algorithm O(m + n log n)

1: _{procedure Dijkstra(seeds, distancemap, path, visited, queue)} 2: for all v = 0 → sizeOf V olume do

3: initializing distancemap, path and visited arrays 4: end for

5: for all n = 0 → numberOf Seeds do

6: distancemap(seeds(n)) ← 0

7: visited(seeds(n)) ← true

8: U pdating the neighbors distance

9: Inserting the voxels index and distance to queue which sorts according to distance

10: end for

11: for all v = 0 → sizeOf V olume do

12: P op out the vertex/voxel with minimum distance f rom the queue

13: M ark that voxel as true, evolve the algorithm

14: end for

15: while true do

16: U pdating the neighbor0s distance

17: Inserting the voxels index and distance to queue which sorts

according to distance

18: P op out the vertex/voxel with minimum distance f rom the queue

19: if not finished then

20: do nothing

21: else

22: pop out f rom queue

23: end if

24: if nothing has been changed to initial values then

25: break

26: end if

27: M ark that voxel as true, evolve the algorithm

28: end while

3.3 An application for abdominal organ segmentation 23

3.3

An application for abdominal organ

segmen-tation

An image processing pipeline is built in MeVisLab. The first module is the source data module where one selects the data set one wants to work with. The second module is a marker-to-mask module which makes is possible to mark the images with the seed points. These two modules are sent as input to the Image Foresting Transform module which calculates the accumulated distances (Connectivity Map) and later also calculating the Label Map which is sent as output number two. In the user interface one has used a module named Overlay which means that one can mark different regions from the segmented image and apply that marked area on the original image so that one has the ability to see which areas of the original image are segmented. Input to this Overlay module is not only the segmented data (Label Map) but also the original. Another module used in the user inter-face is the CalculateVolume module which calculates the volume from a specific marked area in the segmented image. Input to that is only the segmented image (Label Map). For the 3D rendering part of the user interface a module called SoGVRVolumeRenderer is used. This allows the user to visualize the segmented parts with Volume Rendering technique (VRT) and also defining an appropriate transfer function.

The user interface is built in MeVisLab using MDL (MeVisLab Definition Lan-guage) and Javascript as complement which sends the variables between modules correctly. The user interface is running the Dijkstra algorithm by default.

Clicking the tab 3D lets the user view the segmented volume with volume ren-dering. Here, the user is allowed to design a transfer function which handles the opacities and colours to do 3D visualizations of the data.

The user interface is built with the purpose of giving radiologists a program to train with and the possibility to give a quantitative measurement in form of the volume of the organ or the lesion of the segmented data. The interface has two windows, the seed planting area to the left and the marked segmented area viewer to the right. The seed planting area is where the user places the seeds. The seed type depends on which organ or lesion the user chooses under the seed planting window, see figure 3.3. Being consistent in the choice of organs and lesions is the most important when using the user interface. As long as the user places the seeds correctly after clicking the specific organ or lesion the user can choose the organ or lesion of interest in the drop-down menu to the right of the marked segmented area viewer. The specific organ or lesion that is marked should now appear as read in the marked segmented area viewer and the red volume is calculated. If one is not satisfied with the calculation one could simply add or remove seed points and then running the computation again by clicking start. Figure 3.4 is showing the 3D part of the software. A manual describing how the user interface works and how to handle it is found in Appendix A.

24 Implementation

Figure 3.3: The 2D interface part

Chapter 4

Results

All experiments were run on a Mac Pro with a 2 x 2.26 Ghz Quad-Core Intel Xeon processor and 8 GB of memory.

4.1

Computation performance

The computation results are comparisons between the Bellman-Ford and Dijkstra algorithm. The O(n2) implementation of the Dijkstra algorithm was not involved in this comparison due to the long computation time of that version. Three exper-iments were done by choosing an appropriate data set from the abdomen, placing the seed points over the volumes in order to get a segmented volume and also test segmenting with many less seed points. When running the largest data set (Abdomen1) one got the biggest difference in computation time between Dijkstra and Bellman-Ford. The Dijkstra run in 3.66 minutes compared to 5 minutes in Bellman-Ford using the same set of seed points in the volume (Table 4.1). Run-ning the algorithms with smaller volumes and comparing the difference in runRun-ning time is not equally remarkable. The Abdomen2 data set run with a total number of 178 seed points in 1.05 minutes for Dijkstra and 0.83 minutes with Bellman-Ford. Segmenting the Abdomen3 data set yielded a running time of 2.6 minutes with Dijkstra and 2.5 minutes with Bellman-Ford. The total amount of voxels varied between ≈ 6 millions for Abdomen2 to ≈ 23.6 millions for Abdomen1. The seeds when segmenting Abdomen3 was 131 and this was also a MRI data set in contrast to the other data sets which were obtained from CT. Table 4.1 shows information regarding running time, seed planting, size and aquisition techniques used when segmenting volumes with these minimum-cost path search algorithms. The resolution (res in Table 4.1) in the data sets is 16 bits in each data set.

26 Results

Table 4.1: Seeding and computation time

Dataset Abdomen1 Abdomen2 Abdomen3

Size (x*y*z*res): 512*512*90*16 414*343*44*16 448*576*74*16

Num seed voxels: 238 178 131

Time Dijkstra (min): 3.66 1.05 2.6

Time Bellman-Ford (min): 5 0.83 2.5

Aquisition techique: CT CT MRI

4.1.1

Segmentation results - Dijkstra vs Bellman-Ford

(a) Dijkstra (b) Bellman-Ford

Figure 4.1: Comparison of a segmented slice using different algorithm

Visual results showing the difference between the both minimum-cost path search algorithms Dijkstra and Bellman-Ford are found in Figure 4.1. Both algorithms run with the same cost function and use the same set of seed points. One can see small differences at some locations but if one considers the overall picture the change is small compared. For a quantitative measurement, a subtraction between the Dijkstra result and the Bellman-Ford results was done. The total number of voxels difference between segmenting with the both algorithms is 23733 voxels. This is about 0.1 percent of the whole volume. Table 4.2 shows this subtraction statistics.

Table 4.2: Results when subtracting Dijkstra from Bellman-Ford segmented data

Volume Number of Voxels

Difference

Total Amount

Voxels

Percent

4.1 Computation performance 27

4.1.2

Cost functions evaluation - Kidney

In this section, results of doing segmentation and calculation of kidney in a data set when using three different cost functions are compared visually and in quantita-tive measurements. Every result is computed with the exact same number of seed points and same location of the seed points. Figure 4.2 shows visually differences when choosing different cost functions. The seed points have been placed in a way that it should segment the kidneys. Number of seeds points is 181 and Dijkstra algorithm has produced everything in this section. Seeding took approximately about 1 minute to finish. The difference is visually distinguishable. One can see differences in the pictures, while using the cost functions 2.4 and 2.6 which has areas above the right kidney which has the same label as the kidney. The result using cost function 2.5 do not have this area here.

(a) Absolute difference (Eq 2.4)

(b) Exponential of absolute difference (Eq 2.5)

(c) Difference squared Eq (2.6)

Figure 4.2: Comparison of segmented slice using different cost functions

Table 4.3 summrizes details of the quantitative comparison. Using the cost func-tion Absolute difference (eq 2.4) a volume of 0.978 liters was segmented and mea-sured, comparing to the manual segmentation result which was 0.597 liter. Using cost function Difference squared (eq 2.6) one yield approximately 100 ml difference in the volume measurements. Finally using the cost function Exponential of the absolute difference (eq 2.5) which is closer to the manually segmented one. 0.744 liters was measured compared to 0.597 in the manual segmentation.

Table 4.3: Results when calculating volume of kidney with different cost functions

Organ Cost

func-tion

Voxels Total

vox-els

Percent Volume(l)

Kidney Abs diff 261512 23592960 1.1 0.978

Kidney Abs diff

squared

233424 23592960 0.989 0.873

28 Results

Table 4.4: Volume of kidneys manual segmented

Organ Time(min) Volume(l)

Kidney 8 0.597

When first calculating the kidney volume it was segmented with the three different cost functions described in section 2.5. In all three cases, 181 seed points was used at the same location. The seeding took about one minute to finish and the compu-tation time was 3.8 minutes. The total time with seeds planting and compucompu-tation time is thus 4.8 minutes. Manually marking the kidneys in the original data set took 8 minutes. This was done manually in each third slice throughout the kidney marking the kidneys without its vessels or anything else. After this marking of this slices one interpolates the rest of the slices in order to save time. The volume is fairly good measured when looking at the marked area in each slice.

4.1.3

Cost functions evaluation - Liver

The liver seems to be a more difficult organ to segment due to the closeness of the surrounding organs and the similar intensities in those regions. The size of the liver makes it difficult to choose appropriate seed points in the area. Regarding manually segmenting this organ this is very time consuming due to the size. This size of this data set was 512x512x81x16. 65 seed points have been used. While segmenting the liver it is hard to get a good result after the first computation. One often needs to modify the seed points - add more and then check the new results if some improvements were made. In this experiment one needed to do three seed plantings in order to get an acceptable segmentation of the liver. This concern every cost function used. After these three seed plantings, which took a total time of 4 min, fairly good results were achieved. The computation time for every cost function was 3x6 min, about 18 min. The total time including seed planting and computation time was about 22 min for the liver. Notice the computation time of this data set concerning how many seed points is used, regardless of which cost function was used. Before, in table 4.1 we recieved a result where a data set of similar size took approximately 3.5 minutes to segment. Now the data set is smaller, but with about four times less seed points, the computation takes longer time.

4.1 Computation performance 29

(a) Absolute difference (Eq 2.4)(b) Exponential of absolute dif-ference (Eq 2.5)

(c) Difference squared (Eq 2.6)

Figure 4.3: Comparison of segmented slice of liver using different cost functions

Table 4.5: Results when calculating volume of kidney with different cost functions

Organ Cost

func-tion

Voxels Total

vox-els

Percent Volume(l)

Liver Abs diff 556775 21495808 2.59 1.574

Liver Abs diff

squared

585100 21495808 2.71 1.654

Liver Exp diff 572504 21495808 2.66 1.618

Table 4.6: Volume of liver manual segmented

Organ Time(min) Volume(l)

Chapter 5

Discussion

5.1

Computation results

While testing the performance of the Dijkstra and Bellman-Ford algorithms, in-teresting results were obtained. For a smaller data set with less seed points, the Bellman-Ford was faster than Dijkstra. When testing with a larger data set with even more seed points one can notice a significant difference in the performance. A reason for this is that the Bellman-Ford algorithm is scanning the volume back-wards and forback-wards until no change or update is being done in the volume. This means that, if the paths calculated in the image are complicated and take long time, the Bellman-Ford algorithm has to do more scans through the volume, which can lead to a longer computation time. The Dijkstra algorithm works in some-what greedy way. This means that it will traverse all the voxels in the volume anyway and does not depend on how complicated the paths are. The computations are also highly affected by the location of the seed points and how many they are. If one should segment a whole volume of millions voxel from for example three seed points will be difficult for both algorithms. Here, long paths arise, making the computation time for Bellman-Ford longer. Dijkstra will also be affected by this, but not to the same extent.

32 Discussion

5.2

Visual results of Dijkstra and Bellman-Ford

Since the Dijkstra and the Bellman-Ford algorithm solve the same problem, one could expect the same result. Figure 4.1 shows that the results are roughly the same visually. Table 4.2 shows that the difference in voxels is only 0.1 percent when subtracting the volumes of the algorithms when choosing the same cost function. Thus, for practical purposes the both algorithms give the same result both re-garding visual interpretation and when subtracting the voxels in the segmented volumes. This is interesting because they solve the same problem in a different manner, one could expect that the result should be the approximately the same.

5.3

Cost functions and volume calculation

Regarding the different cost functions, a difference in segmentation was found. In this case the two first cost functions 2.4 and 2.6 has some artifacts as one can see in the figure 4.2. There are some areas over the right kidney which has the same label as the kidney itself. Segmenting with the cost function 2.5 did not give this result. Neither does it tend to include the vessels from the kidneys as much as the other two seem to do. That is two reasons why the volume might be overestimated with those cost functions. In the end there was about 150 milliliters difference between the manually marked kidneys and the best approximated volume. The kidneys are relatively easy to segment and one can achieve fairly good results just after one computation.

A possible reason why every cost function calculated a larger volume than the manually segmented volume is that the penalty at edges is not big enough. The cost functions tend to leak across borders and include vessels and other surround-ing organs and tissues when dosurround-ing the calculation. The difference in time was 4.8 min for automated segmentation compared to 8 min for manual marking. One should remember, though, that the placing of seed points is of great importance, and it is possible to get better results with another cost function when placing the seeds better, but that has to do with experience in the segmentation.

The liver is a much harder organ to segment. The closeness of adjacent organs and similar intensities in those gives rise to problems when the algorithm is trying to segment it. The segmentation has greatest problem with the fat and muscles around the liver, but also the heart is quite difficult because of the very smooth transition in intensities. It is almost impossible to segment the liver correctly from just one seed planting. In this example one used three seed plantings and every computation took about six minutes. The total time was about 22 minutes com-pared to 12 minutes when manually marking the organ every 3rd slice and then interpolating. However, the manual marking would take extensively more time if the user marked the liver in every slice. In this study, we run the algorithm with the cost function in equation 2.5 trying to get as good segmentation as possible, and then tested the exact number and placement of the seed points with the other cost functions. According to figure 4.3 the cost function 2.4 has problems with a

5.4 Future work 33

quite homogenous area in the liver. The quantative measurement of the segmented liver volume agreed best with the manual method for cost function 2.4, but when inspecting the visual results one can see that it misses to segment important -crucial and relatively easy parts of the liver. The cost function 2.5 did the best segmentation again. As mentioned earlier the most difficult part, which has to be manually corrected afterwards, is the fat and muscles around the liver.

5.4

Future work

5.4.1

Improving Dijkstra

The performance of the Dijkstra algorithm can be improved further. In [19] sev-eral methods are described to reduce the running time. In the paper structures as Van Emde Boas priority queues, Fibonacci heaps, radix heaps are discussed in detail. However, an implementation of this kind would probably need more time and is out the scope of this thesis. In [19] a bound is achieved with an implemen-tation of a combination of radix heaps and Fibonacci heaps where the bound reads

O(m + nplog(C)). Here, m is the total amount of edges (26*amount of voxels

in 3d), n is the total amount of voxels and C a constant. One should remember, though, when implementing this kind of optimizations that what the theory says does not mean that the algorithm with this solution will be. There are more things to consider when it comes to optimizations, but the choice of algorithm is an im-portant part. One has to test and evaluate these optimizations before knowing how it affects the actual problem.

Another possibility of improving the computation time of the Dijkstra algorithm is by implementing the algorithm on the GPU. A successful approach was done by Kauffmann and Piche in [20]. Due to the good programmability and capability in today’s GPU’s a wide range of applications can be implemented on the GPU. However, due to time limit, this project has focused on implementing it on the CPU.

5.4.2

Improving user interface

In this work, the focus is on the development of the Dijkstra algorithm and com-paring it with Bellman-Ford. The comparisons between the different cost functions is also a big part. The user interface was built in a way including the basic features in order to help a radiologist to segment organs and lesions. This interface can be improved in many ways. Now, the user adds seed points manually by marking the organs. This could be improved by implementing a tool such as a brush that lets the user draw areas in the images. Also a tool like the magic wand would be an improvement. There, the user can use the mouse and click on a specific object and the software will mark the object in three dimensions giving the user quantitive information and also visual where it marks everything in the picture. This makes

34 Discussion

the interface easy to interpret even for less experienced users. The interface now lets the user choose organ or lesion from a drop down menu. This may take longer time in interaction with the program, but due to time limitations of this work, this solution was chosen.

5.4 Future work 35

5.4.3

Cost functions

Further evaluation could include testing other cost functions that might yield different results. It could involve cost functions based on gradients. Another approach is using a tensor cost function that would be interesting to give a more general cost. Then there is a similar function such as Udupa’s affinity function [4], which has quite another approach. In this project three different cases were evaluated. One needs to do a more comprehensive study in order to decide if this can be used in clinic. Also, the results must be validated. Still, the two cases show the ability of the image foresting transform when changing the form of the cost of the path. When using this segmentation tool one should also recall that for using the tool the results is improved with increasing experience of how it works. The goal is to let the user place as few seeds as possible but still segmenting the image well in the specific area. This is difficult and requires one to test and evaluate many types of cost functions. One cost function may work out to be good for segmenting a particular organ or lesion, another cost function work for segmenting other. This needs to be further evaluated and include more statistical measurements.

5.4.4

Conclusion

This study showed that the image foresting transform gave similar segmentation results when using Dijkstra’s algorithm and Bellman-Ford’s algorithm. The time difference between the two algorithms was not that remarkable. However, the segmentation results varied considerably when different cost functions based on intensity difference were chosen.