Enhancement and Visualization of VascularStructures in MRA Images Using Local Structure

(1)

Institutionen för medicinsk teknik

Department of Biomedical Engineering

Examensarbete

Enhancement and Visualization of Vascular

Structures in MRA Images Using Local Structure

Examensarbete utfört i Reglerteknik vid Tekniska högskolan i Linköping

av

Morteza Esmaeili

LiTH-IMT/MASTER-EX--10/003--SE

Linköping 2010

Department of Biomedical Engineering Linköpings tekniska högskola

Linköpings universitet Linköpings universitet

(2)

(3)

Enhancement and Visualization of Vascular

Structures in MRA Images Using Local Structure

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Morteza Esmaeili

LiTH-IMT/MASTER-EX--10/003--SE

Handledare: Mats Andersson

imt, Linköpings universitet

Examinator: Hans Knutsson

imt, Linköpings universitet

(4)

(5)

Avdelning, Institution Division, Department

Division of Medical Informatic Department of Biomedical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2010-006-16 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport

URL för elektronisk version

http://www.imt.liu.se/mi/ http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-ZZZZ ISBN — ISRN LiTH-IMT/MASTER-EX--10/003--SE Serietitel och serienummer

Title of series, numbering

ISSN —

Titel Title

Svensk titel

Enhancement and Visualization of Vascular Structures in MRA Images Using Lo-cal Structure Författare Author Morteza Esmaeili Sammanfattning Abstract

The novel method of this thesis work is based on using quadrature filters to es-timate an orientation tensor and to use the advantage of tensor information to control 3D adaptive filters. The adaptive filters are applied to enhance the Mag-netic Resonance Angiography (MRA) images. The tubular structures are extracted from the volume dataset by using the quadrature filters. The idea of developing adaptive filtering in this thesis work is to enhance the volume dataset and suppress the image noise. Then the output of the adaptive filtering can be a clean dataset for segmentation of blood vessel structures to get appropriate volume visualization. The local tensors are used to create the control tensor which is used to control adaptive filters. By evaluation of the tensor eigenvalues combination, the local structures like tubular structures and stenosis structures are extracted from the dataset. The method has been evaluated with synthetic objects, which are vessel models (for segmentation), and onion like synthetic object (for enhancement). The experimental results are shown on clinical images to validate the proposed method as well.

Nyckelord

Keywords Adaptive Filtering, Tubular Structure, 3D Filter, Enhancement, Local Structure Tensor, Orientation Estimation, Vessel Segmentation

(6)

(7)

Abstract

The novel method of this thesis work is based on using quadrature filters to esti-mate an orientation tensor and to use the advantage of tensor information to con-trol 3D adaptive filters. The adaptive filters are applied to enhance the Magnetic Resonance Angiography (MRA) images. The tubular structures are extracted from the volume dataset by using the quadrature filters. The idea of developing adap-tive filtering in this thesis work is to enhance the volume dataset and suppress the image noise. Then the output of the adaptive filtering can be a clean dataset for segmentation of blood vessel structures to get appropriate volume visualization. The local tensors are used to create the control tensor which is used to control adaptive filters. By evaluation of the tensor eigenvalues combination, the local structures like tubular structures and stenosis structures are extracted from the dataset. The method has been evaluated with synthetic objects, which are vessel models (for segmentation), and onion like synthetic object (for enhancement). The experimental results are shown on clinical images to validate the proposed method as well.

(8)

(9)

Acknowledgments

I would like to thank Prof. Hans Knutsson for giving me the opportunity of work-ing with him on this thesis work. He always has new and valuable ideas and by his professional guidance, he helped me find the correct way.

I would also like to thank my supervisor Dr. Mats Andersson for guiding and supporting me in every step of this thesis work and encouraging me to keep going forward.

Finally, I would like to thank my opponent Yan Xie for reading my report and giving valuable suggestions.

Morteza Esmaeili Linköping, 16 June 2010

(10)

(11)

an almost isotropic tensor, λ1 = 0.67, λ2 = 0.63, λ3 = 0.59. The isotropic case is represented in green sphere (its radius corresponds to the magnitude of the isotropic part), the planar case is repre-sented in red spear through the poles (the length of the spear is proportional to the magnitude of T1), and the plane case is repre-sented in yellow expanded sphere orthogonal to e3direction. see [11] for details. Illustration from J. Wiklund et al. (2006), T-FLASH: Tensor Visualization in Medical Studio. . . 23 3.6 volume visualization of the spherically symmetric synthetic test

data which is used for evaluation. . . 26 3.7 m − f unction with β = 4, σ = 0.2, j = 1, α = 0, 0.2, 0.4, 0.6, 1.

(bottom to top respectively), x ∈ [0 1]. . . 27 3.8 Guassian low-pass filter with standard deviation 1. grid size is 15×15. 28 3.9 µ − f unction with α = 0.4, j = 1, β = 1, 2, 4, 6, 12. (left to right

respectively), x ∈ [0 1], and in bottom µ−f unction with β = 1, j = 1, α = 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 0.9, 0.95.(from top to bottom),

x ∈ [0 1]. . . . 30 3.10 The γ1output (left), which shows the smooth changes on the edges

from bright area to dark (σ = 0.2, α = 0.1, β = 4, j = 1). The

µ2 output (right), which has smooth changes on the edges (α = 0.45, β = 3, j = 1). . . . 31

(14)

2 Contents

3.11 The volume test data used for evaluation, 8 slices (top), The cor-rupted data by white Gaussian noise, with SNR = 0 (bottom). . . 33 3.12 Result from 3D adaptive filtering on the volume test data, 8 slices. 34 3.13 3D visualization of synthetic data with additive Gaussian noise,

SNR = 0,(left) and after adaptive filtering (right). . . 34 3.14 Axial view of the T1-weighted MRA image, slice 90, from top to

bottom original image, corrupted image with SNR = 0, 5, and 10. From left to right noisy image and the result from the 3D adaptive filtering with size 9 × 9 × 9 in spatial domain. . . 36 3.15 Top from left to rigth: The weight function which is constructed

by the control tensor and the result of the enhancement filters to the input which is the white Gaussian noise after 20 iterations. Bottom: the sum of enhancement band-pass filters after 5 iterations in frequency domain, with filter size 7 × 7 × 7 in spatail domain. . 37

4.1 m − f unction with β = 3, α = 0.1, j = 1, δ = 0.1, 0.2, 0.4, 0.6 (top

to bottom respectively), λ2−λ3

λ1 ∈ [0 1]. . . 40

4.2 Second-order 3D variation ellipsoid shape space. The classification is based on eigenvalues of the Hessian matrix. The bright string prototype indicates vessel variations [12, 18]. Illustration from Lin Qingfen,(2003), Enhancement, Detection, and Visualization of 3D Volume Data [18]. . . 42 4.3 Top from left to right: The cross section (Axial cut) of the 3D

tubu-lar structure, The intermediate processing result by low-scale Gaus-sian filter, large-scale filter, and the GVF. Middle: Corresponding gradient vector field representation. Bottom row: Corresponding normalized gradient vector field representation. Illustration from Bauer et al.(2008), A Novel Approach for Detection of Tubular Ob-jects and Its Application to Medical Image Analysis [9]. . . 44 4.4 m − f unction with δ = 0.2, β = 5, α = 0.1, j = 1, λ2−λ3

λ1 ∈ [0 1]. . 45

4.5 The synthetic objects with varying configuration. From top to bot-tom: original dataset, Frangi’s method, Krissian’s method, Bauer’s method [9], new method with ratio PT S1and by using local tensor

eigenvalues, in the last row the same mapping parameters and eigen-values of control tensor are used. For both δ = 0.2, α = 0.1, β = 5, j = 1. From left to right: junction with constant diameter, T-junction with varying diameter, Y-T-junction with constant diameter, Y-junction with varying diameter. . . 46 4.6 The synthetic objects with varying configuration and from different

method response. From top to bottom: original dataset, Frangi’s method, Krissian’s method, Bauer’s method, the proposed method. From left to right: tube with varying diameter, tangenting tubes, and helix. Illustration from C. Bauer et al. (2008), A Novel Ap-proach for Detection of Tubular Objects and Its Application to Med-ical Image Analysis [9].(for the helix object gamma correction = 0.4 is used to be printable) . . . 47

(15)

Contents 3

4.7 (a) The synthetic tubes with varying diameter and contrast with size 60 × 40 × 70 (b) and the result from tubular extraction, mapping function arguments are: Ratio PT S2, δ = 0.2, α = 0.1, β = 3, j =

1. (gamma correction = 0.4 is used to be printable), Right column: cross-section view, slice z = 45. . . 49 4.8 Middle slice of the synthetic T-junc. (a) Tube with additive

Gaus-sian noise, SNR = 0 ,(b) discontinuity appeared by using ratio PT S1.

(c) The result from ratio PT S2. For both situations mapping

func-tion arguments are: δ = 0.2, α = 0.1 , β = 4, j = 1.(gamma correction = 0.4 is used to be printable.) . . . 50

4.9 m − f unction with δ = 0.2, β = 4, α = 0.1, j = 1. √ λ2 2−λ 2 3 λ1 ∈ [0 1]. 51

4.10 Implementation of tubular extraction using PT S2 (left) and PT S3

(right) on synthetic object T-junction (60 × 60 × 60) by the map-ping function arguments: δ = 0.9, α = 0.1, β = 5, j = 1.(gamma correction = 0.5 was used to become printable) . . . 52 4.11 The probability of tubular structure using ratio PT S3. The plot

demonstartes the values in the middle line (yellow line) of the figure 4.10 (b). . . 52 4.12 Visualization of the synthetic objects Y-junction (64 × 64 × 64).

From top to bottom 3D visualization of the corrupted Y-junction object with SNR = 0, slice z = 30 with SNR = 0, slice z = 30 with SNR = 5, and slice z = 30 with SNR = 10. From left to right columns; original image, the tubular extraction result using local structure tensor, and the result using control tensor. (gamma correction = 0.4 was used to become printable) . . . 54 4.13 The MRA of carotid, with volume size of 256 × 256 × 91 (top) the

µ − f unction with parameters: α = 0.15, β = 5, j = 1. (down) . . 55 4.14 The MRA of carotid with volume size of 256 × 256 × 91, the original

image corrupted by Gaussian noise with SNR = 0 (top). The result of tubular extraction without adaptive filtering implementation, the filter size in spatial domain is 13 × 13 × 13 (middle). The result of tubular extraction after adaptive filtering implementation and the filter size in spatial domain is 13 × 13 × 13 (bottom). At the right the corresponding µ − f unctions are shown. . . . 56

4.15 m − f unction with δ = 0.2, α = 0.1, β = 2, j = 1, √

λ2 2−λ23

λ1 ∈ [0 1]. 57

4.16 The MRA image of carotid vascular, (a) The tensor magnitude is averaged by glpwith the spatial size 5×5×5 and standard deviation

σ = 1.5. (b) The tensor magnitude is averaged by glp with spatial

size 15 × 15 × 15 and standard deviation σ = 1.5. The mapping function argumnets are: PT S2, δ = 0.2, α = 0.1, β = 3, j = 1. The

green arrows point to the large vessels and the red arrows point to the small vessels. . . 58

(16)

4 Contents

4.17 The top row shows the 3D visualization of the synthetic data with narrowing in its bifurcation and the possible positions for steno-sis( from left to right), bottom row shows the middle slice of data, and the transfer function. The mapping function argumnets are:

Pstenosis, δ = 0.5, α = 0.1, β = 5, j = 1. . . . 60

List of Tables

3.1 Results from RMS angular measurement . . . 26

4.1 Various structures obtained by Fragie’s method, λk eginevalues, H

= high, L = low, N = noisy, +/- indicates the sign of the eigenvalues. Table from Frangi et al., multiscale vessel enhancement filtering [12]. 43

(17)

Chapter 1

Introduction

1.1 Motivation

Many researches have focused on vessel segmentation of volumetric datasets, which are obtained by Magnetic Resonance Angiography (MRA) or Computed Tomog-raphy AngiogTomog-raphy (CTA). The result of vessel segmentation can be useful for visualization, modeling of the vessel structure, studying its pathologies for diag-nosis, and surgical planning. Several methods have been introduced for vessel segmentation like: Intensity based methods, region growing [16], level-set [17], wa-tershed [13], geometrical model approach [12], and hybrid approach [9] which is the combination of two methods, etc.

Most imaging techniques produce a 3D dataset, so it is desired to process the whole 3D data instead of 2D processing. In order to process 3D image dataset, first enhancement procedure should be developed to suppress the image noise from the volume dataset. The method used in this thesis work is based on using multi-dimensional adaptive filtering. The local structure orientation tensor is estimated by using quadrature filters introduced by Knutsson et al. (1989). By using the local orientation tensor eigen-system, the control tensor will be estimated. The control tensor will be used to control adaptive filtering which provides enhancement processing. So after enhancement, the other processing step can be implemented efficiently.

Considering the eigenvalues of the tensor, the tubular structures can be ex-tracted from an MRA volume. Both enhancemet and segmentation methods which are proposed in this thesis work are based on orientation tensor estimation.

1.2 Objective

The aim of this thesis work is to develop a set of post-processing approaches, which consists of adaptive filtering and segmentation of vascular structures, to get better vessel visualization. So the 3D adaptive filter will be implemnted on a 3D MRA dataset to provide enhanced data volume. Then the tubular structures which

(18)

6 Introduction

introduce vessels in the volume will be extracted. The 3D phase invariant adaptive filters are designed by the combination of a low-pass filter and controlled band-pass filters. The band-band-pass filters control the amount of high-frequency contents in the local structures by using the control tensor. Considering the eigenvalues of the local tensor, the tubular structure in the dataset will be extracted and finally the modified dataset will be visualized.

1.3 Background

This thesis work is presented in three main chapters which are based on filter optimization, estimation of local orientation tensor and adaptive filtering, and segmentation, in 3D space. The filter optimization is based on the mean squared distance estimation, which is mainly combined by previous works from Knutsson et al. (1995, 1999) [7, 10]. The orientation tensor estimation and the synthesis of adaptive filtering are mainly based on the previous works and publications by Knutsson et al. (1989, 1995) [15, 7] , Mats Andersson (1992) [4], and Farnebäck et al. (2002) [1]. In order to investigate a 3D segmentation method, various strategies have been studied but finally the tubular structure extraction is used. The idea was to apply the synthesized adaptive filtering output to extract vascular structures based on using the defined ratios by Knutsson et al. (1989, 1995) [15, 7] and it is compared with the previous works related to tubular extractions by Frangi et al. (1998) [12], Bauer et al. (2008) [9], Lin Q. (2003) [18], and Krissian et al. (1999) [14]. Here the main previous publications and documentations which are used to write this thesis work are listed.

Knutsson H., Granlund G.H. (1995). Signal Processing for Computer Vision. Kluwer Academic, Dordrecht, the Netherlands.

Knutsson H. (1989), Representing Local Structure Using Tensors, The 6th Scan-dinavian Conference on Image Analysis, Oulu, Finland, pages 244-251.

Knutsson H., Andersson M., Wiklund J., Advanced filter design, in Proceedings of the Scandinavian Conference on Image analysis, Greenland, June 1999, SCIA.

Andersson M. (1992), Controllable Multi-dimensional Filters and Models in

Low-level Computer Vision. Ph.D. thesis, Linköping University, Linköping,

Swe-den.

Farnebäck G. (2002), Polynomial Expansion for Orientation and Motion

Esti-mation. Ph.D. thesis, Linköping University, Linköping, Sweden.

Bauer C., Bischof H., A Novel Approach for Detection of Tubular Objects and

Its Application to Medical Image Analysis. In: Delp, Rigoll (eds.) DAGM 2008.

LNCS, vol. 5096, pp. 163-172. Springer, Heidelberg (2008)

Svensson B. (2008), A Multidimensional Filtering Framework with Applications

to Local Structure Analysis and Image Enhancement. Ph.D. thesis, Linköping University, Linköping, Sweden.

1.4 Outline

(19)

1.4 Outline 7

Chapter 2: This chapter consists of the overview of the basic theory behind the various parameters engaged in filter optimization process. It is mentioned how these parameters effect the optimal solution to get the desired frequency response. The method and algorithm which are used to optimize the filter basis are discussed in this chapter.

Chapter 3: The local orientation tensor and the control tensor estimation is discussed in this chapter. It also contains the basics of the tensor properties. The control tensor, the related mapping functions, and generalization of enhancement filters are reviewed. The result of 3D adaptive filtering is evaluated with synthetic and clinical data.

Chapter 4: The characteristics of tubular and stenosis structures and their relationship with the orientation tensor eigenvalues are reviewed in this chapter. 3D segmentation method is based on the geometric properties of the neighborhoods in the dataset. The method is implemented on synthetic and clinical dataset and evaluated in this chapter.

Chapter 5: Conclusion and suggestion for future related works is discussed in this chapter. Also the limitations and possible optimization issues are reviwed.

(20)

(21)

Chapter 2

Filter Optimization

2.1 Introduction

This chapter provides a background about how the filter is optimized to maintain the desired frequency response and spatial locality. It will be discussed how the mean squared distance criterion is used to get the desired frequency response and optimization. Also the input arguments of the filter optimizer function and the optimal characteristic of each argument will be discussed. The implementation of the filters on the input signal in spatial domain can be performed by convolution which is given by:

h[n] = x[n] ∗ y[n] =

+∞ X

m=−∞

x[m]y[n − m] (2.1)

where x is the input signal and y is the filter function, and by considering the Fourier properties it can be defined as:

H(u) = X(u)Y (u) (2.2)

So the other approach to implement the filter on the input signal is the multi-plication of their Fourier transforms. Supposing that h[n] is an impulse response and H(u) is its Fourier transform (equation 2.3, 2.4), it has a continuous form in the Fourier domain and discrete form in the spatial domain [7]. Considering the equations 2.1 and 2.3, the output of the function h[n] is dependent on the interval [−π, π] in the Fourier domain [19].

h(x) = 1/2π

π

Z

−π

H(u) exp (iuTx)du (2.3)

H(u) =X

x

h(x) exp(−iuTx) (2.4)

(22)

10 Filter Optimization

2.2 Representation Spaces

In the filter design process, two representation spaces (domain) are considered; the spatial domain and the Fourier domain [7]. In most cases the spatial domain is given by the sampled input signals [10]. The desired frequency response is considered in the Fourier domain and will be optimized by kernel coefficients. The error, which is based on the mean squared problem, will be minimized respect to ideal filter [7].

2.3 Filter Optimization Components

The filters are optimized using kernel optimization (written by Wiklund J. et al.) procedure which includes the following parameters:

• Frequency ideal function • Frequency weight function • Spatial mask

• Spatial weight function • Spatial ideal function

The output of kernel optimizer is the filter both in the Fourier domain (FD) and the spatial domain (SD). To design the filter, there are several objectives that must be met: to maintain the ideal frequency function, to have small spatial support, to adapt it to the image spectrum and to adapt it to the spectrum of the noise [7]. These cannot all be met at the same time, so we have to compromise to as small error as possible respect to the ideal function.

The different features of filter optimization perform in different domains. To maintain the locality and limit the effective spatial size a spatial ideal function, which is the most compact function (impulse), and a spatial weight function, which is usually a quadratic function (equation 2.5), are used (figure 2.1) [10].

w(x) = |x|n, n > 0 (2.5)

There are adavantages in preserving the spatial locality. The computational load will be decreased by using a small filter [3]. The other important effect of using smaller filters is increasing the resolution of the signal. So by choosing small effective size for filters, they can be more fit to the signal locally [3].

2.4 Distance Measurement

The error is defined as weighted mean squared distance (equation 2.6) which must be minimized to get filter function as close to the ideal filter function as possible.

(23)

2.4 Distance Measurement 11

Figure 2.1. From left to right: the spatial weight function and spatial ideal function

respectively, for filter with spatial size 7 × 7.

min 2=

N

X

n=0

kWn(fn− Bf0)k2 (2.6)

where Wn are the diagonal weighting matrices, B is the Fourier basis matrix, and

fnis the ideal kernel function [7, 10]. In order to minimize the error function,

par-tial derivative respect to the kernel coefficients (f0) can be calculated as following:

min 2= N X n=0 (fn− Bnf0)TWn2(fn− Bnf0) (2.7) ∂2 ∂f0 ≡ ∂ 2 ∂Re(f0) + i ∂ 2 ∂Im(f0) = 0 (2.8) ∂2 ∂f0 ≡ 2Re( N X n=0 BT_nW_n2Bnf0− BnTW 2 nfn) + 2Im( N X n=0 B_nTW_n2Bnf0− BnTW 2 nfn) = 0

The solution [10] of the above equation is given:

N X n=0 BT_nW_n2Bnf0= N X n=0 B_nTW_n2fn) If we define: A = N X n=0 B_nTW_n2Bnf0 and h = N X n=0 B_nTW_n2fn)

(24)

then the optimized kernel coefficient [10] can be estimated as:

f0= A−1h (2.9)

The expected spectrum related to the radius in the Fourier domain, ρ = kuk (nor-malization of spectral energy in radial direction), can be described [7] as following:

S(ρ) = ρα , α < 0 (2.10)

The frequency weight should reflect the image spectrum which contains more important frequencies in natural images. It should be considered that the signif-icant spectral energy is concentrated in the center of the Fourier domain and it will decrease by increasing the ρ [7]. Also if we consider the spectrum energy, the imaging in various scales does not lead to a big difference in energy distribu-tion shape, it just changes the magnitude [7, 10]. The equadistribu-tion 2.10 is a class of functions which fulfill these two properties.

Knutsson et al. suggested the equation 2.11 for weight function to minimize the mean squared error of the kernel function. In order to maintain frequency weight function for our 3D filters the following equation is used and α = −2 is chosen to make sure the expected spectrum will be weighted.

W (u) = ραY

k

cosβ(uk/2) + n (2.11)

where n is the expected signal to noise ratio (SNR) which is chosen 10 here, the

cosine term (β) controls the degree of band limiting of the signal [7, 10]. Also the filters should be insensitive to the local DC-level (average intensity of the neighborhoods of the signal) of the signal which can be maintained by setting the weight W(0) to a very large value (figure 2.2) [3, 10].

(25)

2.4 Distance Measurement 13

Figure 2.2. Log scale of the frequency weight function with a very large value (middle

surface of 3D weight function).

The real part of the filter detects lines and the imaginary part detects edges (figure 2.3). In order to detect lines and edges in all directions, for a 2D image we need at least three different filters and for 3D it should be at least 6 filters. The number of filters should be more than 2N −1 [7] where N is the filter dimention. For example in 2D with two orthogonal filters, we cannot detect structures that are oblique and pass through the origin of these two filters. But with at least three filters this problem will be resolved. For 3D case more than 4 filters are needed but to get symmetric fashion distribution [7], it should be at least 6 in icosahedron volume (it will be shown in later chapter). With more filters, it is possible to increase the quality of the result even more at the cost of computing time.

Figure 2.3. Resulting of even (left) and odd (right) parts of kernels. the grid size is

(26)

2.5 Computational Complexity

The computational time for filtering is one of the critical challenges in multi-dimensional signal processing and causes restriction in choosing filter size in spa-tial and Fourier domains [7]. Since implementation of filters to multi-dimensional signals can be done by convolution in spatial domain or multiplication in Fourier domain, the number of multiplications should be considered. This becomes more critical when we want to design a three dimensional filter in which three dimen-sional convolution is needed and it would face with a large dataset volume. One approach is to design a good convolution function which reduces the computa-tional time (as in this project the conv3 function from spatio-temporal toolbox written by Farnebäck G. is used). Also as the real part of the enhancement filter is considered when applying on the volume dataset, the imaginary part can be ignored and the number of multiplications will be reduced.

Figure 2.4. Frequency response of quadrature filter, the center frequency is 2π/5 and

bandwidth is 2 [octaves].

2.6 Discussion

The suggested mean squared distance (Knutsson et al. 1995) is used to opti-mize the filters respect to ideal frequncy responce. In the optimization process the frequency weighting and localization of the filters in the spatial domain are considered. The spatial size is kept small to preserve locality. The optimization coefficients are chosen manually to get optimal response. This can be performed automatically by developing the coefficient measurement algorithm; like the pre-vious work developed by M. Langer et al., who used Genetic algorithm to design multi-dimensional filters [19].

(27)

Chapter 3

Local Orientation

Estimation in 3D Space

3.1 Introduction

In this chapter the theory and generalization of the local orientation tensor will be reviewed. The algorithm used to represent the local orientation is based on quadrature filters and tensor representation which provides the study of multi-dimensional data. The local orientation tensor can be estimated by using the polar separable quadrature filters [4]. The phase and magnitude of the quadrature filters, which are applied to multi-dimensional data, are used to study the symmetry properties of the signal and estimation of orientation tensors respectively [7]. Also the various structures included in dataset can be studied by tensors information which will be discussed in the next chapter as well.

3.2 Tensor Basics

Some physical quantities can be expressed as scalars like mass, area, etc. and quantities in which direstion is needed are expressed as vector like velocity, flow, etc. But many physical quantities and properties cannot be expressed as scalar or vector like inertia deformation, anisotropic conduction, etc. Such quantities can be expressed by tensor (T) in the form of array of the numerical value [7]. It is described by various orders depending on the complexity and dimensionality of quantities representation. Scalar and vector themselves are zero and first order tensor respectively. Higher order tensor can express the complex properties (ex. 2nd order can express inertia deformation). Here some properties of tensor are given: Scalar Product T • U =X ij tijuij (3.1) 15

(28)

16 Local Orientation Estimation in 3D Space

where, T and U are two tensors, and tij, uij are elements of tensors T and U

re-spectively [7]. The scalar product of T with itself can be defined as norm (λn are

eigenvalues of T). Frobenius Norm kT k2= T • T =X ij t2_ij=X n λ2_n (3.2) Basis Tensors T =X k ckBk, (3.3) Orthogonalbasis ( Bi• Bj = 1 if i = j, Bi• Bj = 0 if i 6= j.

Tensor can also be written as weighted basis where ck are the coefficient of the

basis tensors Bk [7].

Outer product tensors [7] from two vectors a and b.

T = ¯a¯bT (3.4)

The norm of the outer product [7] can be represented as:

kT k2= T • T = (a • a)(b • b) = kak2kbk2 (3.5)

3.3 Orientation Estimation

For the non-constant and real valued simple signal, orientation can be expressed as:

s(x) = g(xTn)ˆ (3.6)

where g is a real valued non-constant function and ˆn is the normalized vector

oriented along maximum variation of the signal [7]. But the classes of simple sig-nals are limited, and it is needed to define more general declaration for orientation of signals [1]. Two various neighborhoods with their energy contribution are shown in figure 3.1. In 3D volume, there are various structures such as line-like, plane-like, isotropic regions (blob-like) [7]. The quadrature filters with various directions is needed to study the local structures. The magnetiude of these filters will be used to define the local structure tensor.

(29)

3.4 Orientation and Tensors 17

Figure 3.1. Two various neighborhoods and their energy contribution in Fourier domain.

Illustration from G.H. Granlund, H. Knutsson, (1995). Signal Processing for Computer Vision[7].

3.4 Orientation and Tensors

Local orientation can be described by 2nd order tensor which should fulfill the following requirements:

• The closing requirement, • The uniform stretch, • The polar separability.

These requirements should be fulfilled in any type of orientation representation.

The Closing Requirement

Discontinuity should be avoided which means all pairs in space, vectors x and -x, should be mapped to the same point in the feature space.

T (x) = T (−x) (3.7)

where x is the vector and T is the tensor which is the map of x [10].

(30)

All orientation changes should be preserved after mapping (signals scale uni-formly). δ ˆT ∝ kδˆxk (3.8)

where x is the vector from the original space, and ˆT (ˆ stands for normalization)

is the mapping of the vector ˆx (normalized) [7].

The Polar Separability

The norm of the tensor (T) is independent of signal (x) direction since magnitude of the orientation x is independent of direction [10].

kT k = f (kxk) (3.9)

One of the mappings which preserves all of the above conditions and is suitable for signals of any dimensionality is defined as:

T ≡ AˆxˆxT (3.10)

where A is a constant (A > 0) and x is a vector pointing to the direction of the orientation [7].

3.5 Quadrature Filters and Estimation of

Orien-tation Tensor

In this section the process from quadrature filter to estimation orientation ten-sor will be reviewed which is based on Knutsson et al. representations in [7, 15]. The local orientation of the N-dimentsional signals can be represented by the orientation tensor as N × N real symmetric matrix [1]. The eigenvalues and corresponding eigenvectors of the tensor indicate the dominant direction of the signal energy. If there is not any dominant direction, then the local neighbour-hood will be represented by an isotropic tensor [1]. Several estimation methods of orientation tensor have been introduced over the past years which differ in some properties like orientation estimation accuracy [23]. Some of these methods are

gradient tensor, polynomial tensor, and quadrature tensor, which is the method

based on the principles of quadrature filters, and can be used to construct orien-tation tensor (for comparison of the various tensor estimation methods see [23]). The adaptive filter consists of low-pass filter (Flp) and high-pass filter (Fhp) in

which Flp+ Fhp = 1. It also eliminates the negative frequency component and

provides phase invariant magnitude [7]. The quadrature filter in Fourier domain is expressed as:

(31)

3.5 Quadrature Filters and Estimation of Orientation Tensor 19

where R(|u|) is radial band-pass function and Dk(ˆu) is directional function

which is zero in half plane [7]. The radial function is used to indicate the size of the structure which the filter should detect.

R(|u|) = exp( 4 B2_ln2ln

2_{(|u| /u}

0)) (3.12)

where B is the bandwidth and u0 is the center frequency [7].

(

Dk(¯u) = (ˆnk• ¯u)2 if ˆnk• ¯u > 0,

0 otherwise. (3.13)

where ˆnk is the direction of high-pass filters and ¯u is the frequency coordinate [7].

As mentioned above, by using quadrature filters, the negative frequency will be eliminated. This can be expressed by the Hermitian property of the Fourier transform of the real image I(¯x):

F {I(¯x)} = I(¯u)

I(¯u) = I∗(−¯u)

where u is the frequency vector and * denotes complex conjugate [4]. Also the energy contribution is even [4] i.e.,

|I(¯u)|2= |I(−¯u)|2= (Re {I(¯u)})2+ (Im {I(¯u)})2

The quadrature filter can be calculated by convolution of the filter with the signal (eq. 3.14). The response will be a complex value with real and imaginary parts (eq. 3.15).

q = f (¯x) ∗ h(¯x) (3.14)

q = qe+ iqo (3.15)

The real part represents the line structure, and the imaginary part represents the edge structure in the image, in the certain orientations [4]. The magnitude of the quadrature filter (eq. 3.16) estimates the phase invarient intensity and considering the phase, (eq. 3.17) it can distinguish edge or line structures [4] (figure 3.2).

|q| =pq2

e+ qo2 (3.16)

θ = arctan(qo/qe) (3.17)

There are two types of energy components which come from even and odd parts of the filters. The even and odd filters provide continuous representation in the spatial domain. The dominant energy contribution from (θ = 0, π) even filters introduces light and dark lines, and the dominant energy from (θ = +/ − π/2) odd filters corresponds to edges [4].

(32)

Figure 3.2. Local phase representation as edge and line detector. Illustration from M.

Andersson. Controllable Multi-dimensional Filters and Models in Low-level Computer Vision, (1992) [4].

3.6 Tensor Construction in 3D

In order to represent the orientation in the 3D space more than four filters are required. In order to achieve uniform and symmetrical distribution fashion at least 6 filters is required in which they can be oriented in the hemi-icosahedrons polygon [7]. The tensor representation of local orientations in the neighborhood is given by: T =   t1 t2 t3 t2 t4 t5 t3 t5 t6  = X k |qk| Mk (3.18) Mk= 5 4nˆknˆ T k − 1 4I (3.19)

where ˆnk are filter orientation vectors, I is identity matrix, and Mk is the tensor

associated with filter k [7]. The 6 filter directions are defined as following:

ˆ n1= c(a, 0, b)T nˆ2= c(−a, 0, b)T ˆ n3= c(b, a, 0)T nˆ4= c(b, −a, 0)T ˆ n5= c(0, b, a)T nˆ6= c(0, b, −a)T (3.20) where a = 2, b = (1 +√5), c = (10 + 2√5)−1/2 [7].

So by knowing these directions (eq. 3.20) the corresponding dual tensor can be obtained as following:

(33)

3.7 Interpretation of Tensor Orientation 21

Figure 3.3. An icosahedrons, the 6 vectors can point to its 12 vertices in symmetrical

fashion. Illustration from G.H. Granlund, H. Knutsson, (1995). Signal Processing for Computer Vision[7]. ˆ M1=   A 0 D 0 B 0 D 0 C   Mˆ2=   A 0 −D 0 B 0 −D 0 C   (3.21) ˆ M3=   C D 0 D A 0 0 0 B   Mˆ4=   C −D 0 −D A 0 0 0 B   (3.22) ˆ M5=   B 0 0 0 C D 0 D A   Mˆ6=   B 0 0 0 C −D 0 −D A   (3.23) A = 5 − √ 5 4(5 +√5) C = 5 + 2√5 2(5 +√5) D = √ 5 4 B = −1 4

Considering the elements of each ˆMk, it can be seen that these filter directions

are related to each other through simple mirroring operation [7].

3.7 Interpretation of Tensor Orientation

There can be complex neighborhoods when considering the 3D dataset. It can be grouped as two regions; isotropic and non-isotropic [7]. The simple neighborhood

(34)

can be represented by tensor rank one but for high dimensionality the higher tensor rank is needed [15]. Generally the tensor can be expressed by its eigenvalues and corresponding eigenvectors.

T =X

k

λkˆekeˆTk (3.24)

where λk are eigenvalues and ˆek are corresponding eigenvectors [7].

3.7.1 2D Representation

In 2D case, the tensor implies an ellipse model of the local energy distribution and the orientation is indicated by the largest eigenvalue which introduces dominant energy distribution. The local tensor relation in 2D is given by:

T = λ1ˆe1eˆT1 + λ2ˆe2eˆT2 (3.25)

where λ1≥ λ2 ≥ 0, and depending on the ratio between the first and the second eigenvalues, the tensor shape can be varied [7]. The local tensor representation can be related to isotropic neighborhood or simple neighborhood (figure 3.4).

Figure 3.4. Tensors with simple neighborhood, λ2 ≈ 0 (left). Approximately simple neighborhood, λ1, λ2 are the eigenvalues and ˆe1, ˆe2, are the corresponding eigenvectors

(middle). Isotropic neighborhood, no typical orientation in which λ1∼= λ2(right).

3.7.2 3D Representation

The orientation tensor in 3D can be expressed as following:

T = λ1eˆ1ˆeT1 + λ2ˆe2eˆT2 + λ3eˆ3eˆT3 (3.26)

where λ1≥ λ2≥ λ3≥ 0, are eigenvalues and ˆekare corresponding eigenvectors [7].

Considering the eigenvalues and corresponding eigenvectors the various structure types can be defined. The dominant direction of the local tensor and the local shape of the neighborhoods structure depends on its eigenvalues. Wiklund et al. [11] have visualized the characteristic of the local tensor as glyph T-flash (figure 3.5). He defined a visualization method to describe the combination of three

(35)

3.7 Interpretation of Tensor Orientation 23

particular cases of tensor T in 3D dimension. Depending on the complexity of the neighborhoods in 3D space, there are three ranks of tensor representation [7].

• Plane case (T1), Tesor Rank 1: This indicates a neighborhood with planar-like prototype and its energy is distributed along a line in the Fourier domain.

• Line case (T2), Tensor Rank 2: This indicates a neighborhood with line-like prototype and its energy is distributed on a plane in the Fourier domain.

• Isotropic case (T3), Tensor Rank 3: No typical orientation, e.g. in the case of noise.

In general tensor representation can be expressed as the combination of above three cases [7]. T = T1+ T2+ T3 where T1= (λ1− λ2)ê1êT1, T2= (λ2− λ3)(ê1eˆT1 + ê2êT2), T3= λ3(ê1êT1 + ê2eˆT2 + ê3eˆT3),

Figure 3.5. Left: T-flash tensor glyph for λ1= 1, λ2= 0.5, λ3= 0.25. Right: an almost

isotropic tensor, λ1 = 0.67, λ2 = 0.63, λ3 = 0.59. The isotropic case is represented in

green sphere (its radius corresponds to the magnitude of the isotropic part), the planar case is represented in red spear through the poles (the length of the spear is proportional to the magnitude of T1), and the plane case is represented in yellow expanded sphere orthogonal to e3 direction. see [11] for details. Illustration from J. Wiklund et al.

(36)

3.8 Eigenvalue Calculation

In order to calculate the eigenvalues of the local tensor, the Cardano’s Formula is used [5]. So the solution of the cubic equation (eq. 3.27) can be used to calculate the eigenvalues from the equation 3.28.

ax3+ bx2+ cx + d = 0 (3.27)

The quadratic term can be removed with substitution x = z - a/3 as following:

z3_{+ pz + q = 0} where p = b −a 2 3 and q = c + (2a3− 9ab) 27 , ∆ = p 3 27+ q2 4 , u =r −q3 2 + √ ∆, v =r −q3 2 − √ ∆, z1= u + v, z2= −(u + v) 2 + (u − v) 2 i √ 3, z3= −(u + v) 2 − (u − v) 2 i √ 3.

When there are 3 real roots then ∆ < 0, which causes some complications in the complex cubic roots for the above equation. The eigenvalues are positive and real since the matrix is symmetric and positive semi-definite [1]. The following scalling has been made for calculation:

z =r −4p 3 y, and 4y 3_{− 3y =} 3q p q −4p 3

By using the identity cos(3α) = 4 cos3_{(α)−3 cos(α), [1] and substitution y = cos(α)} the following equation is obtained:

cos(3α) = 3q p q −4p 3 , β =r −4p 3 ,

(37)

3.9 Evaluation and Error Estimation 25 α = 1/3 arccos(3q pβ), z1= β cos(α), z2= β cos(α − 2π 3 ), z3= β cos(α + 2π 3 ).

In order to compute the eigenvalues algebraically, the equation (3.28) should be solved. Also the above discussion leads to removing the isotropic part of the tensor which is equivalent to removing the quadratic term from the characteristic polynomial [1]. |T − λI| = 0 (3.28) T0 = T −tr(T )I 3 =   t1 t2 t3 t2 t4 t5 t3 t5 t6   (3.29) p = t1t4+ t1t6+ t4t6− (t22+ t23+ t25), q = t1t25+ t4t23+ t6t22− 2t2t3t5− t1t4t6

So we have p and q, by using the above relations z1, z2, z3 will be calculated which are corresponding to λ1, λ2, λ3 respectively, since 0 ≤ α ≤ π/3.

3.9 Evaluation and Error Estimation

In order to evaluate the constructed tensor, a synthetic 3D test volume is used. It is a spherical shells shape with 64 × 64 × 64 dimension and wide frequency range (figure 3.6). Therefore all tensor directions can be evaluated. In order to measure the orientation estimation error, a method which is based on the RMS (root mean squared) angular measurement is used [1, 7]. The evaluation method is based on finding the difference (angular error) between the correct tensor orientation (T1= ˆ

xˆxT) and the estimated one ( T = ê1eˆT1, the estimated orintation corresponds to largest eigenvalue). ∆Φ = arcsin( v u u t( 1 2L L X l=1 xˆˆx T − ê1êT1 2 ) (3.30)

where ˆx is the vector which indicates the expected direction (correct direction)

and ˆe1 is the eigenvector corresponding to the largest eigenvalue of the estimated orientation, L is the number of points in selected landmarks, and ∆Φ is the RMS angular error [1, 7]. The RMS error indicates the euclidean distance between the

(38)

points located in the correct direction and the ones in the estimated direction. In order to avoid the edge effect, the sum of the points within the volume test data (7 < R < 32) are measured. The angular error estimation is done for both 2D and 3D, with kernel size 5 × 5 × 5 for 3D and 15 × 15 × 15 for 2D and various input patterns. The results are given in table below.

Table 3.1. Results from RMS angular measurement

SNR ∆Φ2D ∆Φ3D

∞ 0.55◦ _0.47◦

10 9.06◦ 5.80◦ 5 14.79◦ 10.55◦ 0 23.44◦ 17.59◦

Figure 3.6. volume visualization of the spherically symmetric synthetic test data which

is used for evaluation.

From the above result the local orientation estimation has better accuracy which indicates the processing of 3D dataset with 3D operator would be desirable. This result is obtained before averaging of the tensor output and it would give better estimation if the evaluation is implemented after averaging.

(39)

3.10 Control Tensor Calculation 27

3.10 Control Tensor Calculation

Control tensor has the same eigen system as local tensor T as following:

C = γ1eˆ1eˆT1 + γ2eˆ2ˆeT2 + γ3ˆe3eˆT3 (3.31)

where ˆekare eigenvectors of Tlpand γ1, γ2, andγ3are mappings of the local tensor’s eigenvalues [7]. In order to do mapping of the orientations of local tensor, two functions, m − f unction and µ − f unction are used. These two functions provide some parameters which can be used to adjust the suitable mapping to get the desired control tensor which introduces the adaptive filtering.

3.10.1 Control of High-pass Content

The m − f unction controls the overall high-pass content of the adaptive filter [1]. It declares how much band-pass filter content should be used in the result of the image. m(x, σ; α, β, j) = _xβ xα+β_{+ σ}β 1/j (3.32)

where x is the position-dependent variable and it is measured by the magnitude of local tensor (|T |), σ is intended to be an estimation of the noise, α defines overshoot of the function which can be adjusted to improve the visibility of small structures,

β is the steepness of the function, and j represents the function iteration [7]. Figure

3.7 shows how the m−f unction has been effected by different α values. The largest

Figure 3.7. m − f unction with β = 4, σ = 0.2, j = 1, α = 0, 0.2, 0.4, 0.6, 1. (bottom

to top respectively), x ∈ [0 1].

(40)

of the adaptive filter which is obtained using the m − f unction [7].

γ1= m (kTlpk , σ; α, β, j) (3.33)

where kTlpk is the magnetiude of tensor, and "lp" stands for low-pass filtering.

Tlp= glp∗ T

where T is the local tensor, glp is low-pass Gaussian filter, and * represents the

convolution [7]. In order to map the local tensor to the control tensor, tensor aver-aging is needed. The averaver-aging of the local tensor suppresses the rapid changing in the representation of the local tensors [7]. This is done by convolution of Gaussian filter (size 7 × 7 × 7 and standard deviation 1) with the tensors. So this averaging yields smoother changes in the local tensors [15].

(41)

3.10 Control Tensor Calculation 29

3.10.2 Control of High-pass Content Distribution

By µ−f unction the orientation variation in the band-pass region will be controlled and it controls the shape of adaptive filters [4 ,7]. The output of this function shows the anisotropic energy distribution in the local structure consideration. For example high frequency component, which is represented as white color (high energy) in γ1output, will be represented as black color in µ2output (figure 3.10) [7]. However, the isotropic regions in dataset will be represented as high valued (white color) by µ2, which can introduce noisy regions. So the shape of the tensor field and the shape of the energy distribution will be controlled by this function.

µ(x; α, β, j) =

_{(x(1 − α))}β

(x(1 − α))β_{+ (α(1 − x))}β

1/j

(3.34)

where the parameters α and β control the shape of the energy distribution (fig-ure 3.9). The ratio between eigenvalues of the local tensor, λn/λn−1 control the

mapping process and the shape coefficients, µn, are as following:

µn= µ(

λn

λn−1

; α, β, j) (3.35)

where α defines a value for ratio _λn−1λn in which the magnitude of µ − f unction in set to half, β controls the sharpness of the µ − f unction, and j stands for iteration [7].

(42)

Figure 3.9. µ − f unction with α = 0.4, j = 1, β = 1, 2, 4, 6, 12. (left to right

respectively), x ∈ [0 1], and in bottom µ − f unction with β = 1, j = 1, α = 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 0.9, 0.95.(from top to bottom), x ∈ [0 1].

(43)

3.11 Optimization of Quadrature Filters 31

Figure 3.10. The γ1 output (left), which shows the smooth changes on the edges from bright area to dark (σ = 0.2, α = 0.1, β = 4, j = 1). The µ2 output (right), which has

smooth changes on the edges (α = 0.45, β = 3, j = 1).

3.11 Optimization of Quadrature Filters

Following the discussion in chapter 2, the optimization of quadrature filters is done manually. Various parameters like: center frequency, bandwidth, filters size in spatial domain can be adjusted considering the desired structures in the volume dataset. The center frequency and the bandwidth of the filter can be estimated by the minimum (Lmin) and maximum (Lmax) line-width of the desired structures

in the image as following:

B = ln (Lmax/Lmin) [octaves]

u0= √

uLminuLmax [Hz]

where u0 is the center frequency, uLmin = π/Lmax, and uLmax = π/Lmin [7]. By

calculating the wave length corresponding to the center frequency, the suitable spatial size of the filters can be estimated. The size of the filters in spatial domain should be 1.5-2 times larger than the calculated wave lenght [7].

λ0= 2π

u0

[pixels]

3.12 Computing The Enhancement Filters

Now by having the γ1from equation 3.33 and calculation of µ2and µ3by equation 3.35, the eigen system of control tensor can be constructed.

γn = γ1

N

Y

n=2

(44)

where N here is 3 for 3D tensor estimation [7]. So all unknown eigenvalues of the control tensor (eq. 3.31) can be calculated and the eigenvectors are the same as eigenvectors of the local tensor. Generally, the enhancement filter consists of one isotropic low-pass (LP) filter and six anisotropic band-pass (BP) filters (for 3D case) in the same directions as the quadrature filters. Finally, the enhancement filters can be obtained by adding the LP filter and weighted BP filters as following:

Fenh= FLP +

X

k

WkFBP (3.37)

where Wk is calculated as following:

Wk = C • Mk (3.38)

where C is the control tensor matrix 3 × 3 which contains the eigenvalues of the control tensor and Mk are measured in equations (3.21 - 3.23) [7].

3.13 Implementation and Evaluation

The performance of the 3D adaptive filter is evaluated by both the 3D synthetic dataset and the clinical volume dataset. In order to verify the results from adaptive filtering, a white Gaussian noise with various signal to noise ratio (SNR) is used (eq. 3.39). In this section the results are presented. The 3D synthetic dataset is onion shape with size (64 × 64 × 64) (figure 3.6). The implemented filter has size 5 × 5 × 5, bandwidth 2[octaves], and center frequency of 2π/5.

SN R = 20 log( α σimage

(1 − α) σnoise

) (3.39)

where α is constant and σ is the standard deviation [7]. The 8 slices of original synthetic data and the one corrupted by Gaussian noise, with SNR = 0, is shown in figure 3.11. The result from 3D adaptive filtering, is shown in figure 3.12 which indicates how the high frequency contents and their energy are preserved while the isotropic high frequncy content (noise) is reduced significantly.

(45)

3.13 Implementation and Evaluation 33

Figure 3.11. The volume test data used for evaluation, 8 slices (top), The corrupted

(46)

Figure 3.12. Result from 3D adaptive filtering on the volume test data, 8 slices.

Figure 3.13. 3D visualization of synthetic data with additive Gaussian noise, SNR =

(47)

3.13 Implementation and Evaluation 35

Also to evaluate the 3D adaptive filtering on the clinical dataset, a T1-weighted MRA image of brain with dimension of 240×240×140 is used. The volume data has been corrupted by white Gaussian noise with SNR = 0, 5, and 10 using equation 3.39. In order to filter the current volume clinical data, a filter with size 9 × 9 × 9 in spatial domain, bandwidth 3[octaves], and center frequency of π/3 is used. The axial slice (z = 90) of the corrupted volume data and the enhancement results from adaptive filtering are shown in figure 3.14.

(48)

Figure 3.14. Axial view of the T1-weighted MRA image, slice 90, from top to bottom

original image, corrupted image with SNR = 0, 5, and 10. From left to right noisy image and the result from the 3D adaptive filtering with size 9 × 9 × 9 in spatial domain.

(49)

3.14 Discussion 37

In order to check the controlled tensor and the estimated orientation for high frequency contents, the response of the controlled tensor to the noise as the only input was evaluated. The shape of the enhancement filter after summing up 5 iteration is shown in figure 3.15. First, the control tensor is measured from the synthetic test data (figure 3.6) as the input which produces weight function for the enhancement filters (equation 3.38). Then, the new input data to the filters is constructed as Gaussian noise with mean value of 0 and with the same size as the synthetic data (a cube of the noise 64 × 64 × 64). The filter output is weighted by the constructed weight function. The expected output is round (followed by the weight function), symmetric linear, and should follow the filter orientation (figure 3.15). Its frequency is the same as the center frequency of the enhancement filters.

Figure 3.15. Top from left to rigth: The weight function which is constructed by the

control tensor and the result of the enhancement filters to the input which is the white Gaussian noise after 20 iterations. Bottom: the sum of enhancement band-pass filters after 5 iterations in frequency domain, with filter size 7 × 7 × 7 in spatail domain.

3.14 Discussion

In this chapter a method for designing 3D adaptive filtering and enhancement of 3D volume dataset, by using quadrature filters and considering local orientation

(50)

tensor, is discussed. The magnitude of quadrature filter is used to construct the orientation tensor, which represents the direction and energy components of se-lected neighborhood voxels. Then the controlling process is performed to weight each band-pass filter (they have the same direction as the quadrature filters) and set the different coefficients to get adaptive filtering in multidimensional scale. The main advantages of this method are; computationally efficiency, straight forward parameterization, being extremely robust, and it gives a phase-invariant result [4,7]. The mapping functions are used to construct the control tensor and finally construct the weight function for band-pass filters. The mapping function param-eters are initiated the same for all input data which may not be optimized for all input dataset. One solusion is to construct an interface which can control the mapping property (see the thesis work of Karin Adolfsson [20]). Also there are other methods to construct the local orientation tensor like: the gradient tensor which is proposed by Bigun and Granlund [23], the orientation tensor based on a local polynomial model [1] which is derived by Farnebäck, etc. (for more detail about various orientation estimation method see [21])

(51)

Chapter 4

Segmentation

4.1 Introduction

Magnetic Resonance Imaging (MRI) is a non-invasive diagnostic imaging method without radiation, and pain. It provides excellent contrast for soft tissue imag-ing, e.g. brain tissue, muscular tissue, abdomen, etc. The main parts of an MRI machine are; super conductive big Magnet, RF Amplifier, Gradient Coil, Trans-mitter/Receiver Coils [22].

The dominant proton in the body is hydrogen nuclei. When a patient is placed in the strong magnetic field, the magnetic moment of the hydrogens in his body align to the magnetic field (B0). The resonance frequency of the hydrogen is 42.58 [MHz/T] which is called gyromagnetic ratio. The hydrogens precess along the magnetic field with special frequency called larmor frequency (w0 = γ B0). The transmitter coil transmits the radio frequency pulse (B1) with larmor frequency to excite protons in the selected volume. The protons absorb the energy and then generat the MR signal which will be recieved by the reciver coils. The gradient coil generates three directional (GX, GY, GZ) magnetic fields (10-50 [mT/m]). GZ

and RF bandwidth will contribute together in slice selection. The GX and GY

change the B0 field spatially inside the selected region. This variation introduces a spatial encoding which provides enough information for an image reconstruction [22].

Magnetic Resonance Angiography (MRA) is one of the MRI imaging protocols which is used for vascular imaging. The MRA volume contains high intensity where vascular structures exist and low intensity elsewhere. Various methods are introduced to enhance the vascular structures in the MRA images. There are common problem issues related to MRA image analysis like; discontinuity, small vessels, variation of contrast, and bifurcations.

Vascular extraction and enhancement have numerous applications in medical diagnostics, surgical planning, and visualization, etc. In order to visualize the vascular structure, a set of pre-processing steps are required to suppress non-vessel structures. One part of the pre-processing discussed in the previous section is the enhancement of volume dataset, and the other part is vessel structure extraction

(52)

40 Segmentation

which will be discussed in this chapter.

4.2 Tubular Structure Extraction

The idea is to extract the voxels which belong to the tubular structure prototype in the volume dataset, considering the local structure of the neighborhoods. In order to classify a tubular structure within a neighborhood, the characteristics of possible objects should be estimated. By considering tensor eigenvalues [7], the tubular structure prototype can be classified by following probability:

PT S1=

λ2− λ3

λ1

, (4.1)

Actually the above ratio describes the probability of existence of any tubular struc-ture in the volume dataset. It implies that if there are any tubular strucstruc-tures in the dataset then it should return high value since λ3 is low. In order to classify the volume data by this ratio, the smooth mapping function like the one used to estimate γ1 (m − f unction, equation 3.34), is applied. Instead of using σ as the threshold for noise components, which can be seen as isotropic structures in volume dataset, the new threshold, δ, is defined. The new threshold δ removes the undesired structures in the volume dataset.

m (PT S1, δ, α, β, j), (4.2)

Figure 4.1. m − f unction with β = 3, α = 0.1, j = 1, δ = 0.1, 0.2, 0.4, 0.6 (top to

bottom respectively), λ2−λ3

(53)

4.3 Related Work 41

4.3 Related Work

Some related works have been done to develop a strategy to extract tubular struc-tures and most of them are based on Frangi’s et al. method. In this section, first the three previous strategies are briefly reviewed and then the output of the proposed method is compared with the other methods in the next section. Frangi defines a method based on multi-dimensional line filtering using Gaussian differ-ential [18]. He used the eigenvalues of the 3 × 3 Hessian matrix to measure and evaluate the geometric structures in the volume dataset [18]. So the desired out-put is obtained by considering the eigenvalues of the Hessian matrix [12]. The Hessian matrix is measured by the gradient of the volume data H(x) = ∇ V(x). Considering its eigenvalues, the various structures can be estimated in the volume dataset. T = ( 0 if λ2> 0 or λ3> 0, (1 − exp((RA2) 2α2 )) exp( −(R2 B) 2β2 )(1 − exp( −S2 2c2 )) if else. (4.3) where RA= |λ1| p|λ2λ3| , RB = |λ2| |λ3| , S = q λ2 1+ λ22+ λ23,

The orders of eigenvalues [12] are as following:

|λ1| ≤ |λ2| ≤ |λ3| ≤ 0

RAindicates the blob-like structures, RBcorresponds to the structures between

plane-like and line-like, and S is used to remove the random noise effects [12, 18]. The parameters α, β and c are used to control the sensitivity of the filter to each of the structure indicators (RA, RB, and S respectively). It is assumed the

background is dark and blood vessels are bright in the volume dataset, which is the normal case in MRA and CTA images. The computed eigenvalues |λ1| ≤ |λ2| ≤ |λ3| and the corresponding eigenvectors can be used to detect vessels (figure 4.2). The eigenvector corresponds to λ1indicates the direction along the vessel and the others eigenvectors are perpendicular to the direction of the vessel [18].

By the Frangie’s filtering method, more of the non-tubular structures are sup-pressed and it works well if there is no noise because it is highly sensitive to noise [9, 18].

(54)

42 Segmentation

Figure 4.2. Second-order 3D variation ellipsoid shape space. The classification is based

on eigenvalues of the Hessian matrix. The bright string prototype indicates vessel varia-tions [12, 18]. Illustration from Lin Qingfen,(2003), Enhancement, Detection, and Visu-alization of 3D Volume Data [18].

Enhancement and Visualization of VascularStructures in MRA Images Using Local Structure

Institutionen för medicinsk teknik

Department of Biomedical Engineering

Examensarbete

Enhancement and Visualization of Vascular

Structures in MRA Images Using Local Structure

Enhancement and Visualization of Vascular

Structures in MRA Images Using Local Structure

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Abstract

Acknowledgments

Contents

List of Figures

List of Tables

Chapter 1

Introduction

1.1

Motivation

1.2

Objective

1.3

Background

1.4

Outline

Chapter 2

Filter Optimization

2.1

Introduction

2.2

Representation Spaces

2.3

Filter Optimization Components

2.4

Distance Measurement

2.5

Computational Complexity

2.6

Discussion

Chapter 3

Local Orientation

Estimation in 3D Space

3.1

Introduction

3.2

Tensor Basics

3.3

Orientation Estimation

3.4

Orientation and Tensors

3.5

Quadrature Filters and Estimation of

Orien-tation Tensor

3.6

Tensor Construction in 3D

3.7

Interpretation of Tensor Orientation

3.7.1

2D Representation

3.7.2

3D Representation

3.8

Eigenvalue Calculation

3.9

Evaluation and Error Estimation

3.10

Control Tensor Calculation

3.10.1

Control of High-pass Content

3.10.2

Control of High-pass Content Distribution

3.11

Optimization of Quadrature Filters

3.12

Computing The Enhancement Filters

3.13

Implementation and Evaluation

3.14

Discussion