Ultrasound Imaging Methods for Breast Cancer Detection

(1)

Ultrasound Imaging Methods

for

(2)

(3)

Ultrasound Imaging Methods

for

Breast Cancer Detection

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op donderdag 13 November 2014 om 12.30 uur

door

Neslihan ¨OZMEN

Master of Science in Scientific Computing geboren te Ankara, Turkey.

(4)

en copromotor:

Dr. K.W.A. van Dongen

Samenstelling promotiecommissie:

Rector Magnificus voorzitter

Prof. dr. ir. A. Gisolf Technische Universiteit Delft, promotor Dr. K.W.A. van Dongen Technische Universiteit Delft, copromotor Prof. dr. ir. P.M. van den Berg Technische Universiteit Delft

Prof. dr. H.P. Urbach Technische Universiteit Delft Prof. dr. ing. G. Schmitz Ruhr-Universit¨at Bochum Prof. dr. ir. W. Steenbergen Universiteit Twente

Dr. N. Ruiter Karlsruhe Institute of Technology

ISBN 978-94-6186-392-8

Copyright c 2014 by N. ¨Ozmen

Laboratory of Acoustical Wavefield Imaging, Faculty of Applied Sciences, Delft University of Technology, Delft, The Netherlands.

This research is supported by the Dutch Technology Foundation STW. Typeset by the author with the LA

TEX Documentation System. Cover design: Nalan Liv and Neslihan ¨Ozmen

(5)

The more you struggle to live, the less you live. Give up the notion that you must be sure of what you are doing. Instead, surrender to what is real within you, for that alone is sure... You are above everything distressing. Baruch Spinoza

Your hand can seize today, but not tomorrow; and thoughts of your tomorrow are nothing but desire. Do not waste this breath, if your heart is not crazy, since ‘the rest of your life’ will not last forever. Omar Khayyam

(6)

(7)

Chapter

1

Introduction

Any intelligent fool can make things bigger and more complex. It takes a touch of genius - and a lot of courage - to move in the opposite direction.

E.F. Schumacher (attributed to Albert Einstein) Cancer is a disease in which cells grow out of control. Cancer cells may also invade other tissues and spread to other parts of the body through the blood and lymphatic system. In our body, when cells become old or damaged, they are replaced with new cells. However, sometimes this normal process can be affected by mutations. The DNA of a cell can be damaged and changed, such that it does not die and may produce a new type of cell. These extra cells may form a tissue mass which is called a tumor. Figure 1.1 illustrates an example of controlled and uncontrolled growth of cells. The top row shows the normal cell division, and the normal balance between cell growth and cell loss. Old or damaged cells are normally self-destructed by a mechanism called cell-suicide or apoptosis. The bottom row shows the cancer cell division. During the cancer development, the balance is disrupted and this disruption causes a gradual increase in cell division.

There are many different known tumors that can affect humans, but not all tumors are cancerous. They can be divided in two categories: benign and malignant. Benign tumors are not cancerous, and do not

(12)

(13)

1.2 BREAST IMAGING SYSTEMS WITH ULTRASOUND 3 chance of patients (Jemal et al., 2011). However, in early stages it is too small to detect, and there are no signs and symptoms. When it is diagnosed, it may be too late to treat the cancer, because it can rapidly spread to the other parts of the body. Figure 1.2 shows an example that describes the stages of breast cancer. In the first stage (top left), the tumor is detected less than two cm inside the breast and has not spread outside the breast. In the second stage (top right), the size of the tumor increases, and it has spread to the axillary lymph nodes. In the third stage (bottom left), the tumor has spread to the tissues near the breast (skin, chest-wall or under the arm). In the fourth stage (bottom right), it has spread other organs (e.g., lungs, liver or brain) which is called metastasis (Novartis, 2014).

Scanning the breast periodically is crucial for early diagnosis. Self-examination, X-ray mammography, MRI and ultrasound are some of the common techniques used for breast cancer detection (Kolb et al., 2002). Currently, X-ray mammography is the common modality clinically used for breast imaging. However, it can miss cancers in young women with dense breasts as both the healthy fibrous and glandular tissue, and the cancerous lumps show up white on mammograms (Harvard, 2011). X-ray mammography is also potentially harmful because of the ionizing radiation involved. Besides, some patients find it uncomfortable and painful. Another modality used for breast imaging is MRI, which avoids many of those problems. Since it is expensive and requires a long examination time, it is not useful for scanning large groups. Ultrasound has the potential to be an ideal screening modality for breast cancer detection, because it is relatively inexpensive and requires no ionizing radiation. In addition, ultrasound has the ability to discover small cancers in women whose mammograms are normal, but who have dense breast tissue (Gordon and Goldenberg, 1995; Duric et al., 2007; Schreiman et al., 1984). Based on this fact, at present, breast imaging with ultrasound is gaining interest.

1.2 Breast Imaging Systems with Ultrasound

Clinical studies demonstrate that ultrasound can find more and smaller cancers in dense breast women, where mammography might miss them (Gordon and Goldenberg, 1995; Ying et al., 2012). However, current hand-held ultrasound scanning systems are impractical, because

(14)

Figure 1.2: Stages of breast cancer and how it spreads and forms tumors elsewhere in the body through metastasis (Novartis, 2014).

they are operator-dependent, time-consuming and result in too many false positives. Therefore, many researchers are working on automated whole-breast ultrasound scanning systems.

A well-known scanning system is the Automated Whole-Breast Ultrasound (AWBU) of SonoCine, which was introduced in 2010 by the Sutter North Bay Women’s Health Center, to improve hand-held breast cancer examination (Kelly and Richwald, 2011). AWBU mimics hand-held whole-breast ultrasound, and allows the radiologist to read the images quickly at a convenient time. The Food and Drug Administration (FDA) has stated that this technology may be used as an adjunctive examination to mammography but not as a replacement.

Another system is being developed by Philips, and currently, feasibility of the automated 3D whole-breast ultrasound imaging system is tested in a clinical pilot study. The system can also be merged with Philips’s diffuse

(15)

1.2 BREAST IMAGING SYSTEMS WITH ULTRASOUND 5

Figure 1.3: SoftVue scanning system (Duric et al., 2014).

optical mammography system, to combine ultrasound and optical imaging. The prototype uses three 3D transducers in a fixed geometry to image the whole breast. During the scan procedure there is no interaction with the operator. The added value of combining ultrasound imaging with diffuse optical tomography (DOT) is discussed in Leproux et al. (2010).

Currently, water bath scanners are gaining interest, as they have the potential to give high-quality volumetric images of the breast. This approach is based on ultrasound tomography principles and uses a fixed geometry, so that an automated mechanism can scan the breast and give reproducible images. Since breasts vary in size and shape, water is used as an acoustically coupling medium between the transducers and the skin. In addition, the water bath approach ensures minimal deformation of the breast, while providing a more accurate and complete image of the tissue characteristics. We will discuss in detail three examples of water bath scanners: SoftVue, CVUS and USCT.

• SoftVue: SoftVue, a commercial breast imaging device built in Karmanos Cancer Institute (KCI) together with Delphinus Medical Technologies (Duric et al., 2007, 2013a,b). The SoftVue scanner is pictured in Figure 1.3. It uses a ring-shaped transducer system that

(16)

Figure 1.4: CVUS scanning system (on the left) and transducer arrays (on the right) (Wiskin et al., 2013).

simultaneously generates images from both the transmitted and the reflected signals. The system has 2048 transducer elements with 512 transmitter and 512 receiver channels. The operating frequency is 3 MHz, and the image resolution is 2.5 mm × 0.3 mm × 0.3 mm. Typically 30 to 40 tomograms (image slices) with a slice thickness of 2.5 mm are constructed for each scan. The maximum breast diameter that it can handle is 22 cm. Each examination takes about one minute and produces images for the radiologist in less than 15 minutes. The applied tomographic reconstruction algorithms yield sound-speed, attenuation, and reflection images. Their future research aims to validate SoftVue’s ability to differentiate benign masses from cancer.

• CVUS : Computerized volumetric ultrasound (CVUS) is a clinical prototype built in University of Utah (Andre et al., 2012a; Wiskin et al., 2013), and is shown in Figure 1.4. The system contains a transmitter/receiver pair, and three reflection arrays allowing the system to collect both transmission and reflection data. Hence, two types of imaging algorithms can be used to reconstruct the breast from data collected from CVUS. For transmission mode imaging, data are gathered by transmitting from a specially designed transmitter which transmits an uniquely designed chirp that contains energy in the range of 0.3 MHz − 2 MHz, and is centered at 1.1 MHz. After passing through the water and the breast, the signal is finally collected by the receiver array on the other side. The receiver array contains 1536 elements distributed over eight vertical rows. For reflection mode imaging, three reflection

(17)

1.3 AIM OF THE STUDY 7 transducers are used. These arrays (5 MHz center, 80% bandwidth) focus at 2.5 mm, 4 mm, and 7.5 mm. For their clinical study, they reported that 3D reflection and transmission images were produced for each patient successfully, and each scan took about 8 minutes. Moreover, the patients experienced little or no discomfort during the scanning process.

• USCT : Karlsruhe Institute of Technology (KIT) is currently developing an ultrasound computed tomography (USCT) system, which generates promising 3D images with high-spatial resolution (0.2 mm × 0.2 mm) (Ruiter et al., 2006, 2011). Their current design has 157 transducer arrays, where each array contains four emitters and nine receivers yielding 628 emitters and 1413 receivers in total. The arrays are distributed on the inner surface of a semi-ellipsoidal aperture of 21.5 cm height and with 34.0 cm external diameter. Figure 1.5 shows the 3D USCT II system with new optimized patient bed (Ergonomic patient bed for better immersion depth of the breast). Temperature control and movement mechanics are mounted together with the aperture in the patient bed. Each transmitter generates spherical wave fronts with 2.4 MHz center frequency (approximately 50% bandwidth). The complete sensor system can be rotated and translated, which creates further virtual positions of the ultrasound transducers. In addition, data acquisition time is less than two minutes per breast. It can also provide speed of sound, attenuation, and reflection images.

1.3 Aim of the Study

One of the most important aspects of the different ultrasound imaging systems is the applied imaging method – in the end it is the resulting image which counts, and not how it has been obtained. Consequently, the primary objective of this study is to investigate the applicability of six ultrasound imaging methods for breast cancer detection. The examined methods – Inverse Radon, refraction corrected time-of-flight tomography (Eikonal), synthetic aperture focusing technique, back-propagation, Born inversion and contrast source inversion – are all tested on computed measurement data and show an increase in numerical complexity. Additional attention is paid to contrast source inversion. This full-wave

(18)

Figure 1.5: 3D USCT II scanning system (on the left) and transducer array housing (on the right) (Ruiter et al., 2013).

nonlinear inversion method is tested on both 2D and 3D problems, while in addition, it is investigated how the reconstruction efficiency is affected by reducing either the number of sources or receivers.

To understand the outcome of the above study, and in particular the effect the different acoustic medium parameters (volume density of mass, compressibility and absorption) and the application of the Born approximation have on the wave field, numerical studies are performed. For these studies, the forward problem is solved multiple times using a 3D computer model built from an MRI scan of a real cancerous breast.

Finally, we believe that some preliminary work is needed before testing the methods using real data. In particular, we will study and suggest some steps which may help to design a real scanning system.

1.4 Thesis Outline

The research presented in the thesis is focused on the breast cancer detection problem using ultrasound imaging techniques. We solve the forward problem and model acoustic wave propagation using both full-wave and Born approximation methods. Considering an example scanning setup, we compute the acoustic wave field propagation in a breast that is embedded in water. We investigate the acoustic medium

(19)

1.4 THESIS OUTLINE 9 parameters of different breast tissues and in particular their effect on wave propagation. We analyse imaging and inversion methods by implementing them and looking at the results obtained with synthetic data generated using our forward model. In addition, we study how to shift from synthetic to real measurements. So, we investigate the necessary aspects that should be analyzed before starting to work with the real data sets (e.g., sampling criterion, radiation pattern and transfer functions of the transducers).

Chapter 2 gives the theory of acoustic wave equations starting from the equations of motion (Newton’s law) and deformation (Hooks’s law). The integral equation formulation is also derived in this chapter, which is the key concept in this thesis. Chapter 3 introduces the forward problem (defining a numerical method for modeling wavefield propagation), and investigates different solutions. The chapter first discusses the incident field computation, followed by the description of the full-wave method and the Born approximation. In addition, an analytical solution is given to validate the proposed methods. Chapter 4 discusses various imaging and inversion methods. Inverse Radon, refraction corrected time-of-flight tomography (Eikonal), synthetic aperture focusing technique, Born inversion and contrast source inversion are the methods investigated in this context. These imaging methods are tested using the synthetic data generated by simulations performed in Chapter 3. Comparison of these methods are presented as well. Chapter 5 evaluates some intermediate steps for tieing in synthetic and real measurements to assess feasibility and to design and optimize breast ultrasound scanning systems. The steps presented in this chapter are sampling criterion, perfectly matched layers (PMLs) and system characterizations (i.e., measuring the sound speed of the background medium, calculation of the incident field and the transfer function of the transducers). It also provides suggestions for future experiments. Finally, Chapter 6 closes the thesis with a summary and discussion of the findings of this study, concluding with recommendations for future research.

(20)

(21)

Chapter

2

Theory

I believe that we do not know anything for certain, but everything probably.

Christiaan Huygens

Though ray theory is often used for acoustical imaging, our study is based on wave theory. Ray theory is a high-frequency approximation of wave theory, which means that the medium should vary smoothly compared to the wavelength (Gisolf and Verschuur, 2010). Particularly, the wavelength of the propagating field should be much shorter than any dimension in the volume of interest. If this is the case, the ray will follow a minimum-time path and ray-tracing may be used to describe the wave propagation and to predict the arrival times of reflected and refracted fields. Since it neglects physical phenomena like diffraction and scattering, we will focus on full-wave propagation. Wave theory explaining the acoustic wave propagation in inhomogeneous media will be presented in this chapter. In addition, a frequency-domain integral-equation formulation will be derived, which describes the fundamental model of our study.

(22)

2.1 Acoustic Wave Equations

To derive the acoustic wave equations we first introduce the equation of motion and deformation. The equation of motion is given as

∇p(x, t) + ˙Φ(x, t) = f(x, t), (2.1)

where p(x, t) is the acoustic pressure wavefield, ˙Φ(x, t) is the mass-flow density rate, and f (x, t) is the volume source density of volume force at spatial location x = (x, y, z) and time t. The deformation equation is given as

∇ · v(x, t) − ˙Θ(x, t) = q(x, t), (2.2) where v(x, t) is the particle velocity wavefield, ˙Θ(x, t) is the induced part of the cubic dilatation rate, and q(x, t) is the volume source density of injection rate. The symbol ∇ indicates the nabla operator containing the spatial derivatives and · indicates an inner product.

It is known that, the mass-flow density rate ˙Φ(x, t) depends on the particle velocity v(x, t), and the induced cubic dilatation rate ˙Θ(x, t) depends on the acoustic pressure p(x, t) (Fokkema and van den Berg, 1993; Gisolf and Verschuur, 2010). For inhomogeneous media the relations can be given as

˙

Φ(x, t) = ρ(x)Dtv(x, t), (2.3)

˙

Θ(x, t) = −κ(x)Dt[m(x, t) ∗tp(x, t)] . (2.4)

Here, the volume density of mass ρ(x) and compressibility κ(x) are the constitutive coefficients. Attenuation, i.e., the irreversible transformation of acoustic energy into heat, is included via a casual compliance relaxation function m(x, t) (Huijssen, 2008; Demi et al., 2011). The symbol ∗t indicates a convolution with respect to time, and Dt is the

material time derivative:

Dt=

∂

∂t+ v · ∇. (2.5)

Within the low-velocity approximation, the material time derivative equals Dt=

∂

(23)

2.1 ACOUSTIC WAVE EQUATIONS 13 Using this approximation and substituting equations (2.3) and (2.4) in equations (2.1) and (2.2), acoustic wavefield equations are obtained as

     ∇p(x, t) + ρ(x)_∂t∂v(x, t) = f (x, t), ∇ · v(x, t) + κ(x)_∂t∂ [m(x, t) ∗tp(x, t)] = q(x, t). (2.7)

Applying a temporal Fourier transformation to equation (2.7) yields (

∇ˆp(x) + iωρ(x)ˆv(x) = ˆf(x),

∇ · ˆv(x) + iωκ(x) ˆm(x)ˆp(x) = ˆq(x), (2.8) where ω is the angular frequency. The caret symbol ˆ is used to show that a particular field or quantity is defined in the frequency domain, and is obtained using a temporal Fourier transform, which equals

ˆ u(x) = ∞ Z −∞ u(x, t) exp(−iωt)dt. (2.9)

We can rewrite the wavefield equations shown in equation (2.8), as (

∇ˆp(x) + iωρ0v(x) = ˆˆ f(x) − iω ρ(x) − ρ0 ˆv(x),

∇ · ˆv(x) + iωκ0mˆ0p(x) = ˆˆ q(x) − iω κ(x) ˆm(x) − κ0mˆ0 ˆp(x),

(2.10) where the properties ρ0, κ0 and ˆm0 refer to a homogeneous background

medium. Taking the spatial derivative of the first line, multiplying the second line with iωρ0, and subtracting the second equation from the first

one yields ∇2_{p(x) − ˆγ}ˆ ₀2_{p(x) = −}ˆ hiωρ0q(x) − ∇ · ˆˆ f(x) i − ˆ γ₀2₋ ρ0 ρ(x)ˆγ 2_(x) ˆ p(x) − ∇ ·hiω ρ(x) − ρ0 ˆv(x) i , (2.11)

where ∇2 = ∇·∇ is the Laplace operator, and ˆγ is the complex propagation coefficient, i.e., ˆ γ(x) = iωpκ(x)ρ(x) ˆm(x) = iω c(x)p ˆm(x), (2.12) ˆ γ0 = iωpκ0ρ0mˆ0 = iω c0 p ˆ m0, (2.13)

(24)

with the speed of sound, c0= 1 √_κ 0ρ0 . (2.14)

Assuming that ˆf(x) = 0 in most biomedical applications, and using relation (2.8) for the particle velocity

iωˆv_{(x) = −} 1

ρ(x)∇ˆp(x), (2.15)

the acoustic wave equation can be rewritten as

∇2_{p(x) − ˆγ}ˆ ₀2_{p(x) = − ˆ}ˆ Spr(x) − ˆScs(x), (2.16)

with

ˆ

Spr(x) = iωρ0q(x),ˆ (2.17)

representing the real external source and ˆ Scs(x) = ˆχγ(x)ˆp(x) + ∇ · h χρ(x)∇ˆp(x) i , (2.18)

representing the contrast source accounting for the inhomogeneity of the medium. The contrast functions ˆχγ(x) and χρ(x) which depend on the

medium parameters are equal to ˆ χγ(x) = ˆγ02− ρ0 ρ(x)γˆ 2_(x), _(2.19) χρ(x) = ρ0− ρ(x) ρ(x) . (2.20)

To account for the power law attenuation (Wells, 1975; Duck, 1990) typically observed in biomedical tissue, the complex propagation coefficient ˆγ(x) is defined as

ˆ

γ(x) = ˆα(x) + i ˆβ(x). (2.21)

The attenuation coefficient ˆα(x) and phase coefficient ˆβ(x) are defined as (Huijssen, 2008; Demi et al., 2011)

   ˆ α(x) = a(x) (2π)−b(x)_|ω|b(x), ˆ β(x) = ω c(x) + a(x) (2π) −b(x) _tanhπ 2b(x) i ω |ω|b(x)−1, (2.22)

(25)

2.2 INTEGRAL EQUATION FORMULATION 15 where a(x) and b(x) are medium dependent attenuation parameters. Note that ˆm(x) and the constant ˆm0 follow from equations (2.12) and (2.21) as

p ˆm(x) = c(x) iω ˆγ(x), (2.23) p ˆ m0 = c0 iωγˆ0. (2.24)

2.2 Integral Equation Formulation

Equation (2.16) is known as the Helmholtz equation for homogeneous media and can be solved by convolving the source term, ˆSpr(x) + ˆScs(x),

with the Green’s function. The Green’s function is the solution of the Helmholtz equation, and describes the impulse response of the background medium. The field that is generated from a source and propagates in absence of any contrast is referred to as the incident pressure wavefield. The total (or actual) pressure field ˆp(x) can be seen as a summation of the incident pressure field ˆpinc and the scattered pressure field ˆpsct(x),

ˆ

p(x) = ˆpinc(x) + ˆpsct(x). (2.25) Similar to the domain integral equation formulated by Fokkema and van den Berg (1993), we obtain the domain integral equation

ˆ p(x) = ˆpinc(x) + Z x′∈D ˆ G(x − x′) ˆχγ(x′)ˆp(x′)dV + Z x′∈D ˆ G(x − x′)n_∇′_·χρ(x′)∇′p(xˆ ′) o dV, (2.26)

for all x ∈ R3, where D is a three-dimensional (3D) spatial domain, and ˆ

G(x) is the Green’s function given by ˆ

G(x) = e

−ˆγ0|x|

4π|x| . (2.27)

Furthermore, the incident pressure field ˆpinc(x) is computed by convolving the Green’s function with the source ˆq(x).

(26)

The Green’s function given in equation (2.27) has a singularity at |x| = 0. In our study, problems associated with this singularity are avoided by using a weak form of the Green’s function. By computing the spherical mean, the weak form of the Green’s function (Zwamborn and van den Berg, 1992) is obtained as (see Appendix A)

[ ˆG](x) =                          1 − (1 + ˆγ0∆x₂ )e− ˆγ0 ∆x 2 1 6γˆ0 2_{π (∆x)}3 if |x| = 0, e− ˆγ0|x| 1 3π ˆγ02(∆x) 2 |x| cosh ˆγ0∆x 2 − 1 ˆ γ0∆x₂ sinh ˆγ0 ∆x 2 # if |x| ≥ ∆x, (2.28)

where ∆x = ∆y = ∆z indicates the grid size in the discrete model. The forward and inverse scattering problems can be defined using equation (2.26). The forward (or direct) problem refers to the situation where the unknown total pressure field ˆp(x) is obtained for the known incident pressure field ˆpinc(x) and the known contrast functions ˆχγ(x)

and χρ(x). Meanwhile, the inverse problem refers to the situation where

the unknown contrast functions ˆχγ(x) and χρ(x) are obtained from the

incident pressure field ˆpinc(x), which is known everywhere, and the total pressure field ˆp(x), which is typically measured and hence known at a limited number of locations outside the domain containing the contrasts.

Born approximation

The integral equation (2.26) contains a spatial convolution of the Green’s function with a contrast source, which is a function of the contrast and the total pressure field. For the inverse problem, both quantities are unknown. In situations where the contrasts are small and the scattering is weak, the unknown total field ˆp(x) inside the integrand of equation (2.26) may be approximated by the known incident field ˆpinc(x), yielding the following

(27)

2.2 INTEGRAL EQUATION FORMULATION 17 integral representation ˆ p(x) = ˆpinc(x) + Z x′∈D ˆ G(x − x′) ˆχγ(x′)ˆpinc(x′)dV + Z x′∈D ˆ G(x − x′)n_∇′_·χρ(x′)∇′pˆinc(x′) o dV, (2.29)

for all x ∈ R3. The unknown total field can now be computed directly. This approximation neglects multiple scattering inside the domain of interest. Since no iteration process is required, it is highly efficient with respect to computational time.

(28)

(29)

Chapter

3

Forward Scattering Problem

Science cannot solve the ultimate mystery of nature. And that is because, in the last analysis, we ourselves are a part of the mystery that we are trying to solve.

Max Planck

The forward or direct problem is defined as computing the total wavefield in the domain of interest, in which the acoustic medium parameters and incident wavefield are known. Research in developing new imaging methods which takes advantage of advanced scanning systems, requires computer modelling and simulations as an appropriate, inexpensive, and flexible approach for generating synthetic measurement data. It is also practical for optimizing the design of new scanning systems, to model a specific transducer geometry, to investigate the required number of sources and receivers, or to test if a particular system design meets its predefined requirements. In addition, such a data set can be used to investigate if a particular frequency range in combination with an appropriate imaging or inversion method is suitable to image the breast with sufficient accuracy.

In this study, a full-wave method based on an integral equation formulation is used to model the acoustic wavefield propagation in breasts accurately. This methods includes all phenomena typically

(30)

observed in real measurement data, such as multiple scattering, refraction and diffraction. To reduce the complexity of the full-wave method,

approximations can be made. One of the most well-known

approximations is the Born approximation, which is also implemented for comparison. In this chapter, both the Born and the full-wave methods are explained in detail. In the results section, the accuracy of the full-wave method is confirmed using the analytical solution for a spherical contrast. In addition, the full-wave and the Born approximation methods are tested and compared using a 3D breast model.

3.1 Full-Wave Method

The presented full-wave method allows to model the acoustic pressure wavefield in large spatial domains with arbitrary shaped contrasts. It can be used to compute the propagation and scattering of acoustic wavefields in breasts, and to generate noise-free synthetic measurement data, which still possess all phenomena typically observed in real measurements. The method is based on an integral equation formulation and allows for inhomogeneities in compressibility, mass density and attenuation. Using integral equation (2.26) derived in Chapter 2, the total pressure field can be computed using a breast model for which the acoustic medium parameters are exactly known. Then, equation (2.26) can be reformulated as ˆ pinc(x) = ˆ_{p(x) −} Z x′∈D ˆ G(x − x′) ˆχγ(x′)ˆp(x′)dV − Z x′∈D ˆ G(x − x′)n_∇′_·χρ(x′)∇′p(xˆ ′) o dV, (3.1)

for x ∈ D, where ˆp(x) is the total pressure field, ˆpinc(x) is the incident pressure field, ˆ_{G(x − x}′) is the Green’s function, and ˆχγ(x) and χρ(x)

are the contrast functions. The corresponding integral equation can be formulated as a linear system and solved using iterative methods. In the next section, the iterative methods used in this study are explained in detail.

(31)

3.1 FULL-WAVE METHOD 21

3.1.1 Iterative solution methods

To provide a numerical solution, equation (3.1) can be formulated in operator notation as f = L[u], (3.2) where L_{[u] = u(x) −} Z x′∈D ˆ G(x − x′) ˆχγ(x′)u(x′)dV − Z x′∈D ˆ G(x − x′)n_∇′_·χρ(x′)∇′u(x′) o dV, (3.3)

f is the known incident pressure field ˆpinc(x), u is the unknown total pressure field ˆp(x), and L is the integral operator containing the known contrast functions ˆχγ(x) and χρ(x), and the Green’s function ˆG(x).

Iterative methods are often used to solve this type of equations by minimizing an error functional (van den Berg, 1991). The approximate solution at the nth iteration step can be given as

u0 = 0

un= un−1+ αndn, n≥ 1, (3.4)

where αn is the step size and dn is the update direction. For any updated

solution un, differing from the exact solution u, the residual rn is defined

as

rn= f − L[un], (3.5)

and the normalized error En as

En= krnkD

kfkD

, (3.6)

with the properties En= 0 if un= u, and En= 1 if un= 0, and where k·k2D

represents the L2-norm of a vector. The error here is used as a measure

for the accuracy attained in the iterative scheme. The inner product is defined as

hu(x), v(x)iD=

Z

x∈D

u∗(x)v(x)dx. (3.7)

The symbol∗ represents the complex conjugate of the function. Typically, the iterative process is stopped once the maximum allowable number of

(32)

iterations is reached, or once the normalized error goes below a predefined value ǫ, referred to as the error criterion.

Different choices on how to compute αn and dn yield different

iterative schemes, e.g., Neumann iteration, conjugate gradient, biconjugate gradient methods (Fletcher, 1976; Jacobs, 1986; Kleinman and van den Berg, 1991; Sarkar, 1987). Conjugate gradient (CG) is one of the most well-known methods, which always converges when the operator is self-adjoint and positive definite. However, the integral operator L is not always self-adjoint, and a symmetrization procedure has to be performed. This can be achieved by applying the adjoint of the operator L referred as L†. The adjoint of the operator L is defined via the inner product, i.e.,

hf, L[u]iD= hL†[f], uiD. (3.8)

Assuming the contrast functions, ˆχγ and χρ, reach the boundary of the

domain and using the Gauss theorem, we obtain the following expression for the adjoint operator

L†_{[f] = f (x) − ˆ}χ∗_γ(x) Z x′∈D ˆ G∗(x′_{− x)f(x}′)dV − ∇ ·  χ∗_ρ_{(x) ∇} Z x′∈D ˆ G∗(x′_{− x)f(x}′)dV  . (3.9) Applying the adjoint operator in each iteration increases the computational time. Here, we use a version of the conjugate-gradient algorithm which includes the adjoint operator. The detailed steps of the algorithm are given in Table 3.1.

We also apply the BiCGSTAB (van der Vorst, 1992) which is an efficient variation of the CG method. Since it avoids the implementation of the adjoint operator, it converges faster than the CG scheme. However, when amongst others strong attenuation is taken into account, the BiCGSTAB method sometimes encounteres difficulties with convergence, whereas CG remains convergent. The applied BiCGSTAB algorithm is given in Table 3.2 in detail.

3.2 Born Approximation

To reduce the computational costs needed to model the acoustic wavefield, the forward problem can also be solved within the Born approximation, in

(33)

3.2 BORN APPROXIMATION 23 Table 3.1: The conjugate gradient scheme (L† _{is the adjoint operator)}

CG u0 = 0 r0 = f − L[u0] η1= 0 for n = 1,2,... gn= L†[rn−1] ηn= kgnk2D/ kgn−1k 2 D dn= gn+ ηndn−1 αn= kgnk2D/ kLdnk2D un= un−1+ αndn rn= f − L[un] En= krnkD/ kfkD if En< ǫ or n > nmax stop end

which the unknown total field within the integral is replaced by the known incident field. The resulting integral equation is given in Chapter 2, see equation (2.29), and equals

ˆ p(x) = ˆpinc(x) + Z x′∈D ˆ G(x − x′) ˆχγ(x′)ˆpinc(x′)dV + Z x′∈D ˆ G(x − x′)n_∇′_·χρ(x′)∇′pˆinc(x′) o dV. (3.10)

This approximation provides a computationally efficient solution of the forward problem, but it neglects phenomena such as multiple scattering

(34)

Table 3.2: The BiCGSTAB scheme BiCGSTAB u0 = 0 r0 = f − L[u0] z= r0 ρ0 = α = w0= 1 v0 = p0= 0 for n = 1,2,... ρn= hz, rn−1iD β = ρnα/ρn−1wn−1 pn= rn−1+ β(pn−1− wn−1vn−1) vn= L[pn] α = ρn/hz, vniD s= rn−1− αvn t= L[s] wn= ht, siD/ht, tiD un= un−1+ αpn+ swn rn= s − twn En= krnkD/ kfkD if En< ǫ or n > nmax stop end

(35)

3.3 INCIDENT FIELD COMPUTATION 25 inside the domain of interest. In addition, it neglects phase shifts for waves traveling through an object with a different speed of sound. To investigate the applicability of the Born approximation for our application, the outcome of the full-wave method is compared with the results of the Born approximation. Note that the first step of the conjugate-gradient iteration algorithm yields similar results as applying the Born approximation.

3.3 Incident Field Computation

The incident pressure field can be computed in various ways, and it is generally computed by spatially convolving the primary sources, see equation (2.17) and (2.18), with the Green’s function. In this section, the various incident fields used for the simulations are described.

3.3.1 Plane wave

To validate the accuracy of our forward solution method, the analytical solution for a plane wave scattering at a spherical object is used. Typically, the incident plane wave propagates in the y-direction and equals

ˆ

pinc(x) = ˆse−ˆγ0y_, _(3.11)

where ˆs is the source signature and y is the direction of propagation in the spatial domain. In Appendix B.2, the analytical solution for a plane wave scattering at an acoustically penetrable sphere is derived.

3.3.2 Point source

In general, an injection rate point source is used when computing the wave propagation in a breast model. The resulting incident field is computed by spatially convolving the source term ˆq(x) with the Green’s function ˆG(x),

ˆ pinc(x) = iωρ0 Z x′∈D ˆ G(x − x′)ˆq(x′)dV, (3.12) where ρ0 is the volume density of the mass of the background medium.

(36)

3.3.3 Planar transducer

On the other hand, considering the real scanning systems mentioned previously in Chapter 1, in practice it is more accurate to compute the pressure field generated by vibrating surfaces. The Rayleigh integral method can be used as a flexible approach for any arbitrary geometry and vibrating planar surface. It describes the response to a vibrating planar surface via the pressure field along that surface. The wavefield generated by an infinite source plane S0 can be computed via an infinite

integral over S0 (i.e., over x and y in the z = zS plane). This integral is

known in the seismic literature as a Rayleigh II integral. So, the incident field can be computed in an observation point xA = (xA, yA, zA) via the

Rayleigh integral as ˆ pinc(xA) = zA 2π ∞ Z −∞ ∞ Z −∞ ˆ P (x, y, zS) (1 + ˆγ0∆r)exp(−ˆγ0∆r) (∆r)3 dxdy, (3.13) with ∆r =p(x − xA)2+ (y − yA)2+ (zS− zA)2. (3.14)

Although the Rayleigh integral method is very convenient for any shape of planar transducer, it is computationally very expensive. More efficient methods may be used to reduce the computational load. For predefined transducer geometries faster methods can also be applied to compute the incident pressure fields in the region of interest. Some of these methods are spatial impulse response (Lockwood and Willette, 1973),

fast-nearfield method (McGough, 2004) and time-space

decomposition (Kelly and McGough, 2006). Moreover, some software packages such as FOCUS (McGough and Urban, 2012) or Field II (Jensen, 1999), are based on these methods. Further details can also be found in (Alles, 2012; van Nijhuis, 2012), since these methods were not investigated in detail in this thesis.

3.4 Results

This section presents the simulation results by applying the previously discussed full-wave method. First, we validate the method using an analytical solution for a plane wave scattering at a spherical object. Secondly, we run several simulations using a breast model with

(37)

3.4 RESULTS 27 Table 3.3: Acoustic medium parameters for the analytic solution

Tissue c ρ κ a b

[m/s] [kg/m3] [1/Pa] [dB/m MHzb] [-]

Water 1524 993 4.34e-10 0.2 2.00

Sphere 1485 1020 4.45e-10 44.0 1.01

inhomogeneities in all three medium parameters (i.e., compressibility, mass density and attenuation), separately and combined. Finally, we compare the outcome of the second test with the results obtained with the Born approximation.

3.4.1 Validation with an analytical solution

To test the accuracy of the full-wave method, we use a configuration for which we know the analytical solution: a plane wave scattering at a spherical object. The derivation of the analytical solution can be found in Appendix B. The acoustically penetrable sphere with radius r = 2.5 mm is positioned in the center of a volume of 10 mm × 10 mm × 10 mm consisting of water. The acoustic medium parameters of the sphere are chosen similar to those typically encountered in biomedical applications (Table 3.3). The sphere is illuminated by a Gaussian plane wave that is propagating in the y-direction with a center frequency f0 = 1 MHz and a

bandwidth of 50%. The time span of the simulation is set to 20 µs, with the step size ∆t = 0.16 µs. The computational domain is discretized with a step size ∆x = ∆y = ∆z = 0.125 mm. This corresponds to ten points per wavelength at the center frequency which is necessary for accurate numerical integration. The iterative process is stopped when the normalized error is En≤ ǫ with ǫ = 10−6.

First, we compare the frequency domain results for the center frequency f0. Figure 3.1 presents the real part of the incident and total

pressure fields along the y-axis, obtained with both the analytical and the full-wave method. Both solutions are virtually identical. In addition, the results show that for both solutions the pressure fields crossing the interface of the sphere (indicated by the solid lines) are continuous, and therefore satisfy the acoustic boundary condition.

(38)

−5 −4 −3 −2 −1 0 1 2 3 4 5 −1 −0.5 0 0.5 1 y [mm] Re[p(x=0,y,z=0)] pinc ptot_AS ptot FW

Figure 3.1: Comparison of the frequency domain wavefields obtained with the full-wave method (FW) and the analytical solution (AS). The black solid line is the incident field used in both cases. The red dotted line indicates the AS, the blue dots show the solution obtained with the FW, and the black lines at y = −2.5 mm and y = 2.5 mm indicate the boundaries of the sphere.

transforms. Snapshots of the incident, total and scattered fields at time t = 8 µs in the plane z = 0 mm are presented in Figure 3.2.

The normalized error between the analytical solution PAS(x, t) and

the full-wave method PFW(x, t) is calculated for the incident, total and

scattered fields, as Error1(x, t) = 20 log10 |P FW(x, t) − PAS(x, t)| |PAS(x, t)| , (3.15)

where | · | indicates the absolute value. As can be seen in Figure 3.2, the differences between both solutions are small. For instance, the normalized error in the total field is around −40 dB and occur mostly on the boundaries, which is due to the spatial discretization of the sphere and numerical rounding errors.

3.4.2 Simulations using a breast model

For this study, we use a breast model (see Figure 3.3), built from an MRI scan of a real cancerous breast (Bakker et al., 2009), which is submerged

in water and has acoustic medium parameters as shown in

Table 3.4 (Szabo, 2004; Goss et al., 1980; d’Astous and Foster, 1986; Duck, 1990; Lide, 1999). Those medium parameters may vary according

(39)

3.4 RESULTS 29

Figure 3.2: Snapshots of the incident (top row), scattered (middle row) and total (bottom row) pressure fields in dB relative to 1 Pa at time t = 8 µs in the plane z = 0 mm. The first column shows the results obtained with the full-wave method, and the second column shows the error between the analytical solution and the full-wave method. The white circle shows the contour of the sphere located in the middle of the domain.

(40)

Table 3.4: Acoustic medium parameters for the breast (Szabo, 2004; Goss et al., 1980; d’Astous and Foster, 1986; Duck, 1990; Lide, 1999) (∗ _{assumed values)}

Tissue c ρ κ a b [m/s] [kg/m3] [1/Pa] [dB/m MHzb] [-] Water 1510 995 4.41e-10 0.2 2.00 Muscle 1580 1041 3.85e-10 57.0 1.01∗ Fat 1430 928 5.27e-10 15.8 1.70 Skin 1537 1200 3.53e-10 104.0 1.01∗ Breast 1510 1020 4.30e-10 75.0 1.50 Tumor 1550 1000 4.16e-10 57.0 1.30

to patients and type of tumor. It may also depend on the temperature and frequency of the measurements. Unfortunately, this information is not always available.

During the simulations, the time span considered equals 0.3 ms, and is discretized with a time step ∆t = 1 µs. A point source, located at the point (xs, ys, zs) = (0 mm, −47 mm, 31 mm), generates a Gaussian

modulated field with a center frequency f0 = 0.25 MHz and 50%

bandwidth. The computational domain is set as a volume of 100 mm × 100 mm × 62 mm, which is discretized with a grid size ∆x = ∆y = ∆z = 0.8 mm. To satisfy the requirement for accurate numerical integration, a spatial sampling of seven grid points per center frequency wavelength is chosen. The stopping criterion set by ǫ = 10−6. We choose these parameters to ease the computation and to demonstrate the proof of concepts.

To investigate the effect of each acoustic medium parameter separately, the forward problem is solved four times; one simulation where all medium parameters are set to their appropriate values, and three simulations where only one of the three medium parameters has its appropriate value whereas the remaining parameters are set to values corresponding to the background medium (i.e., water). Consequently, the following cases are considered:

• If we assume that there is only a compressibility contrast we neglect density and attenuation variations (ρ(x) = ρ0 and ˆm(x) = ˆm0), then

(41)

3.4 RESULTS 31

Figure 3.3: Tissue types for the 3D breast model built from an MRI scan of a real cancerous breast (Bakker et al., 2009). The bottom image displays a cross-section of the breast at z = 31 mm. The medium parameters corresponding to the different tissues are given in Table 3.4.

the contrast functions become

χρ(x) = 0, (3.16) ˆ χγ(x) = ˆγ02− ˆγ2(x) = ω2mˆ0 1 c2_(x) − 1 c2 0 . (3.17)

• If we assume that there is only a density contrast we neglect compressibility and attenuation variations (κ(x) = κ0 and

(42)

ˆ

m(x) = ˆm0), then the contrast functions become

χρ(x) =

ρ0− ρ(x)

ρ(x) , (3.18)

ˆ

χγ(x) = 0. (3.19)

• If we assume that there is only attenuation contrast we neglect compressibility and density variations (κ(x) = κ0 and ρ(x) = ρ0),

then the contrast functions become

χρ(x) = 0, (3.20) ˆ χγ(x) = ˆγ₀2− ˆγ2(x) = ω 2 c2₀ ( ˆ_{m(x) − ˆ}m0) . (3.21)

Figure 3.4 shows the scattered pressure fields for the different contrasts in the plane z = 31 mm, which are evaluated at three different frequencies (0.125 MHz, 0.25 MHz and 0.5 MHz). The scattering caused by the different inhomogeneities is clearly visible. The images show that the amount of scattering increases for increasing frequency. In addition, it is shown that the inhomogeneities in compressibility and attenuation have a significant larger effect on the field, compared to the effect of the inhomogeneities in volume density of mass. Finally, the ripples visible inside the breast are spaced roughly one wavelength apart. They are caused by interference of the field inside the breast, and can only be modeled via a full-wave method which allows for multiple scattering.

Time domain results are obtained by applying the inverse fourier transform. Time frames at time t = 0.07 ms and in the plane z = 31 mm of the incident, scattered and total pressure fields are presented in Figure 3.5. All fields are normalized according to the incident field. The total field includes all wave phenomena such as scattering, reflection, refraction and attenuation.

3.4.3 Comparison of the Born approximation and the full-wave method

The effect of applying the Born approximation is investigated using the same set of contrast functions as defined in Section 3.4.2. However, in this case time domain wavefields are compared. Time frames of the total

(43)

3.4 RESULTS 33 x [mm] Compressibility y [mm] 0.125 MHz −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] 0.25 MHz −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] 0.5 MHz −40 −20 0 20 40 −40 −20 0 20 40 x [mm] Density y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] Attenuation y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] All Inhomogenities y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 −85 −80 −75 −70 −65 −60 −55 −50 −45 −40 −35

Figure 3.4: Scattered pressure field in the plane z = 31 mm in dB

normalized with respect to the maximum absolute value of the incident field of the entire spatial computational domain. The different columns show results for different frequencies (from left to right 0.125 MHz, 0.25 MHz and 0.5 MHz), and the different rows show results for different contrasts (from top to bottom compressibility, density, attenuation, and all three combined).

(44)

x [mm] y [mm] Incident field −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] Scattered field −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] Total field −40 −20 0 20 40 −40 −20 0 20 40 −90 −85 −80 −75 −70 −65 −60 −55 −50 −45

Figure 3.5: Snapshots of the incident (left), scattered (middle) and total (right) pressure fields in dB at time t = 0.07 ms and in the plane z = 31 mm. They are all normalized according to maximum absolute value of the incident pressure field. pressure fields modeled with the Born approximation and with the full-wave method at t = 0.07 ms and z = 31 mm, for different contrasts are shown in Figure 3.6. The figure shows differences between both sets of solutions, which are attributed to, among others, the multiple scattering. It also shows that the spatially varying compressibility and attenuation causes significant scattering for both the Born approximation and the full-wave solution.

A-scans of incident and total fields are calculated at the observation point (xo, yo, zo) = (0 mm, 47 mm, 31 mm), and presented in Figure 3.7

and Figure 3.8 for the Born approximation and the full-wave method, respectively. All A-scans are normalized with respect to the maximum of the incident field and interpolated by zero-padding in the frequency domain. Finally, Figure 3.9 presents the A-scans, at the same observation point, of the incident and total fields (with all inhomogeneities present) for results obtained with the Born approximation and using the full-wave method. Application of the Born approximation again yields a significant increase of pressure amplitude. A normalized error is computed to quantify the difference between the Born approximation and the full-wave method including all inhomogeneities at the plane z = 31 mm, i.e.,

Error2 =

kPFW(x, y, t) − PBorn(x, y, t)k_{(x,y)∈D,t∈R}

kPinc_{(x, y, t)k}

(x,y)∈D,t∈R

(45)

3.4 RESULTS 35 x [mm] y [mm] Full−wave −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] Compressibility y [mm] Born −40 −20 0 20 40 −40 −20 0 20 40 x [mm] Density y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] Attenuation y [mm] −40 −20 0 20 40 −40 −20 0 20 40 x [mm] All Inhomogenities y [mm] −40 −20 0 20 40 −40 −20 0 20 40 −85 −80 −75 −70 −65 −60 −55 −50 −45 −40

Figure 3.6: Snapshots of the total wavefield obtained with the Born

approximation (left) and using the full-wave method (right), at time t = 0.07 ms and in the plane z = 31 mm. Both rows show the total fields obtained for different contrast functions (from top to bottom: compressibility, density, attenuation, and all three combined). All fields are in dB, and are normalized with respect to maximum absolute value of the incident field computed in the entire spatial computational domain for the compressibility contrast case.

(46)

The normalized error between the two methods yields Error2 = 0.38. The

unexpected amplitude increase in the A-scans and the normalized error indicate that the Born approximation can yield spurious results.

3.5 Summary

To meet the demand for a method that can model the acoustic wavefield propagation in breasts accurately, we presented a 3D full-wave method based on a frequency domain integral equation formulation. The method allows for continuous inhomogeneities in compressibility, volume density of mass and attenuation. In order to solve this integral equation, we use conjugate gradient based iterative methods which update the unknown total field in each iteration by minimizing a cost functional.

The accuracy of the full-wave method is confirmed using the analytical expression for a plane wave scattered at a spherical penetrable object. Results obtained via the full-wave method and the analytical expression are in close agreement with each other. Moreover, the normalized error is everywhere well below −20 dB. The minor differences observed are mainly a consequence of the spatial discretization and the stopping criterion En,

i.e., the normalized error used for the iterative scheme.

A breast model obtained from an MRI scan of a cancerous breast is used for several simulations. Inhomogeneities belonging to the different medium parameters are investigated separately and combined, and results obtained with the full-wave method are compared with results obtained with the Born approximation. The results indicate that the Born approximation leads to errors especially in the presence of compressibility contrasts, which in some cases even cause a significant increase in field amplitude, where we expect a decrease. A second disadvantage of applying the Born approximation is that phase shifts caused by variations in the speed of sound are not handled correctly. Moreover, from a theoretical point of view it is clear that within the presented framework the scattered field obtained with the Born approximation only travel with a speed of sound of the background medium. Most of the imaging methods neglect attenuation, due to the computational complexity. Nevertheless, spatially dependent attenuation and mass density have effects on the pressure field, and considering these effects during imaging may be useful for accuracy.

(47)

3.5 SUMMARY 37 0.065 0.07 0.075 0.08 0.085 0.09 0.095 0.1 0.105 0.11 0.115 −2 −1 0 1 2 Born t [ms] p(t) Incident Compressibilty Density Attenuation 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.5 1 1.5 2 2.5 f [MHz] |p(f)|

Figure 3.7: The Born approximation results of the A-scans at the point

(xo, yo, zo) = (0 mm, 47 mm, 31 mm), normalized with respect to the incident field. The upper figure shows the time domain results and the lower figure shows the frequency domain results. The solid black line is the incident field. The dashed blue, dashed-dotted green, and dotted red lines are the total fields for compressibility, density, and attenuation contrasts, respectively.

0.065 0.07 0.075 0.08 0.085 0.09 0.095 0.1 0.105 0.11 0.115 −1 −0.5 0 0.5 1 Full−wave t [ms] p(t) Incident Compressibilty Density Attenuation 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.2 0.4 0.6 0.8 1 f [MHz] |p(f)|

(48)

0.065 0.07 0.075 0.08 0.085 0.09 0.095 0.1 0.105 0.11 0.115 −2 −1 0 1 2 All inhomogenities t [ms] p(t) Incident Born Full−wave 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.5 1 1.5 2 2.5 f [MHz] |p(f)|

Figure 3.9: Results obtained using the Born approximation and the full-wave method with all inhomogeneities. The A-scans are at the point (xo, yo, zo) = (0 mm, 47 mm, 31 mm), normalized with respect to the incident field. The upper figure shows the time domain results and the lower figure shows the frequency domain results.

In future, the full-wave method presented can contribute to the development of new imaging methods, because it allows one to generate noise-free synthetic data, which still possess all phenomena typically observed in real measurement data, such as multiple scattering, refraction and diffraction. Such a data set can be used to investigate if a particular frequency range in combination with a particular imaging or inversion method is suitable for imaging the breast with sufficient accuracy. In addition, it can be a useful tool during the design of new scanning systems. It may be used to model a particular transducer geometry, to investigate the required number of sources/receivers, or to test if a particular system design meets its predefined requirements.

(49)

Chapter

4

Imaging and Inversion

Aerodynamically, the bumble bee shouldn’t be able to fly, but the bumble bee doesn’t know it so it goes on flying anyway.

Mary Kay Ash

Ultrasound is widely used in medical imaging such as fetal, cardiac, orthopedic, abdominal, genitourinal, and breast imaging (Duric et al., 2005; Hill and Bamber, 2004). Specifically, it can be utilized for breast cancer screening of high-risk women, young women or women with dense breasts. It has the potential to provide accurate characterization of the breast tissue. Sonography or B-mode imaging is the most accepted technique based on pulse-echo principles. It uses reflection data to produce grey-scale images. However, acoustic wave propagation in heterogeneous media is not simply described by reflections only. It is demonstrated that transmission data for ultrasound tomography provides higher quality images (Filipik, 2009; Hill and Bamber, 2004). The feasibility of ultrasound computer tomography in breast imaging was first recognized by Greenleef (Greenleaf et al., 1975, 1978) and Glover (Glover and Sharp, 1977; Glover, 1977). Later, researchers in the field of medical imaging attempted to build automated scanning systems, which aim to contribute to the detection of breast cancer (Andre et al., 2012b; Duric et al., 2007; Gemmeke and Ruiter, 2007). Some of currently available

(50)

whole breast ultrasound scanning systems are introduced in Chapter 1. The general principles of ultrasound tomographic imaging techniques are explained in this chapter. Several image reconstruction techniques are discussed in detail and tested using a two-dimensional (2D) synthetic example. In principle, reconstructing an image of the acoustic medium parameters from the acoustic field is a nonlinear inverse problem. Therefore, the full-wave inverse problem based on the integral equation methods will also be discussed. Finally, a full-wave inversion method will be investigated in some detail and tested with a three-dimensional (3D) synthetic example. Additionally, we assume that there is no density contrast (ρ(x) = ρ), and no dispersion in the data ( ˆm(x) = ˆm = 1), for all imaging methods discussed in this chapter.

4.1 Imaging with Ultrasound Tomography

Although ultrasound imaging is useful for breast cancer detection, conventional handheld ultrasound techniques have major drawbacks. First of all, imaging capabilities are limited when using the straight-ray approximation, or assuming a 2D geometry which acquires only back-scattered data. Furthermore, it is difficult to generate reproducible images because of the operator dependency of the examination. Application of tomographic imaging techniques have the potential to improve the image quality and may allow ultrasound to become a useful tool in breast cancer diagnosis (Duric et al., 2005; Mueller et al., 1979). Tomographic image reconstruction methods use reflection and transmission algorithms. In general, these methods can be divided into three groups: ray-based tomography, diffraction tomography and full-wave inversion.

4.1.1 Ray-based tomography

Initial attempts show that it is reasonable to assume straight-ray propagation of the acoustic wavefield in soft tissues which only shows a small variation in the speed of sound (Lavarello and Hesford, 2013). In particular breast imaging is an important application of the straight-ray approximation. Since it neglects refraction, diffraction and scattering phenomena this method only produces low-resolution images. Reflection tomography is formulated by (Norton and Linzer, 1979b,a) using data collected in pulse-echo mode (Duric et al., 2007; Ruiter et al., 2006).

(51)

4.1 IMAGING WITH ULTRASOUND TOMOGRAPHY 41 Later, in order to increase the image quality and to include refraction effects in the ray-based methods, refraction-corrected tomography was introduced (Li et al., 2010; Norton and Linzer, 1982). The method employs the Eikonal equation, which is based on a high frequency approximation. Thus, it assumes the wavelength to be small compared to the size of the scatterer. In addition, it requires to measure a large number of ray paths using different source/receiver pairs. Another difficulty when reconstructing refraction-corrected tomography images is that it requires a first estimate of the spatial distribution of the index of refraction.

4.1.2 Diffraction tomography

Since the size of the scattering structures inside the breast are in the order of one wavelength, diffraction correction techniques need to be considered for sub-wavelength resolution. Diffraction tomography, usually referred to linearized methods, first introduced by (Mueller, 1980), can produce sound speed and attenuation images. It basically solves wave equations, rather than geometrical ray models, using linearizion techniques such as the Born and Rytov approximations (Devaney, 1981, 1982b, 1983).

Diffraction tomography systems can provide promising image quality, but they mainly suffer from instability. So, if the stopping criterion is not defined appropriately, the algorithm may diverge, and reconstructed image can be blurred (e.g., with wrong amplitude, missing / wrong-sized structure). In addition, these methods assume that there are only small variations in acoustic parameters, and therefore only working under weakly scattering conditions (Slaney et al., 1984).

4.1.3 Full-wave inversion

Nonlinear inversion techniques which use the full-wave equation as a starting point, can be used to account for effects such as refraction, diffraction and multiple scattering. Inversion using the wave theory is described as computing the pressure wavefield and the medium parameters in a domain surrounding the breast (Wiskin et al., 2012). Because of the medium dependency of the total field, this approach represents a nonlinear inverse problem which may be solved via an iterative determination of both the medium properties and the total pressure field. Using the full-wave equation for 3D breast imaging is

(52)

computationally very challenging because modeling the nonlinear inversion is time-consuming, requires definition of many model parameters, and might need many iterations before converging to a reasonable solution.

Several inversion methods have been developed and many of them are summarized in (Lavarello and Hesford, 2013), such as the alternating variables methods, Newton-type methods, conjugate-gradient methods, Kaczmarz-like inversion, eigenfunction methods, and the T-matrix formulations. Wiskin et al. (Andre et al., 2013; Wiskin et al., 2013) have successively shown the feasibility of the 3D implementation of a conjugate gradient based inversion algorithm for the Techniscan Medical Systems (now owned by CVUS LLC, Salt Lake City, UT) scanner. This inversion approach uses an approximate factorization of the Helmholtz wave equation and requires a very fast solution method for the forward problem. Zhang et al. (Zhang et al., 2002, 2004) also investigated different methods to estimate the conjugate gradient directions and different ways to include regularization for improved robustness. Contrast source inversion (CSI) method, developed by van den Berg and Kleinman (van den Berg and Kleinman, 1997), is a very successful variant of the conjugate gradient methods and it will be presented here how it can be applied to breast imaging. Before describing the inversion method, the integral equation formulation of the inverse problem is explained in the next section.

4.2 Inverse Problem

In scattering theory, the total pressure field, recorded by the receivers, is equal to the summation of the incident field ˆpinc(x), generated by the sources and propagating in the homogeneous background medium, and the scattered field ˆpsct(x), hence

ˆ

p(x) = ˆpinc(x) + ˆpsct(x). (4.1) To describe the nonlinear inverse problem, we consider a model which is given by a bounded, simply connected, inhomogeneous object domain D in an unbounded homogeneous background medium. The scattering object of which contrast and location are unknown is placed in the object domain

(53)

4.3 METHODS 43

D_{. The sources and receivers are located in a domain (or on a curve) S} surrounding D (See Figure 4.1).

Assuming the density to be constant and the medium to be lossless, the integral equation formulation of the total pressure wavefield ˆp(x) in D is given in the frequency domain as

ˆ psct(x) = ω2 Z x′∈D ˆ G(x − x′)χ(x′)ˆp(x′)dV. (4.2) Here, the contrast function χ(x′) equals to

χ(x′) = 1 c2_(x′₎−

1

c2₀. (4.3)

During imaging, the unknown contrast function in the domain D is computed using the known total field, which is typically measured by a set of receivers located in the spatial domain S, and the known incident field, which may be computed easily as the sources generating the field are known as well. Solving the inverse problem for a limited number of measurements is a challenging task, as both unknowns, the contrast function and the total field, appear inside the integral. Because of the χ(x′) dependency of the total field, this problem represents a nonlinear inverse problem. Various methods exist to solve the inverse problem approximately or fully, to compute an image of the breast from the measured pressure field. In the next section, some of the solution methods will be discussed in detail.

4.3 Methods

Several imaging methods are applied on different medical applications varying from ray based methods (some of them similar to the ones used for X-ray mammography) (Kak and Slaney, 1988; Greenleaf, 1983) via doppler based (Greenleaf et al., 1987) up to (full-) wave methods (Duric et al., 2007; Jirik et al., 2012; Hesse and Schmitz, 2012; Simonetti et al., 2007b; van Dongen and Wright, 2006, 2007; Hesse et al., 2013). An overview of the six different imaging methods applied in this thesis is given in this section. Of the six methods, two are operating within the ray approximation or high-frequency limit. The first one is the inverse Radon (IR) method, which is based on the Radon transform (Kak and

(54)

Slaney, 1988; van Dongen and Wright, 2006). In literature, this method is also referred to as image reconstruction from projections or filtered back-projection (Ammari, 2008; Devaney, 1982a). The second one, is the refraction corrected time-of-flight tomography . This method allows for bent rays which satisfy the Eikonal equation (Greenleaf, 1983; Dapp, 2014). The Eikonal equation is the infinite frequency limit of the wave equation and is able to account for refraction. The resulting nonlinear inverse problem is solved by an adapted minimization algorithm, which uses total variation for regularization. The main differences with the back-projection is that it uses a regularization term. The remaining four methods are wave based. The first wave based method considered is the synthetic aperture focusing technique (SAFT), which is also known in literature as time-domain migration or delay-and-sum imaging. This method is typically used to compute reflectivity images (Cobbold, 2007). Next, back-propagation (BP) is tested. This single step frequency-domain method is based on the reverse process of the forward propagation, i.e., computing the total field for a known source distribution (van Dongen and Wright, 2006; Ammari, 2008; Devaney, 1982a). This method yields similar results as obtained in the first iteration step of the Born inversion (BI) method, which is tested as well. With this method, the inverse problem is formulated in the frequency domain via an integral equation formulation. This integral equation is first linearized by applying the Born approximation, and successively solved for the unknown contrast using an iterative scheme (van Dongen and Wright, 2006; Abubakar et al., 2004). Finally, a nonlinear inversion method contrast source inversion (CSI) is investigated. This iterative method operates beyond the Born approximation and solves the full nonlinear inverse problem by reconstructing the contrast and contrast sources simultaneously (van Dongen and Wright, 2007; van den Berg and Abubakar, 2001).

A schematic representation of an ideal scanning system is shown in Figure 4.1. Sources and receivers are located inside the domain S at xs and xr respectively. The domain S, referred to as the data domain, is enclosing the region of interest D, referred to as the object domain. Finally, the subscript j = 1, 2, ..., J is used to denote a unique source/receiver combination.

(55)

4.3 METHODS 45

Figure 4.1: Schematic representation of the scanning system. The sources at xs and receivers at xr

are located in the data domain S, which encloses the object domain D.

4.3.1 Inverse Radon (IR)

Inverse Radon transform is mainly used to reconstruct the image of the cross section of an object from projection data. Radon transform represents the projections of the object and is defined by line integrals (integral of some parameter of the object along a straight line).

Based on the assumption that the pressure field travels along straight lines or rays, and that the presence of inhomogeneities in the speed of sound only results in a phase shift of the measured signal (i.e., a change in travel time), the forward problem may be formulated as a Radon transform. Considering the travel time for a wave propagating from a source to a receiver, the corresponding Radon transform in a single plane (x, y) can

Ultrasound Imaging Methods for Breast Cancer Detection

Ultrasound Imaging Methods

for

Ultrasound Imaging Methods

for

Breast Cancer Detection

Contents

Chapter

1

Introduction

1.2

Breast Imaging Systems with Ultrasound

1.3

Aim of the Study

1.4

Thesis Outline

Chapter

2

Theory

2.1

Acoustic Wave Equations

2.2

Integral Equation Formulation

Born approximation

Chapter

3

Forward Scattering Problem

3.1

Full-Wave Method

3.2

Born Approximation

3.3

Incident Field Computation

3.4

Results

3.5

Summary

Chapter

4

Imaging and Inversion

4.1

Imaging with Ultrasound Tomography

4.2

Inverse Problem

4.3

Methods