Segmentation of Carotid Arteries from 3D and 4D Ultrasound Images

(1)

Segmentation of Carotid Arteries

from 3D and 4D Ultrasound

Images

Per Mattsson

Andreas Eriksson

LiTH-ISY-EX-3279-2002 Department of Electrical Engineering

Linköpings Universitet, SE-581 83, Linköping, Sweden

(2)

Printed by:

UniTryck, Linköping, Sweden Distributed by:

Linköpings Universitet

Department of Electrical Engineering SE-581 83, Linköping, Sweden

No part of this publication may be reproduced, stored in a retrieval system, or be transmitted, in any form or by any means, electronic, mechanic, photocopying, recording, or otherwise, without prior permission of the authers.

(3)

iii

Abstract

This thesis presents a 3D semi-automatic segmentation technique for extracting the lumen surface of the Carotid arteries including the bifurcation from 3D and 4D ultrasound images.

Ultrasound images are inherently noisy. Therefore, to aid the inspection of the acquired data, an adaptive edge preserving filtering technique is used to reduce the general high noise level. The segmentation process starts with edge detection with a recursive and separable 3D Monga-Deriche-Canny operator. To reduce the computation time needed for the segmentation process, a seeded region growing technique is used to make an initial model of the artery. The final segmentation is based on the inflatable balloon model, which deforms the initial model to fit the ultrasound data. The balloon model is implemented with the finite element method.

The segmentation technique produces 3D models that are intended as pre-planning tools for surgeons. The results from a healthy person are satisfactory and the results from a patient with stenosis seem rather promising. A novel 4D model of wall motion of the Carotid vessels has also been obtained. From this model, 3D compliance measures can easily be obtained.

Keywords: Stroke, Carotid arteries, 3D Ultrasound, Adaptive Filtering, Segmentation, Edge Detection, Monga-Deriche-Canny, Seeded Region Growing, Inflatable Balloon Model, Finite Element Method, Triangulation, 4D wall motion, 3D compliance

(4)

(5)

v

Sammanfattning

Detta arbete presenterar en 3D halv-automatisk segmenteringsteknik för att extrahera halsartärens inre kärlväggen samt dess delning från 3D och 4D ultraljudsbilder.

Ultraljudsbilder är av naturen väldigt brusiga. För att förbättra den visuella tolkningen av insamlade data används en adaptiv kantbevarande filtreringsteknik som sänker den generellt höga brusnivån. Den egentliga segmenteringsprocessen börjar med en kantdetektering med en rekursiv och separabel 3D Monga-Deriche-Canny operator. För att minska tidsåtgången för segmenteringsprocessen används en teknik med frösåddbaserad regionstillväxt, vars resultat utgör en initialmodell av artären. Den slutliga segmenteringen av artären är baserad på ballongmodellen som deformerar och expanderar initialmodellen till att bli optimalt anpassad till de kanter som kan upptäckas i insamlade ultraljudsdata. Ballongmodellen är implementerad med finita element metoden.

Segmenteringstekniken skapar 3D modeller av artärerna som ytterst är avsedda som planeringsverktyg för kirurger. Resultaten från en frisk person förefaller helt tillfredställande medan resultaten från en patient med stenos kan sägas vara lovande. Dessutom har för första gången en 4D modell av väggrörelser hos halsartären tagits fram. Från denna modell kan sedan 3D modeller och kvantifiering av flexibilitet erhållas.

(6)

(7)

vii

© 2002 Per Mattsson and Andreas Eriksson. 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 owners. The publisher has taken technical and administrative measures to assure authenticity, security and accessibility. According to intellectual property law the authors has the right to be mentioned when their 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/.

(8)

(9)

ix

Preface

This master’s thesis is the last requirement for our Master of Science Degree from the Linköping Institute of Technology. It was written during the fall of 2001, at Karolinska Institutet, Division of Medical Engineering in Stockhom in collaboration with Huddinge University Hospital in Stockholm, Sweden and VingMed Sound AS, Norway. Professor Per-Erik Danielsson, Image Processing Lab, Deptartment of Electrical Engineering at Linköping University has been our examiner and the project was supported by MD, PhD and Associate professor Lars-Åke Brodin working at Huddinge University Hospital.

Acknowledgements

This work had not been possible without the help from the following people that we would like to thank.

First of all we would like to thank Professor Per-Erik Danielsson, not only for being the examiner of our project, but also for all the help he has given us and the interest he has shown. He introduced us to the interesting field of image processing and gave us the burning interest we have today. His knowledge and ideas has been invaluable to us.

Thanks to Professor Håkan Elmqvist at the Department of Medical Engineering, Royal Institute of Technology in Stockholm for letting us do our work at his department and for all the support he has given us.

(10)

Preface x

Thanks to our supervisor Lars-Åke Brodin. He came up with the idea for the project and helped us with the contact with VingMed Sound AS.

We also wish to thank MSc Qingfen Lin, MSc Katarina Flood and PhD Maria Magnusson-Seger at the Image Processing Laboratory for their support and enthusiasm in the beginning of our project. We are very grateful for your ideas and comments in the initial phase of the project.

Thanks to Licentiate in Technology Björn Johansson at the Computer Vision Laboratory at Linköping Institute of Technology for explaining adaptive filtering to us.

Thanks to Britta Lind and Rita Balsano for helping us with the clinical part of our work. We spent many hours together performing ultrasound examinations and they taught us a lot about how to interpretate the ultrasound images of the Carotid arteries.

Thanks to Tomas Jogestrand, clinical director of the Physiology Department at Huddinge University Hospital, and Jacek Nowak, Associate professor, for showing a great interest in our project and for helping us with the layout of the interface. They also helped during the image acquisition.

Thanks to Tim McInerney, Assistant Professor in the School of Computer Science Department of Math, Physics and Computer Science at Ryerson Polytechnic University and an Adjunct Professor in the Department of Computer Science, University of Toronto, for helping us with our questions concerning the balloon model and providing us with good references to literature.

Thanks to Wim Lenferink at numerical analysis, department of mathematics at Linköping Institute of Technology for helping us with the numerical aspects of the finite element method.

We are also most grateful to Professor Anders Eriksson, from the Structural Mechanics Research Group at the Royal Institute of Technology for his incredible enthusiasm for our project and for all the help he gave us concerning the finite element method. We spent several hours in his office

(11)

Acknowledgements xi

asking him questions about this method and he always took the time to answer them.

Of course, the time here at the Department of Medical Technology would not have been as enjoyable as it was without all our friends here, you all know who you are. Special thanks to MD, Phd student Magnus Fredby, for interesting discussions and for all ideas and comments about our project, and Phd student Camilla Storaa, for helping us with all the computer problems.

It has been very interesting working with this project and we are very grateful gaining so many new contacts and friends. Now, four and a half years of studying are over. At the time of writing the report the following addresses could be used to contact us:

Per Mattsson

Andreas Eriksson

(12)

(13)

xiii

2.3.1 Harmonic imaging ... 21 2.3.2 Compound scanning ... 22 2.4 Image acquisition... 23 3 Adaptive Filtering 27 3.1 Nagao-Ye filtering ... 28 3.2 Granlund-Knutsson filtering ... 30 3.3 Results ... 34 4 Edge detection 37 5 Approximation with seeded region growing 47 6 The balloon model and the finite element method 51 6.1 Overview ... 51

(14)

Contents xiv

6.2 The inflating balloon model... 52

6.3 The finite element method ... 56

6.3.1 Discretization of the structure... 59

6.3.2 Approximation with shape functions... 60

6.3.3 Analysis of thin plates ... 62

6.3.3.1 The triangular membrane element... 63

6.3.3.2 The triangular bending element... 68

6.3.3.3 The shell element – Combining the membrane and bending element ... 73

6.3.4 Assembling the total stiffness matrix ... 76

6.3.5 Computing nodal forces ... 77

6.3.6 Dynamic finite element representation... 80

6.3.7 Nonlinear FEM... 82

7 Implementation 85 7.1 The user interface ... 85

7.2 Implementation of the program ... 87

7.2.1 Filtering module ... 89

7.2.2 Edge detection module ... 89

7.2.3 Approximation module... 89

7.2.4 Force calculations module ... 89

7.2.5 Triangulation module ... 91

7.2.6 FEM module... 94

7.2.7 User interface ... 95

8 Experimental results 97 8.1 Overview of the segmentation process ... 97

8.2 A healthy person ... 101

8.3 A patient with stenosis... 106

8.4 Execution times ... 108

9 Concluding remarks and future directions 111 Appendix A Estimation of orientation 115 A.1 Tensors ... 119

A.2 The orientation tensor ... 124

A.3 Orientation estimation using lognormal filters ... 125

A.4 Interpretation of the 3D orientation tensor... 130

Appendix B Control of filter shape 135

(15)

xv

Appendix D Shape functions 143

Appendix E Deriving the mass and damping matrices 145

Appendix F The Runge-Kutta method 149

Index 151 Bibliography 157

(16)

(17)

xvii

List of Figures

1.1 Structure of the artery ... 3

1.2 The left and right Carotid arteries...4

1.3 An angiography image of a Carotid stenosis ... 5

1.4 The surgical procedure of removing plaque from the Carotid artery. ... 5

1.5 Velocity profile for a Newtonian fluid...6

1.6 Velocity profile of a Newtonian fluid. ...7

1.7 Velocity profile of a Casson fluid. ...7

1.8 Velocity and pressure profile through a stenosis ...9

1.9 Places that are prone to be exposed by plaque...9

1.10 The procedure to segment a 3D model from an ultrasound examination... 12

2.1 Reflection and refraction of ultrasound at an interface between muscle and bone. ... 17

2.2 False echoes because of reverberations. ...19

2.3 False echo because of strongly reflecting surface not perpendicular to the transducer... 20

2.4 Spatial distortion due to materials with different speed of sound. ...21

2.5 Comparison between conventional and harmonic imaging. ...22

2.6 Comparison between conventional image and compound image ...23

2.7 During the merging, the frames are assumed to be parallel. But with freehand acquisition the frames become tilted reducing the quality of the 3D volume. This is corrected for if a position detecting devise is used... 24

2.8 The principle behind time-gating. A number of volumes are created, each consisting of frames taken from the same position in the heart cycle. The time-gated volumes are not as dense but do not contain motion artefacts... 25

2.9 Comparison between 3D- and 4D-acquisition. Since the 4D volume only consists of 15 frames the level of detail in one direction is very low, causing the blurred appearance. In 3D-acqusition mode, all the acquired frames are used to produce the 3D-volum. Therefore the level of detail is much higher at the cost of motion artefacts. ... 25

(18)

List of Figures xviii

3.1 The seven subregions used in Nagao-Ye filtering with 6-connectivity... 28 3.2 The two steps of Nagao-Ye filtering. First the mean and variance

is calculated for each 3x3x3 subregion. Then, the mean of the subregion with the lowest variance around each voxel is taken as the output value... 29 3.3 Overview of the Granlund-Knutsson filtering technique. The input

signal is filtered with a set of lognormal quadrature filters with different directions to produce orientation tensors T. These tensors are then lowpass filtered into TLP and mapped into control tensors

C. The control tensors are used to do a weighted summation of the output from a set of shift invariant filters to produce the output signal... 31 3.4 Comparison of Gaussian smoothing and the two adaptive filtering

techniques. (a) shows a slice through the unprocessed volume. (b) is the result after filtering with a 5x5x5 Gaussian kernel. (c) is the result after filtering with the Nagao-Ye method (6-connectivity and one iteration). (d) is the result after Granlund-Knutsson filtering ... 35 4.1 The isotropic rotational-invariant Laplacian operator...40 4.2 The Laplacian results in two spikes when detecting a line. One

spike is negative and one is positive. Therefore, to detect the line, the Laplacian operator looks for places in the image with zero-crossings ... 41 4.3 The two images (b) and (c) are results of the thresholding from (a).

The edges are displayed as lines of perfect zero-crossings of the Laplacian operator. (b) has a higher threshold than (c) resulting in less detected edges in (b) than in (c). (b) has large gaps in the contour while (c) has false edges indicated by the red arrow. From these images a hysteresis is made to produce curves that are more closed as in (d). ... 42 4.4 The three slices which will be analysed... 44 4.5 The Carotid Communis. (a) is the adaptive filtered data and (b) is

the result from the edge detection. ... 45 4.6 The bifurcation of the Carotid artery. ...45 4.7 Carotid Interna and Externa. ...45 5.1 A gap in the contour where the region growing might “leak out”

into the neighbouring blood vessel ... 48 5.2 A slice of a volume before and after down-sampling with a factor

of four. ... 48 5.3 Result of the region growing in two slices of a volume... 48

(19)

xix

5.4 The approximation shown as a triangular mesh. Before and after

reduction of the number of triangles. ... 50

6.1 The potential field around the Carotid artery for a slice. ...55

6.2 Flowchart of the finite element method ...58

6.3 Tessellating the model into a set of triangular elements. ...59

6.4 Normal and inplane loads. ...62

6.5 The local coordinate system for the membrane element... 63

6.6 The local coordinate system for the bending plate and its nine degrees of freedom... 68

6.7 Forces and moments in a plate...70

6.8 The relationship between a triangular element in the xy-system and the reference element in the xh-system. ...78

6.9 Sampling points used in quadrature formulas of table 1...80

6.10 Comparison between linear and nonlinear solutions. The structure is a flat plane with two of its edges clamped. The forces are applied at the middle of the right edge and are indicated by the arrows. In (a) the solutions are almost identical and the linear solution only slightly overestimates the deformations. In (b) the force is four times larger then in (a) and the linear solution severely overestimates the deformation and does not capture the forces that pulls the edge inwards to the left since inplane forces starts to appear. ...83

7.1 The program after opening a new examination. By clicking and dragging the mouse in the three orthogonal cut-planes, the user can crop the imported volume ...86

7.2 Result from cropping and filtering the volume. ...86

7.3 Result from the approximation with seeded region growing. ...86

7.4 Resulting model of the Carotid artery. The model is shown in the bottom-right view and as a red contour in each orthogonal view. ...88

7.5 3D-model of the Carotid artery that can be rotated by using the mouse...88

7.6 The modules of the program and the flow of information between them. The functions performed by each module are listed together with the language used for the implementation. ...88

7.7 Flowchart of the program. The inner-loop of the program is the “ODE-solver” which repeatedly calls “Compute nodal force vector” and, if needed, recomputes the triangular mesh and its elements. ...90

7.8 Geometric transformations. A creation or inversion is performed when an edge is too long, while a fusion is performed for edges that are to short. ...91

(20)

List of Figures xx

7.9 The relations between the Mesh, Node, Edge and Triangle objects. ...93 7.10 Each Edge object has references to the Nodes U, V, SU and SV. It

also has references to the two triangles SU and SV...93 8.1 The steps in the segmentation process. (a) is a slice from the

ultrasound volume. (b) is the result after filtering with the Granlund-Knutsson method. (c) is the result after approximation with seeded region growing. (d) is the result from the edge detection with the Monga-Deriche-Canny detector used to produce the potential field displayed in (e) and (f). The white curve in (e) is the cross-section of the approximation from (c). The white curve in (f) displays the final segmentation result after the inflatable balloon model. The final segmentation is shown as a red curve on the original data in (g)...98 8.2 The approximation and the final segmentation...99 8.3 Part of the final segmentation and a slice of the potential field

displayed as a surface. ...99 8.4 The segmentation near a false edge. (a) displays the initial model

for the inflatable balloon model. During the segmentation process, the inflation force acting on the model pushes the model out of the potential valley (b). The final segmentation is shown in (c)...99 8.5 A slice from the patient with stenosis (a). The segmentation is bent

around the plaque like a banana. A part of the neighbouring vein is showing at the top left corner. Notice how the low image quality makes the edge detection uncertain and how the potential field (b) has both many gaps and false edges. In cases like this, it is hard to say what the correct shape is...99 8.6 A 3D model of the Carotid artery and its bifurcation seen from two

different angles. ...101 8.7 3D model consisting of 17322 elements ...103 8.8 4D model consisting of 18724 elements. Since the 4D volume

consists of fewer slices, the level of detail in one direction is very low. ...103 8.9 (a) shows the artery during diastole. This model is colour-coded

with respect to the expansion or compression from the previous systole. (b) shows the artery during systole. ...105 8.10 Wall-motion of a cross-section of Carotid Communis. The vertical

axle in (a) corresponds to time and the horizontal axle corresponds to the measured diameter in (b). A comparison between two cross-sections of the model at systole and diastole are shown in (b). In this image it can be seen that the deformation of the artery is mainly in one direction. ...105

(21)

xxi

8.11 3D model consisting of 17228 elements, of a Carotid artery with

stenosis...107

8.12 A 3D model segmented from the first time step of a 4D volume from a patient with stenosis. The model contains 19424 elements. ...107

A.1 Illustration of the principle behind adaptive filtering using steerable filters. (a) is a black and white picture with added white noise. Filtering with a square averaging kernel gives the picture in (b). The noise is reduced but the transition between the black and the white area is not as sharp as in (a). In (c) a narrow averaging kernel parallel to the edge has been used. The noise is suppressed as in (b) while the transition is as sharp as in (a). ...115

A.2 The difference between direction and orientation...117

A.3 Examples of simple neighbourhoods with only one orientation ...118

A.4 Representing orientation with double angle vectors. ...118

A.5 A lognormal filter function. ...128

A.6 (b) is the local 2D DFT of the area inside the square in (a). (c) shows a vertical and horizontal slice of the 2D DFT of an horizontal edge such as the edge in (a). 40 different areas have been used for the estimation. The blue curve in (d) is an estimation of the SNR of an edge based on a comparison between the vertical and horizontal slice in (c) The red curve is the frequency characteristic of the lognormal filter chosen for orientation estimation...129

A.7 Planar structure. The Fourier transform is concentrated along a line given by the first eigenvector. This eigenvector is normal to the planar structure. ...130

A.8 Line structure. The Fourier transform is concentrated to a plane spanned by the two first eigenvectors. The third eigenvector is oriented as the line in the spatial domain. ...131

A.9 Isotropic case. No dominant orientation is present and all eigenvalues are equal. ...132

A.10 A slice through a 3D-volume (a) and the eigenvectors (b). The eigenvectors indicate the dominant orientation. The greyscale levels of the arrows represent the size of the eigenvalue with dark arrows correspond to large eigenvalues. The image in (c) displays the measures P1, P2 and P3 (plane case, line case and isotropic case) as red, green and blue respectively. (d) is the same as (c) but after the tensor lowpass filtering...134

B.1 Plot of the m-function. a is the value of r for which the function equals 0,5. b determines the slope of the function...137

(22)

List of Figures xxii

B.2 m T

(

_lp ,s ,d ,b

)

as a function of T_lp

d b

for four values of s : 0.05 (topmost curve), 0.1, 0.15, 0.2. =0, =2 ...139 B.3 m T

(

_lp,s ,d ,b

)

as a function of T_lp

b

for five values of d: 0 (bottom curve), 0.25, 0.5, 0.75, 1. =0.05, =2 ...139 s

(23)

xxiii

Notations

Below follows a list of notations used in this thesis. Less frequently used symbols, and symbols that have different meaning in different contexts, are defined where they are used.

The following notations are used for mathematical entities:

· t Scalars (lowercase letter in italics) · t Vectors (lowercase letter in boldface) · T Matrices (uppercase letters in boldface) · t(x) Functions (lowercase letter)

The following notations are used for mathematical operations:

· tˆ Normalized vector

T

A

· Matrix and vector transpose · t Absolute value of real numbers · t Matrix or vector norm

b a×

· Pointwise multiplication of elements ·

(

g*f

)( )

x Convolution

q& g Ñ

· The time derivative of q

(24)

(25)

1

1 Introduction and Background

1.1 Motivation

Stroke has a great social and economic impact on modern society and is now the subject of substantial research efforts. It is the third leading cause of death in the western world and the primary cause of adult disability. There is therefore an urgent need for better techniques to find patients at risk of stroke, examine them and deliver guidelines for the choice of treatment. Finding new imaging methods and improving the detection capabilities of existing ones might improve the life quality of many patients.

The Carotid arteries are the major vessels that supply the head with blood and are quite often subject to build-up of plaque. Pieces of the plaque may break loose and follow the blood flow into the capillaries of the brain and cause blockage. This kind of stroke accounts for approximately 75 percent of all strokes.

(26)

2 Chapter 1 Introduction and Background

· Digital subtraction angiography is the present standard imaging technique for monitoring the Carotid arteries. This is an invasive technique and therefore poses a certain risk to the patient. · 2D Doppler ultrasound is a cheap and non-invasive technique

that is becoming a useful tool in assessing disease severity of the Carotid arteries. In the presence of plaque, changes in the blood flow through the arteries can be seen. However, it does not provide a sufficient view of the surface of the plaque and it also suffers from high intra- and inter-operator variability.

Today patients having stroke symptoms are first examined with ultrasound for a first diagnose to select the patients who would benefit from operation. These patients are then sent to the radiology department to be examined with angiography, which is the current standard method for producing images that are used when planning the operation. Instead, if a trustworthy model of the Carotid artery could be obtained from of the ultrasound examination directly, the patients do not have to undergo the angiography examination, saving time and money for the hospital and a considerable risk for the patient.

Consequently, it should be worthwhile to investigate the possibilities to use truly three-dimensional ultrasound imaging and image analysis for the task of assessing disease severity of the Carotid artery [1].

1.2 Stroke and the Carotid system

Blood vessels form a closed system of tubes that carries blood away from the heart, transports it to the tissues of the body and than returns it to the heart. The

(27)

1.2 Stroke and the Carotid system 3

wall of an artery has three coats, with the innermost coat being the endothelium (Figure 1.1). The endothelium is a continuous layer of cells that line the inner surface of the entire cardiovascular system. Normally, the only tissue that blood comes in contact with is the endothelium. The hollow center through which blood flows is called the lumen.

Figure 1.1 Structure of the artery [2].

The Carotid system is located in the neck and contains the Carotid

communis, which branches into Carotid interna and Carotid externa (Figure 1.2). The left Carotid communis artery is the second branch from the arch of the aorta. It divides into Carotid interna and externa the same way as the right Carotid communis. The Carotid interna artery supplies blood to structures inside the skull like most of the cerebrum of the brain. It also supplies blood to the eyeballs, ears, and external nose. The general distribution of the Carotid externa artery is to structures external to the skull.

(28)

4 Chapter 1 Introduction and Background Carotis Interna Carotis Externa Carotis Communis Carotis Interna Carotis Externa Carotis Communis

Figure 1.2 The left and right Carotid arteries [3].

A stroke occurs when the blood supply to parts of the brain is suddenly interrupted, ischemic stroke, or when a blood vessel in the brain bursts, spilling blood into the spaces surrounding brain cells, hemorrhagic stroke. For all types of medical intervention, treatment must be given immediately, as neuronal death progresses quickly after the onset of symptom. The most common form of stroke is ischemic stroke, which occur when an artery becomes blocked, suddenly decreasing or stopping the blood flow. Blood clots are the most common cause of artery blockage and brain infarction. Ischemic strokes can also be caused by stenosis. This kind of stroke accounts for approximately 75 percent of all strokes [4].

Carotid stenosis refers to the blockage and narrowing of the Carotid artery. This blockage is caused by fatty build-up called plaque and is also referred to as arteriosclerosis. This fatty material accumulates in the inner lining of blood vessels and results in narrowing and irregularity of the artery as seen in Figure 1.3. The complex structure of the plaque cause heterogeneous stress-strain states, which may ultimately lead to rupture and the formation of blood clots that can flow up to the brain and cause ischemic stroke.

(29)

1.2 Stroke and the Carotid system 5

Figure 1.3 An angiography image of a Carotid stenosis. The line indicates where the lumen wall is supposed to be.

The decision to treat narrowing of the Carotid artery is not always straightforward. The potential benefit of the surgery must be weighed against the risks of the surgery. The degree of stenosis of the Carotid artery and the presence or absence of symptoms are some of the important factors to consider when making this decision. The degree of stenosis is measured as the difference between the largest and smallest area of the artery in relation to the largest area.

Carotid artery An incision is made to open the Carotid artery. Plaque is removed. Then the repaired artery is closed. Carotid artery An incision is made to open the Carotid artery. Plaque is removed. Then the repaired artery is closed.

Figure 1.4 The surgical procedure of removing plaque from the Carotid artery.

Compared to medical therapy alone, surgery has been found highly beneficial for persons who have already had a stroke or experienced the warning

(30)

signs of a stroke and have a severe degree of stenosis of 70-99%. Figure 1.4 show the procedure of removing plaque from the Carotid artery. For a degree of stenosis of less then 30%, medical therapy is preferred. For a degree of stenosis between 30 and 70%, the best therapy has not been determined, since the risk/benefit ratio varies between the conditions of the patients. Patients that are at a high risk for a surgical procedure may be placed on medications to inhibit their blood from clotting.

1.3 Fluid mechanics and elasticity in vessels

Blood consists of both plasma and blood cells and is therefore not a homogeneous fluid. Pure plasma is a Newtonian fluid, which means that the shear stress is proportional to the velocity gradient of the fluid [5].

0 y ) ( y u dy du 0 y ) ( y u dy du

Figure 1.5 Velocity profile for a Newtonian fluid.

(31)

1.3 Fluid mechanics and elasticity in vessels 7

dy du

m

t = (1.1)

where µ is the coefficient of viscosity. Blood is a non-newtonian fluid and do not in general follow equation (1.1). But in larger vessels like the Carotid system, where the diameter of the vessel is much larger than the diameter of a blood cell, the blood can be modelled as a homogeneous fluid, which is an important characteristic of the blood. Its viscosity reduces as the shear stress increases and therefore, the blood can be considered as a Casson fluid. This means that it follows Casson’s viscosity model, which is more accurate than the Newtonian model for blood flow at low shear rates.

At relatively low flow velocities in large vessels, the flow can be considered to be steady and laminar. The flow is called laminar when it flows in parallel layers that don’t mix. The opposite of laminar flow is turbulent flow, which is swirling. For Newtonian fluids, the velocity profile for laminar flow in a tube is parabolic as seen in Figure 1.6. The velocity is zero at the wall and reaches a maximum at the centreline. The velocity profile for laminar flow of a Casson fluid is rather similar as seen in Figure 1.7.

(32)

For Casson fluids, the core flow in the centre of the tube moves as a rigid body. The flow rate is proportional to the radius in fourth power and the resistance is thereby inversely proportional to the radius in fourth power, which is why small alteration in the radius gives a large effect on the resistance of the flow. In a small vessel, a larger portion of the blood is near the wall, where it is slowed down by frictional forces and therefore the Casson fluid model is necessary. In bigger vessels the frictional forces interfere less and the flow velocity is higher and therefore the simpler Newtonian model might be sufficient.

1.3.1 Flow through stenosis and bifurcations

An increased amount of energy is required to pump the blood through the narrowing at a stenosis. The constriction results in an abnormally high pressure gradient across the opening.

Figure 1.8 is a model of the actual flow thorough a stenosis. To the left is a nozzle and it should be noticed that the nozzle flow has a favourable gradient. The nozzle flow never separates from the endothelium, nor does the throat flow. But in the expanding area an adverted gradient gives the flow low velocity and increased pressure. If the expanding area angle is too large, it will even result in backflow as to the far right in Figure 1.8. Energy is dissipated by viscous forces and by turbulence. The process of conversion of potential energy to kinetic energy in the throat is resulting in an increased velocity. Most of the energy losses occur as the jet expands again to fill the downstream as the flow decelerates which will heat the tissue. This can eventually result in additional damage on the lumen.

(33)

Figure 1.8 Velocity and pressure profile through a stenosis [6].

A bifurcation of a vessel will also change the velocity profile of the flow. In the bifurcation, the nearly parabolic profile will be separated into two. Each of the two branches will have a flow profile, which is wry against the innermost walls of the bifurcation. The resistance along the walls will eventually reduce the velocity of the flow near the wall, ones again resulting in a Casson velocity profile. The same modelling can be used at a bend of a vessel. These places, which have high affliction and thereby more damage, are prone to be exposed by plaque (Figure 1.9).

(34)

1.3.2 Elasticity in vessels

It is the gradient of pressure that is the force that makes the blood flow through the vessels. The arterial system can be considered as an elastic chamber, which represses the blood volume in systole and delivers a suppressed, almost constant flow to the capillaries in diastole. The circulation system is a network of elastic tubes with complex geometry and non-linear elastic properties. It is therefore a very complicated system to model. For practical reasons it can be sufficient to describe the mechanics of the vessel wall as a function of pressure and diameter.

Many medical conditions, such as arteriosclerosis, change the mechanical properties of the vessel walls, which is why physicians study parameters like

compliance and distensibility [7]. The build-up of plaque makes the wall stiffer

which decreases the compliance. Today, measurements of vessel compliance are two-dimensional.

It is believed that plaque rupture is due to functional impairment of the vessel caused by physical damage and that mechanical features of the arteriosclerotic lesion significantly influence its probability of rupture. Recognition of specific features that increase vulnerability to rupture would allow prediction of plaque that are most likely to cause acute events.

Locally, the elasticity of vessels can be characterised by the parameter Cl, local arterial compliance, which is defined as the ratio between the change in area and the change in pressure [8]. In the 2D case the vessel is assumed to be circular and the change in diameter is uniform.

dp dA

(35)

For small changes in the diameter Dd, in relation to the initial diameter d, the local compliance can be simplified to

p d d p d d d Cl D D × » D -D + = 4 ₂ ) ( 2 2 p p p (1.3)

The distensibility coefficient DC is a mechanical property of the vessel

wall. It is the relative compliance coefficient, which means the relative change in volume due to a change in the local pressure, and is defined as

A C p A A p V V D_C = l D D = D D = (1.4)

The fact that arteries are not perfectly circular and the wall motion is not uniform makes it interesting to study the 3D compliance. We present a novel 3D model of compliance that gives more detailed information since vessels are truly 3D. We define the local 3D compliance at a point on the lumen surface as

p r C_l _D ¶ ¶ = 3 , (1.5)

where ¶r is the translation distance perpendicular to the surface of the vessel. A

3D equivalent to the 2D distensibility does not exist since no assumption about the geometry is made.

(36)

1.4 Objectives and overview

This master’s thesis is a contribution to the field of non-invasive methods for 3D semi-automatic segmentation of vessels, where the Carotid artery with its bifurcation is our main focus. This work was done with the following objectives:

· Develop a 3D semi-automatic segmentation technique for extracting

the lumen surface of the Carotid vessels including the bifurcation from 3D ultrasound examinations.

· Produce a 4D model of the Carotid arteries wall motion.

· Implement an interactive user interface to make it easy for the

physician to use the segmentation technique.

· Produce a 3D model of the compliance of the Carotid vessels to be

used in stroke research and improved assessment of Carotid disease.

· Evaluate the complete method.

Figure 1.10 shows an overview of the segmentation process.

Ultrasound examination

Rawdata Filtered result Initial model

The balloon

model and FEM The resulting 3D model

Planning tool for surgeons Ultrasound

examination

Rawdata Filtered result Initial model

The balloon

model and FEM The resulting 3D model

Planning tool for surgeons

Figure 1.10 The procedure to segment a 3D model from an ultrasound examination. The initial model might also be generated directly from raw data.

(37)

1.4 Objectives and overview 13

The thesis has been divided into the following chapters:

In Chapter 2 – Ultrasound physics and image acquisition, we describe

how images are produced from ultrasound, discuss artefacts that can arise and present techniques that can increase the image quality. Finally, we explain how the ultrasound images were collected and converted into 3D and 4D volumes.

In Chapter 3 - Adaptive Filtering, the filtering of our 3D volume of

ultrasound data is described. The filtering is done since ultrasound images inherently suffer from a high level of noise. A comparison between Nagao-Ye filtering and “Granlund-Knutsson” filtering is performed.

In Chapter 4 - Edge detection, the edge detector used is discussed. Edge

detection is the first step in deriving a 3D potential field, which is used later in the process. We have chosen the Monga-Deriche-Canny detector.

Chapter 5 - Approximation with seeded region growing, describes a

method to obtain an initial model for the final segmentation described below. This reduces the computation time needed for the final segmentation.

Chapter 6 - The balloon model and the Finite Element Method, is the

main chapter of this thesis. Here the final method to produce the segmentation is explained. The theory behind the finite element method is presented here.

Chapter 7 - Implementation, presents an overview of the implementation.

The user interface is described with both screenshots and a short tutorial.

Chapter 8 - Result, contains discussions and reflexions on the outcome of

this thesis. First a short overview in both text and images of the different steps of the segmentation method is presented. Thereafter, images and models from a healthy person and from a patient with stenosis are shown and analysed. Finally, the execution time of the program is discussed.

Chapter 9 - Concluding remarks and future directions, contains

conclusions and some suggestions for future work and extensions of the 3D semi-automatic segmentation technique presented in this thesis.

(38)

(39)

15

2 Ultrasound physics and

image acquisition

2.1 Ultrasound physics

2.1.1 Acoustic impedance

One way to describe acoustic impedance is by analysing the wave equation. This equation tells us that a concentration of pressure (or rather a pressure gradient) will give rise to an acceleration opposite to the pressure gradient, until the difference is equalised [9, 10, 11]. If the partial differential equation is solved, the propagation velocity of the sound wave will be determined by

rk 1

=

c (2.1)

where c is the velocity of the sound wave, ρ is the density and κ is the compressibility of the material. The acoustic impedance Z that is defined as the product between ρ, which is the density of a material, and c, which is the speed of sound in it, can now be expressed as

(40)

16 Chapter 2 Ultrasound physics and image acquisition k r r =× = c Z (2.2)

The reflection of the wave is determined by the difference in acoustic impedance between two materials, which means that what really determines the reflection are the compressibility and the density of the materials. The elasticity of the material makes it want to return to its original shape after being compressed. The acoustic impedance is a measure of resistance to oscillation of the particles in the given material.

2.1.2 Reflection and refraction

When sound waves passes a intersection between two material with different acoustical properties, one part of the sound wave will be reflected while the other part will be transmitted through the interface (Figure 2.1). It is the reflecting portion of the beam that produces ultrasound images. The amount of reflection is determined by the angle of incidence between the sound beam and the reflecting surface. In medical ultrasound, in which the same transducer both transmits and receives ultrasound, almost no reflected sound will be detected if the ultrasound strikes the patient’s surface at an angle of more than 3º from perpendicular (Figure 2.1). Transmitted sound contributes nothing to image formation, but transmission must be strong enough to produce echoes at deeper levels.

(41)

2.1Ultrasound physics 17

Figure 2.1 Reflection and refraction of ultrasound at an interface between muscle and bone.

The reflected portion at tissue interfaces depends on the difference of the two tissues acoustic impedance and the angle of incidence of the beam. If Z1 and Z2

are the acoustic impedances of two different materials and the beam strikes the surface at a nearly right angle, then the reflection fraction R is defined as

2 2 1 2 2 1 ) ( ) ( Z Z Z Z R + -= (2.3)

2.1.3 Attenuation

When a sound wave travels through a material, its energy is attenuated by absorption and scattering exponentially with the depth of travel. Absorption is a result of frictional forces that oppose the motion of the particles in the medium. Energy absorbed from the beam is converted into heat. Scattering is another form of energy loss. When the beam hits a particle smaller than the beam, sound waves bounce in all directions.

Higher frequency results in increased attenuation and therefore smaller depth of penetration. This means that lower frequencies must be used to detect structures deep inside the body. On the other hand, increasing the frequency

(42)

18 Chapter 2 Ultrasound physics and image acquisition

results in increased resolution. Evidently, the proper frequency is a compromise between the resolution and the ability to propagate the energy into the tissue.

2.2 Artefacts in ultrasound imaging

2.2.1 General Comments

Ultrasound images usually suffer from several types of artefacts and low signal-to-noise ratio, both of which limit the image quality. Imaging artefacts appear from inevitable physical processes that weaken the correspondence between the obtained image to be examined and the actual cross-section. In image formation it is assumed that sound travels in straight lines, with a constant velocity and constant attenuation and is reflected only once from each intersection. Neither of these assumptions holds exactly. Some of the most frequent artefacts in ultrasound imaging are summarised below.

2.2.2 Speckle noise

When sound encounters a structure that is much smaller than the wavelength of the beam, the sound is scattered almost equally in all directions. Only a small part will be reflected back to the transducer and detected. The part of the beam reflected in other directions will continue to bounce between interfaces. Interference between the waves scattered from many small structures produces a characteristic pattern called speckle noise. The echo pattern is random and unrelated to the actual pattern of scatters within the structure but it may be sufficiently characteristic to assist in tissue differentiation.

(43)

2.2Artefacts in ultrasound imaging 19

Speckle can be modelled as multiplicative noise. The speckle noise is independent of time and can therefore not be reduced by time averaging. It can be reduced either by compound scanning (see section 2.3.2 Compound scanning) or by increasing the frequency. Unfortunately, by increasing the frequency the penetration ability of the beam is decreased. This can be somewhat compensated for by harmonic imaging (see section 2.3.1 Harmonic imaging).

2.2.3 False echoes

There are two main types of false echoes, namely reverberation and ghost

images from non-perpendicular surfaces. Reverberation echoes occur when the

returning echoes from the second reflection bounces between two interfaces. The returning echo reflects on the back surface of the interfaces where the first echo reflected and initiates a third echo. The transducer interprets this echo as another object that are displayed as a spurious distant structure. Thus three surfaces are displayed, with the third being a reverberation image (Figure 2.2).

1st echo 2nd echo Reverberated echo

1st echo 2nd echo Reverberated

echo

Figure 2.2 False echoes because of reverberations (true object to the left and returned picture to the right).

(44)

A strongly reflecting non-perpendicular surface can reflect an echo from a structure that is not lying in the direction of the ultrasound beam. This echo will be displayed behind the surface as in Figure 2.3.

Figure 2.3 False echo because of strongly reflecting surface not perpendicular to the transducer (true object to the left and returned picture to the right).

2.2.4 Shadowing

The proportions of energy reflected and transmitted depend on the acoustic impedances of the two tissues at the interface. The greater the difference, the greater the percentage reflected. At tissue-air interface, less than 0,1% of the beam is transmitted. This means that almost no energy is left for imaging behind the interface. Gas-filled organs therefore cast a shadow, and structures underneath cannot be displayed. This is the reason why special gel is used between the transmitter and the body so that no air will be present in the very first part of the pathway.

(45)

2.3Improving the image quality 21

2.2.5 Refraction artefacts

In ultrasound imaging the speed of sound in the body is assumed to be the same as for water, which is 1540 m/s. But for some tissues the speed of sound is significantly higher such as for bone (4000 m/s). This will distort the image. Structures behind tissue with high speed of sound will appear to close to the ultrasound source. If only a part of the structure is hidden behind the refraction will cause a distortion in the contour (Figure 2.4).

c3= 2c2

c2

c1

Ultrasound source

Figure 2.4 Spatial distortion due to materials with different speed of sound.

2.3 Improving the image quality

2.3.1 Harmonic imaging

As mentioned above, the resolution improves and the scattering decreases with increased frequency. However, the tissue attenuation increases which makes the

(46)

penetration depth shorter. Harmonic imaging is a method to take advantage of increased frequency and at the same time retain the penetration.

When the sound wave passes through the tissue it will be distorted which causes overtones to appear in the propagating ultrasonic wave. Since the overtones have a higher frequency they will give better resolution and less scattering. Unfortunately there is a trade-off between the scattering of the tone and its energy content and thereby its detectability. Today, only the second overtone is used and therefore the method is often called Second Harmonic Imaging. By filtering the reflected signal with a bandpass-filter so that the second overtone is preserved, an image with higher resolution and lesser scattering can be obtained (Figure 2.5).

Conventional

Harmonic

Conventional

Harmonic

Figure 2.5 Comparison between conventional and harmonic imaging.

2.3.2 Compound scanning

One usual problem with ultrasound images is that only interfaces orthogonal to the beam gives sufficient signal-to noise ratio. The aim of compound scanning is to produce better images by averaging several separate images. If a number of images taken from different angles are averaged the SNR (signal-to-noise ratio)

(47)

2.4 Image acquisition 23

will increase and the speckle is reduced. Another positive effect is that compound scanning is less sensitive to the angle of incidence of the transducer.

Figure 2.6 Comparison between conventional image (to the left) and compound image (to the right).

2.4 Image acquisition

The data were acquired using a System FiVe digital ultrasound scanner (GE VingMed Ultrasound, Horten, Norway) using B-mode configuration. In B-mode (Brightness-mode) the ultrasound system shows slices of the 3D-volume where the intensity corresponds to the strength of the returned echoes. A 10 MHz linear array transducer was translated by hand over the area of the neck where the Carotid artery is located. The translation time was about 15-20 seconds per acquisition and during that time about 400 transversal 2D frames were collected. The size of a complete acquisition is in the order of 50-200 MB. Time-gating is employed and the collected data may be interpreted as a 4D data volume

(48)

f(x,y,z,t) with x, y, z, and time as the four independent variables. More about time-gating follows below.

When the examination is finished, the collection of 2D frames is merged into 3D volumes using a program called EchoPAC 3D (GE VingMed Ultrasound). Once the 3D volume is generated, there are a number of possibilities that enhance the viewing potential of the 2D frames comprising the volume. Since the acquisition is made by hand, the frames will unavoidable become tilted in relation to each other as seen in Figure 2.7. Inaccurate probe movements reduce the quality of the 3D volume since the frames are assumed to be parallel.

By attaching a position detecting devise to the probe, errors due to freehand movement can be corrected for. When using the position-detecting device, a magnetic field is used to record the position and orientation of the frames. The positions and orientations are then used during the merging of the frames. Unfortunately, such a working position detection device was not available.

Figure 2.7 During the merging, the frames are assumed to be parallel as in the left figure. But with freehand acquisition the frames become tilted as in the right figure reducing the quality of the 3D volume. This is corrected for if a position detecting devise is used.

During the examination the pressure of the blood from each heartbeat will deform the vessel wall and therefore the shape of the vessel will be different during the systole compared to the diastole. This produces a problem during the

(49)

2.4 Image acquisition 25

Timegating No timegating

Figure 2.8 The principle behind time-gating. A number of volumes are created, each consisting of frames taken from the same position in the heart cycle. The time-gated volumes are not as dense but do not contain motion artefacts.

Figure 2.9 Comparison between 3D- (left) and 4D-acquisition (right). Since the 4D volume only consists of 15 frames the level of detail in one direction is very low, causing the blurred appearance. In 3D-acqusition mode, all the acquired frames are used to produce the 3D-volum. Therefore the level of detail is much higher at the cost of motion artefacts.

(50)

merging of the 2D frames into 3D volumes. To compensate for this, the ECG signal from the patient was recorded to measure the heart rate and to be able to correspond the pulse to the geometrical changes of the vessel. Each cardiac cycle is divided into a number of time steps. The frames recorded at a given time step in the cardiac cycle are put in the same volume according to the patients ECG signal, so that each time step results in a volume (time-gating). The result will thus be a number of 3D volumes of different time steps from the ultrasound acquisition as in Figure 2.8. This eliminates the movement artefacts. In 4D acquisition mode the EchoPAC program divides the cardiac cycle into 26 time steps, each consisting of frames taken from the same position in the heart cycle. Therefore each 3D volume will consist of only 15 frames. In total this typically amounts to 26x15x512x512 = 100 Mb. A comparison between 3D- and 4D-acquisition mode is given in Figure 2.9. To get good 4D results more frames would be needed. Unfortunately, the memory size in the ultrasound equipment prevented collection of additional frames.

3D-acquisition has been used for all the results in this work except for the animations of wall-motion and compliance.

(51)

27

3 Adaptive Filtering

Because of the high level of noise in ultrasound images, it is necessary to filter the images before further processing. Filtering is always a trade-off between noise suppression and loss of information, something that physicians are often very concerned about. It is therefore attractive to keep as much of the important information as possible.

The filtering should be performed in 3D to achieve good noise suppression while preserving the sharpness. If a 2D kernel would be used, the kernel would have to be wider to give the same noise suppression and would therefore give a blurrier image. Since the original data is actually 4D (3D + time) it would be preferable to perform 4D filtering. However, this would take considerable time on present workstations and require a large amount of memory and we have therefore not investigated this further.

We have tested two adaptive filtering methods for edge preserving noise suppression: The Nagao-Ye method [12] and the Granlund-Knutsson method [13, 14, 15].

(52)

28 Chapter 3 Adaptive Filtering

3.1 Nagao-Ye filtering

A simple method for edge preserving and noise suppression is the Nagao-Ye method. It works by studying a number of subregions around each pixel and setting the value of that pixel to the mean of the subregion with the lowest variance. Nagao first proposed the method and Ye later changed the subregions to reduce the computations needed.

In the 3D case 6, 18 or 26 connectivity can be used, resulting in 7, 19 or 27 subregions. It was found that 6-connectivity worked best for our applications and the subregions for 6-connectivity is shown in Figure 3.1.

Figure 3.1 The seven subregions used in Nagao-Ye filtering with 6-connectivity.

The method can be described by

å

W = i f n i 1 (a,b,c) m (3.1)

(

)

[

å

W -= i i i2 _n1 f a,b,c m 2 s

]

(3.2) (3.3)

[

m k s s s s =min ₀, ₁,K, (3.4) ) , , ( ) , , (x y z x y z fnagao-ye =mk

(53)

3.2 Granlund-Knutsson filtering 29

where:

n is the number of voxels in each subregion (9 in the 3D-case).

m is the number of subregions (7 for 6-connectivity).

m 2 x,

It can be seen that seven means and variations are used for calculating each voxel. However, each mean- and variation-value are used seven times (for

6-connectivity) which means that and s can first be computed for each

3x3x3 subregion and stored in the result volumes and s .

Then the result is obtained by (3.3) and (3.4). The process is

illustrated by Figure 3.2. i m ) z 2 i

(

x,y,z

)

(

y,z

)

, , (x y fnagao-ye 3x3x3 f(x,y,z) m(x,y,z) s2_(x,y,z) Equation (1.1) (1.2) Equation (1.3) (1.4) fnagao-ye(x,y,z) 3x3x3 f(x,y,z) m(x,y,z) s2_(x,y,z) Equation (1.1) (1.2) Equation (1.3) (1.4) fnagao-ye(x,y,z)

Figure 3.2 The two steps of Nagao-Ye filtering. First the mean and variance is calculated for each 3x3x3 subregion. Then, the mean of the subregion with the lowest variance around each voxel is taken as the output value.

The process can be iterated to obtain more noise reduction. Regular or gaussian

mean can be used for computing . Tests were made that showed that

the best results were obtained with 6-connectivity and one iteration. Gaussian or regular mean gave no apparent difference. A result is shown in Figure 3.4.

(

x,y,z

(54)

3.2 Granlund-Knutsson filtering

Granlund and Knutsson have proposed an adaptive filtering technique based on orientation estimation using tensors. The basic principle behind this technique is to adapt the shape and orientation of the filter so that smoothing is only performed parallel to edges, not across them, thus preserving the sharpness of the edges. The technique can be applied to signals of any dimension. In this section,

only a brief description is given. A longer description is given in Appendix A -

Estimation of orientation, Appendix B - Control of filter shape and Appendix C - Synthesis of the adaptive filter.

In the first step, local structures and orientations in the image are estimated using a set of quadrature filters with different directions. The estimated orientations are used in the second step to compute the desired shape and orientation of the filter. In the third step, the output is generated as the weighted output from a set of shift invariant filters. One ordinary lowpass filter is used together with a number of highpass filters with different orientations. By weighting the output from each filter it is possible to adapt the filter to the input. The process is outlined in Figure 3.3.

The adaptive filter is constructed from shift-invariant fixed filters Flp and Fk according to (3.5)

( )

+

( )

=

( )

+

_å

( )

= k k k lp hp lp F F c F F F(u,C) r u,C r u

where u is the N-dimensional frequency, r= u is the radial frequency and C

is the set of control tensors. lp F Fhp hp F k F ck

is a fixed non-adaptive spherical lowpass filter while is a

high-pass filter that adapts to the signal, as mentioned earlier. consists of six

(55)

Estimation of orientation Control of filter shape

Set of quadrature filters LP Tensor mapping ( j) mTlp,s,a,b, (r,a,b,j) m å = k k k q T M Construct orientation tensor T Tlp C 1 M 2 M 3 M 1 q 2 q 3 q Input signal

Synthesis of the adaptive filter

Calculate weighting coefficients k F k c=C·M, 1 , F M 2 , F M 3 , F M F4 4 , F M · · · · S 1 c c2 c3 c4 Set of shift invariant filters Flpand Fk Output signal FLP F3 F2 F1 C Input signal

Estimation of orientation Control of filter shape

Set of quadrature filters LP Tensor mapping ( j) mTlp,s,a,b, (r,a,b,j) m å = k k k q T M Construct orientation tensor å = k k k q T M Construct orientation tensor T Tlp C 1 M 2 M 3 M 1 q 2 q 3 q Input signal

Synthesis of the adaptive filter

Calculate weighting coefficients k F k c=C·M, 1 , F M 2 , F M 3 , F M F4 4 , F M · · · · S 1 c c2 c3 c4 Set of shift invariant filters Flpand Fk Output signal FLP F3 F2 F1 C Input signal

Figure 3.3 Overview of the Granlund-Knutsson filtering technique. The input signal is filtered with a set of lognormal quadrature filters with different directions to produce orientation tensors T. These tensors are then lowpass filtered into TLP and mapped into control tensors C. The control tensors are used

to do a weighted summation of the output from a set of shift invariant filters to produce the output signal. (The number of quadrature filters in the first step and the number of fixed filters in the third steps, depends on the dimension of the signal. The number of filters used here might be used in the 2D case).

(56)

k

C

for each corresponding filter. depends on the orientation of each filter and on

the control tensor which depends on the estimated orientation. c

To construct the control tensors C, an estimation of local orientation is first done by six lognormal quadrature filters with the same directions as the

highpass filters Fk above. A representation using tensors are used to describe

orientation according to (3.6) T T T 3 3 3 2 2 2 1 1 1eˆeˆ eˆ eˆ eˆ eˆ T=l +l +l

where êi are normalized eigenvectors with corresponding eigenvalues li. Each êi

represent an orientation and li can be thought of as the strength of that

orientation. The orientation tensor T is constructed by

=

å

=

å

(

-

)

k T k k k k k k q q M n_ˆ n_ˆ ₅1I T (3.7)

where qk is the output from quadrature filter k, is the dual orientation tensor

for filter k and n is the direction of quadrature filter k.I is the unity tensor.

k M k ˆ 1 ˆ ˆe₂ 3 ˆ

Before the orientation tensors are remapped into control tensors, a lowpass filtering of the orientation tensor field is performed to get good estimates of orientation even at low SNR. It should be noted that this lowpass filtering does not reduce the sensitivity to high frequencies in the original volume. It only suppresses rapid changes in orientation.

Each orientation tensors can be said to belong to one of three cases:

· The plane case (l1>>l2»l3, T»l1eˆ e1ˆT1 ): This case corresponds to a neighbourhood that is approximately planar, i.e. the signal is constant on parallel planes. e is normal to the plan while the eigenvectors and e span the plane.

(57)

·

The line case (l1 » l2 >> l3, T» l1(ˆe1eˆ1T +eˆ2eˆT2) 3

ˆe ˆ1 ˆ2

): This case

corresponds to a neighbourhood that is approximately constant on

parallel lines. The line is parallel to ± and e and e are

perpendicular to the line.

· The isotropic case (l1»l2»l3, T»l1I): This case corresponds to an approximately isotropic neighbourhood, meaning that there exist energy in the neighbourhood but in no typical orientation, e.g. in the case of noise.

In reality few orientation tensors can be said to belong to only one case and therefore the following three measures are used to measure how close a tensor is to each case (a value of one represents a perfect match).

case isotropic The P case line The P case plane The P 1 3 3 1 3 2 2 1 2 1 1 l l l l ll l l = -= -= (3.8) where 0£Pk£1 and P1+P2+P3=1.

The control tensor C is constructed by remapping the lowpass filtered orientation tensor Tlp. Often C is simply an adaptively scaled version of Tlp.

The weighting coefficients ck for each highpass filter Fk are computed as

(3.9)

k F k

(58)

where MF,k is the dual orientation tensor for filter Fk. The filter output is

produced as a summation of the output from each shift invariant highpass filter Fk weighted by ck.

3.3 Results

A 100x100x100 volume took approximately two minutes to filter using the Granlund-Knutsson technique on an 800 MHz Pentium III system. The Nagao-Ye method was approximately 30 times faster than the Granlund-Knutsson technique.

A comparison between the original image and images filtered with Gaussian, Nagao-Ye and Granlund-Knutsson filtering is shown in Figure 3.4.

The noise suppression in the gauss filtered image in Figure 3.4 (b) and the image

obtained by Granlund-Knutsson filtering (d) are approximately equal in areas

without any dominant structures (e.g. the dark blue areas). However, the Granlund-Knutsson filtering preserves more details of the edges and does not have the lowpass filtered appearance that (b) has.

The results from the Nagao-Ye filtering have a “blocky” appearance and the edges are more jagged. Physicians that we talked to preferred the images produced by the Granlund-Knutsson technique since the Nagao-Ye method altered the appearance of different types of tissue. This makes it harder for physicians to distinguish between plaque and other kinds of tissue when using the program. Therefore, the Granlund-Knutsson technique has been used in the remaining part of this work even though the Nagao-Ye method is considerable faster.