An iterative reconstruction algorithm for quantitative tissue decomposition using DECT

(1)

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

An iterative reconstruction algorithm for

quantitative tissue decomposition using DECT

Examensarbete utfört i Datorseende vid Tekniska högskolan vid Linköpings universitet

av

Oscar Grandell

LiTH-ISY-EX--12/4617--SE

Linköping 2012

Department of Electrical Engineering Linköpings tekniska högskola Linköpings universitet Linköpings universitet SE-581 83 Linköping, Sweden 581 83 Linköping

(2)

(3)

An iterative reconstruction algorithm for

quantitative tissue decomposition using DECT

Examensarbete utfört i Datorseende

vid Tekniska högskolan i Linköping

av

Oscar Grandell

LiTH-ISY-EX--12/4617--SE

Handledare: Maria Magnusson

Datorseende, ISY och Radiofysik, IMH, Linköpings universitet Alexandr Malusek

Radiofysik, IMH, Linköping universitet Gudrun Alm Carlsson

Radiofysik, IMH, Linköping universitet

Examinator: Maria Magnusson

Datorseende, ISY och Radiofysik, IMH, Linköpings universitet

(4)

(5)

Avdelning, Institution

Division, Department Division of Computer Vision Department of Electrical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2012-09-10 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.cvl.isy.liu.se/ http://www.ep.liu.se ISBN — ISRN LiTH-ISY-EX--12/4617--SE

Serietitel och serienummer

Title of series, numbering

ISSN

—

Titel

Title

En iterativ rekonstruktionsalgoritm för kvantitativ vävnadsklassificering via DECT

An iterative reconstruction algorithm for quantitative tissue decomposition using DECT Författare Author Oscar Grandell Sammanfattning Abstract

The introduction of dual energy CT, DECT, in the field of medical healthcare has made it possible to extract more information of the scanned objects. This in turn has the potential to improve the accuracy in radiation therapy dose plan-ning. One problem that remains before successful material decomposition can be achieved however, is the presence of beam hardening and scatter artifacts that arise in a scan. Methods currently in clinical use for removal of beam hardening often bias the CT numbers. Hence, the possibility for an appropriate tissue de-composition is limited.

Here a method for successful decomposition as well as removal of the beam hard-ening artifact is presented. The method uses effective linear attenuations for the five base materials, water, protein, adipose, cortical bone and marrow, to perform the decomposition on reconstructed simulated data. This is performed inside an iterative loop together with the polychromatic x-ray spectra to remove the beam hardening

Nyckelord

(6)

(7)

Abstract

The introduction of dual energy CT, DECT, in the field of medical healthcare has made it possible to extract more information of the scanned objects. This in turn has the potential to improve the accuracy in radiation therapy dose plan-ning. One problem that remains before successful material decomposition can be achieved however, is the presence of beam hardening and scatter artifacts that arise in a scan. Methods currently in clinical use for removal of beam hardening often bias the CT numbers. Hence, the possibility for an appropriate tissue de-composition is limited.

Here a method for successful decomposition as well as removal of the beam hard-ening artifact is presented. The method uses effective linear attenuations for the five base materials, water, protein, adipose, cortical bone and marrow, to perform the decomposition on reconstructed simulated data. This is performed inside an iterative loop together with the polychromatic x-ray spectra to remove the beam hardening

Sammanfattning

Introduktionen av dual energy CT, DECT, inom medicinsk sjukvård har gjort det möjligt att utvinna mer information om de skannade objekten, vilket har poten-tial att förbättra nogrannheten inom dosplanering vid strålterapi. Ett problem som dock återstår innan man framgångsrikt kan genomföra vävnadsklassificering är förekomsten av strålhärdnings- och spridningsartefakter som uppstår i en scan-ning. Nuvarande metoder som används i kliniska tillämpingar för avlägsnande av strålhärdning ändrar ofta CT-numren. Därmed begränsas möjligheten för en till-fredsställande uppdelning.

Här presenteras en metod för framgångsrik uppdelning såväl som avlägsnande av strålhärdningsartefakter. Metoden använder de effektiva linjära dämpningskoeffi-cienterna för de fem basmaterialen, vatten, protein, fett, kortikalt ben och märg, för uppdelning av den rekonstruerade simulationsdatan. Detta görs innuti en ite-rativ loop tillsammans med det polykromatiska röntgenspektrat för att åtgärda strålhärdningen.

(8)

(9)

Acknowledgments

I would like to start by giving a great thanks to my examinator Maria Magnusson for the opportunity to take part in this exciting research project that this thesis has been a part of. Also a great thanks to Alexandr Malusek and Gudrun Alm Carlsson whom both are part of the research group I have had the opportunity to work with. Thanks for all input and comments during these weeks!

I would also like to thank Karl Stierstorfer from Siemens Health-care for giving me the possibility to work with drasim and providing much of the data that has been used in this thesis.

I would also like to thank all my friends at the university whom made these my best years ever, and a special thanks to my girlfriend Emilia. Finally my parents, Doris and Tommy, and my sister, Ida, deserves the biggest of thanks for the love and support you have given me all my life.

Oscar Grandell

Linköping. 22 August, 2012.

(10)

(11)

Introduction

As a cooperation between the Division of Computer Vision, Department of Elec-trical Engineering, Linköping University and Division of Radiation Physics, De-partment of Medical and Health sciences, Linköping University a novel method for tissue decomposition has been developed [1], [2] and [3]. By using the extra infor-mation given by the dual-energy CT (DECT), compared to single-energy CT, the possibility to improve the quantitative tissue classification increases. This could lead to improved accuracy in radiation therapy dose planning, which is the main focus regarding this work.

Some problems that remains however is the effect of beam hardening and scatter effects in the reconstructed image. Current methods for removal of the beam hard-ening often bias the CT numbers in the reconstructed image. The new method currently being developed in Linköping however, has the potential to remove the artifacts while keeping the CT numbers intact.

The implementation presented in this thesis is based upon a dual energy CT, working at tube voltages of 80 kV and 140 kV. At 140 kV an extra tin filter is used to make the spectrum harder and thus reduce the overlap of spectra. The method is however not dependent on these voltages and can be adjusted to others.

1.1 Purpose and goal

The beam hardening artifact, appearing as shades and streaks in Fig. 1.1, that arises when using polychromatic radiation in CT scanning, can alter the CT val-ues to the extent that they can not be segmented and decomposed correctly. For example using the tissue classification in Section 2.10 and the three material de-composition (3MD) and two material dede-composition plus density (2MD+ρ) in Sections 2.12 and 2.13, respectively, on an image with strong beam hardening might classify bone as soft tissue, then trying to decompose it into soft tissue base materials.

The three material decomposition method has successfully been implemented for images containing only soft tissue [1], [2] and [3]. The main purpose of this thesis is to extend that method to include bone tissue. From that follows an

(16)

(a) Polychromatic CT scan (b) Monochromatic CT scan

Figure 1.1: Example of beam hardening artifacts in a polychromatic CT scan compared to a monochromatic CT scan.

tation of the tissue classification to allow evaluation on simulated data including bone and soft tissue. Later, evaluation on data containing noise will be done as well.

1.2 Problem formulation

As mentioned earlier, the main problem when performing tissue decomposition is the beam hardening artifact, Fig. 1.1. The streaks resulting from it can alter the CT values to the extent that a tissue classification can fail. By extension the decomposition will be into the wrong types of base materials, for example bone tissue might be decomposed into soft tissue base components.

This thesis will extend the method presented in [1], [2] and [3] to include tissue classification and decomposition of bone tissue as well as soft tissue. This in-cludes the introduction of a tissue classifier and separate decomposition for the soft and bone tissue. The method for the soft tissue will be the three material decomposition, and for bone tissue the two materials plus density decomposition.

1.3 Approach

The idea of the iterative algorithm is to estimate via a computer simulation the difference between polychromatic and monochromatic projections. Then add this to the original CT scan and in that way remove the beam hardening artifact, see Fig. 1.2. In truth this addition is done for the projection data, Fig. 1.2 however shows it for reconstructed images which hopefully gives a clearer visual aid to the reader. If the iterative algorithm removes the beam hardening fully, the application of the 3MD and 2MD+ρ on the reconstructed volume should render a correct decomposition.

The idea behind using the 2MD+ρ instead of the 3MD on bone tissue, for example cortical bone, red and yellow marrow, is that a mixture of bone with high mass

(17)

1.4 Related research 5

density and bone marrow with low mass density creates a mixture with a mass density that differs much from the original values of both components. In other words the relationship between the densities of the individual components and the mixtures density is not linear. However, with the calculated density that is the result of the 2MD plus density, correct projections can be calculated.

(a) An image profile showing the linear attenuation coef-ficient (LAC) of a polychro-matic CT reconstruction.

(b) An image profile showing the LAC of a monochromatic CT reconstruction.

(c) Difference between the poly- and monochromatic CT reconstructions, added to the original CT scan.

Figure 1.2: Illustration of the iterative algorithm approach.

1.4 Related research

The research in the field is divided into two areas, the first being methods for beam hardening correction and the other methods for tissue decomposition. Two methods for beam hardening correction are presented by Fuchs [4] and De Man [5]. Fuchs uses a method, called iterative beam hardening correction (IBHC), partly similar to the one used in this thesis. The IBHC is based on correction via filtered back projection. The reconstructed volume is supposed to be composed by a number of basic substances, and that each voxel contains a mixture of two nearby basic substances. However, this method only uses single energy CT, hence each mixture can only be a mixture of two basic substances. The method devel-oped by De Man is called iterative maximum-likelihood polychromatic algorithm for CT (IMPACT). Like Fuchs method it is supposed that each voxel contains a mixture of two nearby basic substances. The attenuation value µ(E) is divided into two components, the photoelectric effect Φ(E) and Compton scattering Θ(E) which gives µ(E) = φΦ(E) + θΘ(E) and only two partial volumes, φ and θ, need to be processed.

Another approach for beam hardening correction is presented by Van Gompel in [6]. Here it is presumed that the spectrum can be represented with a small and predefined number of energy bins, so that no prior knowledge of the spectrum is needed. The object is first segmented, assigning each voxel to one of the mate-rials. Then an estimation of the amplitude of each energy bin is made and used to estimate the attenuation coefficient of each material. The parameter estima-tion is treated as a minimizaestima-tion problem, minimizing the difference between the observed values and the polychromatic projections calculated from the estimated

(18)

parameters.

For tissue decomposition the work in [7] and [8], uses a somewhat similar approach as in this thesis. But it’s not an iterative method and instead of using an effective energy as in this thesis they approach the decomposition with different effective densities for each basic substance.

(19)

Chapter 2

Theory

In this chapter the foundation theory upon which this thesis stands will be de-scribed. In Sections 2.1 to 2.6 the algorithms of computed tomography theory is described and in Section 2.11 the beam hardening artifact is presented. These sec-tions are based on the work in [9], [10], and [11]. Finally in Secsec-tions 2.12 and 2.13 the three material decomposition and two material plus density decomposition will be described, which are based on the work in [1], [2], [3] and [12].

2.1 General projection generation

When generating projections, what is really calculated is line integrals through a 2D or 3D volume. By using the coordinate system in Fig. 2.1, an arbitrary line through the object, represented by a 2D function f (x, y), can be described by

x cos θ + y sin θ = ρ. (2.1)

A line integral can be calculated as

Pθ(ρ) =

Z

(θ,ρ)line

f (x, y)ds. (2.2)

By using a dirac delta function (2.2) can be rewritten as

Pθ(ρ) = ∞ Z −∞ ∞ Z −∞

f (x, y)δ(x cos θ + y sin θ − ρ)dxdy. (2.3)

This is also known as the Radon transform. By combining a set of line integrals taken at different ρ but with constant θ, a projection is formed. The projections can be calculated using different settings, the simplest is collection of parallel rays known as parallel projections, Fig. 2.2a. A more common approach nowadays, especially in CT, is to use fan-beam projections, Fig. 2.2b. Here a single source sending out x-rays in a fan pattern is used, thus called fan-beam projections.

(20)

Figure 2.1: An object f (x, y) and its projection Pθ(ρ) at angle θ.

2.2 Projection generation in CT

The projections generated by the CT scan is actually measurements of the atten-uation of the material. In short a set of beams are sent through the object and an attenuation profile is measured on the other side. From the result obtained from this single projection there is no way of determining any major attributes of the scanned object, but by repeating from different angles and then projecting this measured attenuation profile over an empty image grid the object can be recon-structed. This technique is more commonly known as backprojection, which will be described in detail in Section 2.3.

(a) Parallel. (b) Fan-beam.

(21)

2.2 Projection generation in CT 9

2.2.1 Monochromatic projections

The beam attenuation in the material is mainly due to three different effects: the photoelectric effect, the Rayleigh scattering and the Compton scattering. The photoelectric absorption consists of a photon transferring all its energy to a tightly bound inner electron in an atom, while the Compton scattering is interaction between loosely bound or free electrons, making the photon change direction and loose some of its energy to the electron. The Rayleigh scattering is a function of the electric polarizability of the particles, and like the Compton scattering causes the photons to change direction.

Figure 2.3: An x-ray tube sending out monochromatic photons through a thin slice of an unknown object.

Now consider the setup in Fig. 2.3, here N monochromatic photons travel through an unknown substance of width ∆x and hit the detector. Due to the attenuation,

N − ∆N will hit the detector. This can be described as in (2.4), where the attenuation coefficient µ is the probability per unit distance that the photon will be removed from the pencil beam, due to the combined effect of photoelectric absorption, Rayleigh and Compton scattering.

∆N

N

1

∆x = −µ (2.4)

If ∆x goes to zero in (2.4) a differential equation is obtained, 1

NdN = −µdx, (2.5)

which can be solved as

N Z N0 1 NdN = −µ x Z 0 dx ⇒ lnN N0 = −µx ⇒ N (x) = N0e −µx_, _(2.6)

(22)

where N (x) is the number of photons at distance x. Here the assumption that

µ is constant has been made, however is this not true most of the time. Instead

consider Fig. 2.2, where a single x-ray beam travels in the s-direction through a function, f (x, y), describing the attenuation of the material. It can then be described as Z ray µ(x, y)ds = lnNs Nd , (2.7)

where Nsis the number of photons emitted from the source and Nd the number

of photons detected.

2.2.2 Polychromatic projections

In the real world monochromatic x-ray sources are hard to find and x-ray spectra are used instead, see Fig. 2.4. Since the attenuation coefficient is dependent on the energy, (2.7) has to be rewritten to

Nd = Z Sin(E)e− R µ(x,y,E)ds dE, (2.8)

where Sin(E) is the incident photon number density for energy level E.

Figure 2.4: Example of a typical X-ray spectrum.

2.3 Backprojection

As mentioned in Section 2.2, the gathered data are smeared back over the image grid to reconstruct the object. This can mathematically be expressed as proposed in [10],

fθk(x, y) = Pθk(ρ) = Pθk(x cos θk+ y sin θk), (2.9)

where fθk(x, y) is the result of backprojection at a certain angle, θk. Since this is

true for all angles θ (2.9) can be rewritten as

(23)

2.4 The Fourier-Slice Theorem 11

The final image is then formed by integrating over all projection angles,

f (x, y) = 2π

Z

0

Pθ(x cos θ + y sin θ)dθ. (2.11)

This will however result in a blurred image which is only an approximation of the true image.

2.4 The Fourier-Slice Theorem

The Fourier-slice theorem states the relationship between the 1D Fourier trans-form of a projection and the 2D Fourier transtrans-form of the region from which the projection was taken. The derivation here is similar to the one in [10]. This rela-tion is the basis for reconstrucrela-tion without blurring. The 1D Fourier transform of a projection can be written as

Sθ(ω) =

∞

Z

−∞

Pθ(ρ)e−j2πωρdρ. (2.12)

By inserting (2.3) in (2.12) the 1D Fourier transform can be written as

Sθ(ω) = ∞ Z −∞ ∞ Z −∞ ∞ Z −∞

f (x, y)δ(x cos θ + y sin θ − ρ)e−j2πωρdxdydρ (2.13)

= ∞ Z −∞ ∞ Z −∞ f (x, y)[ ∞ Z −∞

δ(x cos θ + y sin θ − ρ)e−j2πωρdρ]dxdy (2.14)

= ∞ Z −∞ ∞ Z −∞

f (x, y)e−j2πω(x cos θ+y sin θ)dxdy, (2.15)

and by using u = ω cos θ and v = ω sin θ the expression becomes

Sθ(ω) = [ ∞ Z −∞ ∞ Z −∞

f (x, y)e−j2π(ux+vy)dxdy]u=ω cos θ;v=ω sin θ, (2.16)

which is the 2D Fourier transform of f (x, y). The Fourier-slice theorem can be summed up as in (2.17), which states that Fourier transform of a projection at a certain angle, θ, is located along a line at the same angle of the 2D Fourier transform of the region from which it was obtained, see also Fig. 2.5.

(24)

Figure 2.5: The Fourier slice theorem.

2.5 Filtered backprojection

Filtered backprojection is the most common way of reconstruction within com-puted tomography, but as stated in Section 2.3 normal backprojection only recon-structs the object approximately and a blurred image is obtained. One way to deal with this is to filter the projection data and use the properties of the Fourier slice theorem described in Section 2.4.

2.5.1 Parallel-beam filtered backprojection

Here follows a derivation of the filtered backprojection for a parallel-beam gen-erated projection based on [11]. The object function f (x, y) can be obtained by applying the inverse Fourier transform on the Fourier transform of the object,

f (x, y) = ∞ Z −∞ ∞ Z −∞

F (u, v)ej2π(ux+vy)dudv. (2.18)

A substitution from the rectangular coordinate system, (u, v), in the Fourier do-main to a polar coordinate system, (ω, θ) via,

u = ω cos θ, v = ω sin θ, dudv = ωdωdθ, (2.19)

gives the Fourier transform of a polar function as

f (x, y) = 2π Z 0 ∞ Z 0

(25)

2.5 Filtered backprojection 13

The integral in (2.20) can be split into two by considering 0 ≤ θ < π and π ≤ θ < 2π which gives f (x, y) = π Z 0 ∞ Z 0

F (ω, θ)ej2πω(x cos θ+y sin θ)ωdωdθ

+ π Z 0 ∞ Z 0

F (ω, θ + π)ej2πω[x cos(θ+π)+y sin(θ+π)]ωdωdθ. (2.21)

When using the Fourier property in (2.22) and combining it with (2.21), the ex-pression in (2.23) is obtained. F (ω, θ + π) = F (−ω, θ) (2.22) f (x, y) = π Z 0   ∞ Z −∞

F (ω, θ)|ω|ej2πω(x cos θ+y sin θ)dω



dθ (2.23)

The equation above can be simplified by setting ρ = x cos θ + y sin θ. The 2D Fourier transform F (ω, θ) can be exchanged to the Fourier transform of the pro-jection at angle θ, Sθ(ω), to get

f (x, y) = π Z 0   ∞ Z −∞ Sθ(ω)|ω|ej2πωρdω  dθ, (2.24)

which also can be expressed as

f (x, y) = π Z 0 Qθ(ρ)dθ = π Z 0 Qθ(x cos θ + y sin θ)dθ, (2.25) where Qθ(ρ) = ∞ Z −∞ Sθ(ω)|ω|ej2πωρdω. (2.26)

The last equation represents a filtering operation, where |ω| is a 1D filter in the Fourier domain called a ramp filter, and Qθ(ρ) is the filtered projection. Equation

(2.25) states that the object can be reconstructed by summation of all these filtered projections.

2.5.2 Fan-beam filtered backprojection

Nowadays not many parallel-beam CT scanners exist, this due to the fact the scan times become very long as the x-ray source has to linearly scan the object at all angles. Instead modern CT scanners use an x-ray point source that emits a fan-shaped beam, Fig. 2.2b, or even more modern using multiple fan-beams

(26)

at different pitch angles. These are referred to as multi-slice spiral CT or helical cone-beam scanners. As the beams no longer travel parallel to each other, the filtered backprojection has to be changed as well.

There are two types of fan-beam projections depending upon if the projections are sampled at equiangular or equispaced intervals. Here follows a derivation, similar as in [11], of how the filtered backprojection is altered in the case of equiangular rays.

Figure 2.6: Equiangular geometry.

First consider the setup in Fig. 2.6 and let Rβ(γ) denote the projection obtained

from a fan-beam. β is then the angle the source makes with a reference axis y, and the angle γ is the location of a ray within the fan in relation to the beam that hits the center of the detector. If the projection data would have been generated along a set of parallel rays and belonged to a parallel projection Pθ(ρ) it would

have to satisfy the conditions

θ = β + γ, ρ = D sin γ, (2.27)

where D is the distance from the source to the origin. The object f (x, y) can then be reconstructed as f (x, y) = π Z 0 ρm Z −ρm Pθ(ρ)h(x cos θ + y sin θ − ρ)dρdθ, (2.28)

where ρm is the value of ρ for which Pθ(ρ) = 0 with |ρ| > ρm in all projections.

Equation (2.28) only requires parallel projections over an 180◦ interval, for recon-struction from fan-beam data, however, data from a 360◦ is necessary and the

(27)

2.5 Filtered backprojection 15

equation has to be rewritten as

f (x, y) = 1 2 2π Z 0 ρm Z −ρm Pθ(ρ)h(x cos θ + y sin θ − ρ)dρdθ. (2.29)

Derivation of the algorithm becomes easier when a point (x, y) is expressed in polar coordinates (r, φ) where

x = r cos φ, y = r sin φ. (2.30)

Then (2.29) can be written as

f (r, φ) = 1 2 2π Z 0 ρm Z −ρm Pθ(ρ)h(r cos(θ − φ) − ρ)dρdθ. (2.31)

Using (2.27), this integral can be expressed in terms of β and γ as

f (r, φ) = 1 2 2π−γ Z −γ sin−1ρm/D Z − sin−1_ρ m/D

Pβ+γ(D sin γ)h(r cos(β +γ −φ)−D sin γ)D cos γdγdβ.

(2.32) Here a couple of observations can be made. Firstly the interval of −γ to 2π − γ for β cover the entire 360◦ range that is needed, and since all functions of β are 2π-periodic the limits can be changed to 0 to 2π. The term sin−1(ρm/D) is equal

to the outermost ray that hits the detector, which is equal to γmso the limits can

be replaced with −γm to γm. Last, the expression Pβ+γ(D sin γ) corresponds to

a ray integral in the parallel projection data Pθ(ρ), and its correspondence in the

fan-beam geometry is Rβ(γ). This results in that (2.32) can be rewritten as

f (r, φ) = 1 2 2π Z 0 γm Z −γm

Rβ(γ)h(r cos(β + γ − φ) − D sin γ)D cos γdγdβ. (2.33)

This equation can be further simplified, for proof see [11], by using the following relations from Fig. 2.6

r cos(β + γ − φ) − D sin γ = L sin(γ0− γ), (2.34)

where γ0= tan−1 r cos(β − φ) D + r sin(β − φ), (2.35) and h(L sin γ) = _γ L sin γ 2 h(γ). (2.36)

(28)

This allows (2.33) to be written as f (r, φ) = 2π Z 0 1 L2 γm Z −γm Rβ(γ)g(γ0− γ)D cos γdγdβ, (2.37) where g(γ) = 1 2 _γ sin γ 2 h(γ). (2.38)

Equation (2.37) can be interpreted as a weighted filtered backprojection, since it can be written as f (r, φ) = 2π Z 0 1 L2Qβ(γ 0_)dβ, _(2.39) where Qβ(γ) = R0β(γ) ∗ g(γ), (2.40) and R0_β(γ) = Rβ(γ)D cos γ. (2.41)

2.6 Fan-beam to parallel-beam rebinning

Another approach is to transform the fan-beam data to parallel data, this process is also known as rebinning. Once the data has been rebinned we can use parallel filtered backprojection to reconstruct. This approach has a few advantages com-pared to the weighted fan-beam filtered backprojection, with the greatest being that it is not as computationally heavy as the fan-beam filtered backprojection. Also, the angular scanning interval can be reduced. Another advantage in this thesis is that the projection calculations in the iterative algorithm, see Section 3.2, becomes easier to calculate. Here follows a derivation of how the rebinning can be done, similar as in [11]. As was seen in Section 2.5.2 the relationship between fan-beam projections and parallel projections are given by (2.27), namely

θ = β + γ, ρ = D sin γ, (2.42) or equivalently β = θ − γ, γ = sin−1ρ D . (2.43)

So by using interpolation and the relationship in (2.43) it is possible to rebin the fan-beam data to parallel data and use parallel filtered backprojection. Another observation to be made is that parallel projections that are 180◦ apart are mirror images of each other,

Pθ(ρ) = Pθ+π(−ρ). (2.44)

This implies that to reconstruct the object, a fully 360◦interval with the fan-beam projections is no longer necessary, only an interval that lets us rebin to a 180◦

(29)

2.6 Fan-beam to parallel-beam rebinning 17

interval. The relation in (2.44) also implies that an object can be reconstructed as long as we have

θ0≤ θ < θ0+ π. (2.45)

By extending this to the beam case, two ray integrals represented by the fan-beam angles (β1, γ1) and (β2, γ2) are identical if

β1− γ1= β2− γ2+ π, (2.46)

where

γ1= −γ2, (2.47)

for full proof see [11]. This means that a scanning interval for the fan-beam of only 180◦ is not enough, the scan has to be over 180◦+ 2γminterval to be able to

rebin to a 180◦ interval without losing data, Fig. 2.8 compared to Fig. 2.7. Here

γm is the angle from the central detector element to the outermost element, see

Fig. 2.6.

(a) Fan-beam data. (b) Parallel data.

Figure 2.7: Rebinning of fan-beam data, shaded, to parallel with too few angles. The area A will be unnecessary when rebinned, and area B will be missing.

(30)

(a) Fan-beam data. (b) Parallel data.

Figure 2.8: Rebinning of fan-beam data, shaded, to parallel data. By scanning the extra 2γm a full 180◦ interval can be obtained.

2.7 Forward projection

In the iterative algorithm, see Section 3.2, both monochromatic and polychromatic projections are calculated. This can be done in different ways, a common way is to use Joseph’s Method [13], which is presented here.

The method used by Joseph, assumes that the object function, f (x, y), is slowly varying so that linear interpolation can be used between the pixels and that the image consists of N × N pixels in a Cartesian grid. Then consider an exactly straight ray K specified as

y(x) = − cot(θ) + y0, (2.48)

or

x(y) = −y tan(θ) + x0, (2.49)

where θ is the angle to the y-axis and y0/x0 is the cross-point with the

y-axis/x-axis. The desired line integral P (K) is dependent on whether the rays direction lies closer to the x or y axis and can be written as

P (K) =

( ₁

| sin θ|R f (x, y(x))dx for | sin θ| ≥ 1 √ 2 1

| cos θ|R f (x(y), y)dy for | cos θ| ≥ 1 √ 2

. (2.50)

These can then be approximated by a Riemann sum and the x-direction integral becomes P (K) = 1 | sin θ| "N −1 X n=2 Fn,n0+ λ_n(F_n,n0₊₁− F_n,n0) + T₁+ T_N # , (2.51)

(31)

2.8 Monochromatic and polychromatic projections 19

where T1 and TN represents the first and last pixels along each ray, Fn,n0 =

f (xn, yn0) and λ_n is defined as

λn= y(xn) − n0, n0= integer part of y(xn), (2.52)

by rewriting the inner part of (2.51)

Fn,n0+ λ_n(F_n,n0₊₁− F_n,n0) = (1 − λ_n)F_n,n0+ λ_nF_n,n0+1, (2.53)

the linear interpolation becomes even more clear.

2.8 Monochromatic and polychromatic projections

The monoenergetic projection PmE at effective energy E1for base material m with

linear attenuation coefficient µmE can be calculated with Joseph’s method which

gives,

PmE1 = µmElm. (2.54)

The polyenergetic projection is also calculated with Joseph’s method, and by using,

PU1 = ln  Sin S0 = − ln S0 Sin , (2.55)

where Sin is the incident-photon intensity and S0 is the exiting-photon intensity

given by Sin= Emax X 0 EN (E), (2.56) S0= Emax X 0 EN (E)e− PM m=1µm(E)lm = Emax X 0 EN (E)e− PM m=1σm(E)ρmlm , (2.57) where E is the photon energy, N (E) is the number of photons, µm(E) is the linear

attenuation coefficient and ρm the density of material m and lm is the length of

intersection with material m calculated by using (2.54). Since µm(E) = ρmσm(E),

where σm(E) is the mass attenuation coefficient, the calculated density given by

the decomposition methods 3MD or 2MD+ρ, Section 2.12 and 2.13, should be used instead of the predefined density of the pure material. This is especially necessary for bone decomposition with 2MD+ρ since the mixtures density is not linearly dependent on the components densities. For the 3MD method ρ is given by (2.71) and for the 2MD+ρ, ρ is given by

ρ = _w 1 1 ρ1 + w2 ρ2 . (2.58)

(32)

2.9 Partial volume

One effect that often occurs when forward projection is preformed is the partial volume effect. It is an effect that happens when summation a long a ray that passes near edges between two areas with large differences for their linear attenuation coefficients, for example water and air. What happens is that the interpolation described in Section 2.7 gives a value between the two values which in turn gives a bad result in the 3MD or 2MD+ρ decomposition.

Since this problem only occurs between water and air for this study, the problem is solved so that the decomposed elements along the intersection between water and air are recalculated as the mean of pixels, with a value above a threshold, inside an area of 5 × 5 pixels around the pixel in question.

2.10 Tissue classification

In CT images the tissues are represented by the linear attenuation coefficient, in DECT two such are obtained. By using these, a linear attenuation diagram can be created, see Fig. 2.9a, where regions for lung, soft and bone tissue are indicated.

(a) Linear attenuation diagram. (b) Thresholding of linear attenuation dia-gram.

Figure 2.9: Linear attenuation diagram and simple thresholding.

Different approaches can then be used to classify and segment the tissue. Here two different methods are presented, one being a simple threshold and the other being a more advanced method called mixture of gaussians.

2.10.1 Simple threshold

A common and simple approach to classify the different tissue is to use a threshold as in Fig. 2.9b. Here the current pixel value (linear attenuation coefficients for 80 keV) is simply compared to a predefined set of classifiers and classified as the correct kind of tissue.

(33)

2.11 Beam hardening 21

2.10.2 Mixture of gaussians

Another method is to use the known correct values of the linear attenuation coef-ficients, pj = (µ80, µ140)T, and assign each sample to the closest, in other words

using a nearest neighbor approach. As seen in Fig. 2.10a the classification will with high probability fail. Instead the method can be further extended in a way

(a) Nearest neighbour method. (b) Mixture of gaussians method.

Figure 2.10: Advanced methods for tissue segmentation.

that each correct value is also assigned a Gaussian covariance matrix, Cj. By

using an iterative semi supervised training method the mean of all pixels assigned to each class and the convergence matrix adapts to the spreading of each class. This allows for using other parameters than the Euclidean distance. In fact we classify each pixel according to

arg max j f (xi; pj, Cj) = 1 2π|Cj| 1 2 e−12(xi−pj)TC−1j (xi−pj)_, _(2.59)

j here represent the class. This gives a result as Fig. 2.10b.

2.11 Beam hardening

The beam hardening artifact takes the shape of a cupping of the image or streaks in it. This is due to the interaction between the material and the x-rays, which oc-cur in three different ways, coherent scattering or Raylight scattering, incoherent scattering or Compton scattering and the photoelectric effect. The linear atten-uation coefficient, µ[m−1], measured by the CT scanner can be described as a summation of these, i.e.

µ(E) = ρ(µmCo(E) + µmIn(E) + µmP h(E)), (2.60)

where ρ[kg/m3_{] is the density of the material and µ}

m[m2/kg] is the mass

atten-uation coefficient of the interaction process. A typical appearance of the linear attenuation coefficient with respect to the energy level is shown in Fig. 2.11.

(34)

Figure 2.11: Example of an attenuation curve.

(a) Entering spectrum. (b) Exiting spectrum.

Figure 2.12: Comparison between entering and exiting spectrum.

When a polychromatic beam in the form of an x-ray spectra, Fig. 2.12a, passes the object, the photons with lower energy has a higher probability of being absorbed, see Fig. 2.11 and [9]. This leads to the beam having a more narrow or harder spectra when exiting the object, Fig. 2.12b.

2.12 Three material decomposition

The three material decomposition method, 3MD, presented here is based on [1], [3] and [12].

Assume that a material is built up of three separate components, with mass at-tenuation coefficients µ1

ρ1,

µ2

ρ2 and

µ3

ρ3 respectively. µi is the linear attenuation

coefficient and ρi is the density. The law of mixtures can be applied on the mass

attenuation coefficients if we assume that the mass is conserved, which is true in most cases. This is stated in (2.61) and (2.62).

µ(E) ρ = w1 µ1(E) ρ1 + w2 µ2(E) ρ2 + w3 µ3(E) ρ3 (2.61) w1+ w2+ w3= 1 (2.62)

(35)

2.12 Three material decomposition 23

By measuring the mass attenuation at two different energy levels, E1 and E2, a

set of equations is acquired,

µ(E1) ρ = w1 µ1(E1) ρ1 + w2 µ2(E1) ρ2 + w3 µ3(E1) ρ3 , (2.63) µ(E2) ρ = w1 µ1(E2) ρ1 + w2 µ2(E2) ρ2 + w3 µ3(E2) ρ3 . (2.64)

By using (2.62), the equations can be rewritten as

µ(E1) ρ = w1 µ1(E1) ρ1 + w2 µ2(E1) ρ2 + (1 − w1− w2) µ3(E1) ρ3 , (2.65) µ(E2) ρ = w1 µ1(E2) ρ1 + w2 µ2(E2) ρ2 + (1 − w1− w2) µ3(E2) ρ3 , (2.66) which can also be written as

µ(E1) ρ = µ3(E1) ρ3 + w1  µ1(E1) ρ1 −µ3(E1) ρ3 + w2  µ2(E1) ρ2 −µ3(E1) ρ3 , (2.67) µ(E2) ρ = µ3(E2) ρ3 + w1  µ1(E2) ρ1 − µ3(E2) ρ3 + w2  µ2(E2) ρ2 − µ3(E2) ρ3 , (2.68)

or on a more compact form as

~ M = ~M3+ w1h ~M1− ~M3 i + w2h ~M2− ~M3 i , (2.69) where ~ M = "_µ(E 1) ρ µ(E1) ρ # and M~i= "_µ i(E1) ρi µi(E1) ρi # . (2.70)

This is also illustrated in Fig. 2.13.

Considering the density, ρ, in (2.63) it is an unknown parameter. If the assumption that the volume of the mixture is equal to the sum of the volumes of each individual component is made, the density can be written as

ρ =m V = m m1 ρ1 + m2 ρ2 + m3 ρ3 = _w 1 1 ρ1 + w2 ρ2 + w3 ρ3 . (2.71)

Combining (2.67) and (2.68) with (2.71) gives

µ(E1)−µ3(E1) ρ3 µ(E2)−µ3(E2) ρ3 ! + Mw1 w2 =0 0 , (2.72) where M = "_µ(E₁_)−µ₁_(E₁₎ ρ1 − µ(E1)−µ3(E1) ρ3 µ(E1)−µ2(E1) ρ2 − µ(E1)−µ3(E1) ρ3 µ(E2)−µ1(E2) ρ1 − µ(E2)−µ3(E2) ρ3 µ(E2)−µ2(E2) ρ2 − µ(E2)−µ3(E2) ρ3 # . (2.73)

Solving (2.72) gives the weight fractions w1 and w2 and then w3 can then be

(36)

Figure 2.13: Vector diagram for three material decomposition (3MD).

2.13 Two material plus density decomposition

The two material plus density decomposition, 2MD+ρ, assumes similar to the three material decomposition, that a material is built up by individual compo-nents with mass attenuation coefficients µ1(E)

ρ1 and

µ2(E)

ρ2 , where µi is the linear

attenuation coefficient and ρi is the density of each component. A mixture of the

two materials can then be described according to (2.74). Also assume that the mass of the individual components are equal to the total mass, which in terms of mass fractions, wigives (2.75).

µ(E) ρ = w1 µ1(E) ρ1 + w2 µ2(E) ρ2 (2.74) w1+ w2= 1 (2.75)

By measuring at two different energy levels, E1 and E2, the following equation

system is obtained: µ(E1) ρ = w1 µ1(E1) ρ1 + w2 µ2(E1) ρ2 , (2.76) µ(E2) ρ = w1 µ1(E2) ρ1 + w2 µ2(E2) ρ2 , (2.77)

with (2.75) the equations can be rewritten as:

µ(E1) ρ = w1 µ1(E1) ρ1 + (1 − w1) µ2(E1) ρ2 , (2.78) µ(E2) ρ = w1 µ1(E2) ρ1 + (1 − w1) µ2(E2) ρ2 , (2.79)

(37)

2.13 Two material plus density decomposition 25 or equally µ(E1) ρ = µ2(E1) ρ2 + w1  µ1(E1) ρ1 − µ2(E1) ρ2 , (2.80) µ(E2) ρ = µ2(E2) ρ2 + w1  µ1(E2) ρ1 −µ2(E2) ρ2 , (2.81)

which can be rewritten with the same notation as in (2.70) to

~ M =

~ L

ρ = ~M2+ w1( ~M1− ~M2), (2.82)

where ~L are the linear attenuation coefficients measured by the CT-scanner at

energy levels E1 and E2. Equation (2.82) is illustrated in Fig. 2.14. Since only

two materials make up the mixture and the assumption in (2.75) has been made, the density ρ can be left as a free parameter and the equations can be rewritten as µ2(E1) ρ2 µ2(E2) ρ2 ! + Mw11 ρ =0 0 , (2.83) where M = µ1(E1) ρ1 − µ2(E1) ρ2 −µ(E1) µ1(E2) ρ1 − µ2(E2) ρ2 −µ(E2) ! . (2.84)

which if solved gives w1 and _ρ1 and w2 can be obtained from (2.75).

Figure 2.14: Vector diagram for two material plus density decomposition (2MD+ρ).

(38)

(39)

Chapter 3

Methods

In this chapter a simple beam hardening correction algorithm and the iterative algorithm will be presented, as well as the modifications of them that had to be done. Last, the three phantoms that were used in the evaluation process are presented.

3.1 Simple water beam hardening correction

The simple water beam hardening correction, SWBHC, see for example [11], is based upon the fact that for a monochromatic beam the intensity of the exiting pencil beam can be calculated as (3.1), see Section 2.2.1,

Sexit= Sine−

R

µ(E0,l)dl_, _(3.1)

which in the case of a homogenous phantom (e.g. a water phantom) can be written as

Sexit= Sine−µ(E0)l, (3.2)

while in the case of a polychromatic spectrum the relationship between the incident spectrum and exiting spectrum can be described as in (2.8),

Sexit= Z Sin(E)e− R µ(E,l)dl dE. (3.3)

However, for both cases the projection data for image reconstruction is calculated as Pmono= − ln  Sexit Sin Ppoly= − ln  Sexit Sin , (3.4) where Sin= Emax Z 0 ESin(E)dE. (3.5) 27

(40)

For the monochromatic case this implies a linear relationship between the mea-sured data and the linear attenuation coefficient, by combining (3.1) and (3.4), (3.6) is obtained, Pmono= ln _S in Sexit = µ(E0)l. (3.6)

For the polychromatic case however, this is not true, and a nonlinear relationship is obtained due to the beam hardening. This can be seen in Fig. 3.1, where the obtained projection data ln Sin

Sexit

has been calculated at different thicknesses, l.

Figure 3.1: Projection data, ln Sin

Sexit

, for the monochromatic and polychromatic cases plotted against the thickness, l.

However, with the knowledge of this relationship a simple correction can be made. The approach is to via the polychromatic projection data find the width of the object, l0, and with the use of it find the correct monochromatic projection value.

First an effective linear attenuation coefficient is calculated as the energy-fluence weighted linear attenuation coefficient, for example according to [2], as

µmE = REmax 0 EN (E)µ(E)dE REmax 0 EN (E)dE , (3.7)

where E is the energy, N (E) is the relative intensity of photons and µ(E) is the linear attenuation coefficient.

A normal way to present a CT scan is by Hounsfield numbers, these are calculated according to (3.8). In (3.8) µw is the effective linear attenuation coefficient for

water and corresponds to µmE in (3.7).

HU = µ − µw µw

(41)

3.1 Simple water beam hardening correction 29

µmE is then used to compute the monochromatic projection curve µmEl, see Fig.

3.2. Now the trick is to find the thickness, l, via the polychromatic curve. Once l is found the monochromatic projection value is calculated as

Pmono= µmEl, (3.9)

and a correction factor, ccorr, can be calculated as

ccorr= Pmono Ppoly = µmEl ln Sin Sexit  . (3.10)

Figure 3.2: Graphical interpretation of the simple water beam hardening correction

The correction factor can then be used to correct the polychromatic projection by a simple multiplication, Ppolyccorr= ln _S in Sexit _µ mEl ln Sin Sexit  = µmEl = Pmono. (3.11)

This method works perfectly for a homogeneus substance and fairly well for objects containing similar tissues, for example other soft tissues. However, for bone tissue, which has different attenuation and density, the beam hardening will alter the numbers to much for it to correct them in a similar way as for the chosen material (e.g. water).

(42)

3.2 The iterative algorithm

To remove the beam hardening artifact fully, and to decompose the material an iterative algorithm based on [2] is used. The iterative algorithm is shown in Fig. 3.3 and described in detailed here:

Figure 3.3: The iterative reconstruction algorithm

• Two measured parallel projections, PM,U 1and PM,U 2, for two different

spec-tra corresponding to two different tube voltages, U1 and U2, are input. All

other data are initially set to 0.

• PM,U 1and PM,U 2are reconstructed via filtered backprojection into the

vol-umes µ1 and µ2. The volumes now contain linear attenuation coefficients

approximately similar to the effective energies E1and E2of the spectra for

tube voltages U1 and U2. This corresponds to 1 in Fig. 3.3. The effective

energies E1and E2were found by calculating µmE, (3.7), for water and then

using the attenuation curve, Fig. 2.11, for water.

• The reconstructed volumes, µ1 and µ2, are classified via one of the tissue

classifications in Section 2.10 as either soft tissue or bone tissue. This cor-responds to 2 in Fig. 3.3.

• The soft tissue are then transmitted to the three material decomposition, Section 2.12, and the bone tissue are transmitted to the two material plus density decomposition, 2.13 with a few extensions, see Section 3.2.1. The result is the decomposed volume µC.

• Then monoenergetic projections PE1 and PE2 at the energies E1 and E2,

and polyenergetic projections PU1 and PU2for spectra U1and U2are forward

projected using Joseph’s method, see Section 2.7. This corresponds to 3 and 4 in Fig. 3.3.

(43)

3.2 The iterative algorithm 31

• The monoenergetic projections are added to, and the polyenergetic projec-tions are subtracted from the measured parallel projecprojec-tions, corresponding to 5 in Fig. 3.3 PM0_,U 1 = PM,U1+ PE1− PU1, (3.12) PM0_,U 2 = PM,U2+ PE2− PU2. (3.13) • PM0_,U

1 and PM0,U2 are then submitted to the next iteration. After a certain

number of iterations the beam hardening artifacts are removed and accu-rate mass fractions of the base materials and density for the bone tissue is obtained.

3.2.1 Extensions to the 2MD plus density

To increase the robustness of the 2MD plus density decomposition, see Section 2.13, an extra limitation was introduced

0 ≤ w1≤ 1, (3.14)

0 ≤ w2≤ 1. (3.15)

This limits the weight fraction of each tissue to be between 0 and 1, which seems fully correct since no tissue can be more than 100 % pure and definitely not have a negative weight fraction. With the introduction of these the possibility of decomposition compensation is removed as well, see Example 3.1.

Example 3.1: Decomposition compensation

Suppose that a pixel contains 25 % cortical bone of its total mass, ρ = 1.92 ,and 75% marrow, ρ = 1.00. The marrow used here is a mixture of 50 % red and 50% yellow marrow.

Energy µm(E1) for cortical bone [m2/g] µm(E2) for marrow [m2/g]

49.9 42.83 21.83

88.5 20.30 17.47

The values of the pixel can then be calculated via,

µ(Ei) ρ = w1 µ1(Ei) ρ1 + (1 − w1) µ2(Ei) ρ2 , where µj(Ei)

ρj = µm,j(Ei). This will result in the pixel value of

µ(E1)

ρ = 27.08 for

E1= 49.9 ekV and µ(E2)

ρ = 18.18 for E2= 88.5 ekV. Now solving the equations,

µ(E1) ρ = w1 µ1(E1) ρ1 + (1 − w1) µ2(E1) ρ2 , µ(E2) ρ = w1 µ1(E2) ρ1 + (1 − w1) µ2(E2) ρ2 ,

(44)

or rewritten as µ2(E1) ρ2 µ2(E2) ρ2 ! + µ1(E1) ρ1 − µ2(E1) ρ2 −µ(E1) µ1(E2) ρ1 − µ2(E2) ρ2 −µ(E2) ! w1 1 ρ =0 0 ,

gives w1= 0.5 and since w2 = 1 − w1 = 0.5. However if the value varies a little,

for example due to tissue misclassification or noise, the optimal solution might be

w1= 21.84 which gives w2= −20.84 without the extra conditions. The two extra

limitations are then needed for finding the optimal solution within the correct range.

3.3 Simple water beam hardening correction and

iterative algorithm

As a final method the two previous methods will be combined. That is, in the 0:th iteration the simple water beam hardening correction will be applied before reconstruction and tissue classification. The hope is that the tissue will earlier be correctly classified and decomposed, hence saving iterations in the loop and minimizing the calculations and time consumption.

3.4 Experimental setup

In this section the setup for the simulations is described. First the geometry is presented and then the calculations of the linear attenuation coefficients and effective energies that were used during the simulations are given.

(45)

3.4 Experimental setup 33

3.4.1 Simulation geometry

The geometry of the simulations for studies 1-3, preformed using take, is presented in Fig. 3.4, where α is the fan-beam angle of 5.8◦, D is the source isocenter distance equal to 1 meter. As can be seen a flat detector array of length L = 0.2018m was used, the detector elements were 0.7883 mm, which gives a total of 256 elements. A total of 280 projections were taken at the angular interval, 0◦ ≤ θ ≤ 207.2◦_,

which gives an angular increment of 0.74◦.

Two different spectra were used, one for 80 kV and one for 140 kV with extra tin filtration, these were both provided by Siemens.

3.4.2 Mean and standard deviation calculations

The mean of the region is calculated as

µ = 1

N − 1

X

i∈Rj

wi, (3.16)

and the standard deviation is calculated over the same area as

σ = s 1 N − 1 X i∈Rj (wi− µ)2. (3.17)

N is the number of pixels in each region, Rj, and wi is the weight value for a

certain pixel. It should also be noted that consideration about the partial volume effect has been taken into account in the form that the pixels between two regions are not part of any region, Rj.

3.4.3 Linear attenuation coefficients and effective energies

The linear attenuation coefficients for each material were calculated via (2.60),

µ(E) = ρ(µmCo(E) + µmIn(E) + µmP h(E)),

where the mass attenuation coefficients, µmCo(E), µmIn(E), and µmP h(E), where

taken from [14].

The effective linear attenuation coefficients, µmE, is calculated via the effective

energy of each spectrum. First the effective linear attenuation coefficient for water was calculated as,

µmE = REmax 0 EN (E)µ(E)dE REmax 0 EN (E)dE ,

in accordance with (3.7). The corresponding energy level, Eef f, is then received

from the linear attenuation coefficient curve for water. Using Eef f in the

µ(E)-curve for the other material, the effective attenuation for protein, adipose, cortical bone and marrow is obtained, see Table 3.1. These values are also presented in the linear attenuation coefficient (LAC) diagrams in Figure 3.5-3.6. Note that adipose and water are rather close, which will lead to problems as will be seen in Chapter 4.

(46)

Figure 3.5: LAC diagram for three material decomposition (3MD)

(47)

3.5 Phantoms 35

Table 3.1: Effective energies and linear attenuation coefficients for the two spectra.

X-ray spectrum 80kV Sn140kV ρ Eef f [keV] 49.9 88.5

-µmE for water [1/cm] 0.2271 0.1774 1.00

µmE for protein [1/cm] 0.2815 0.2270 1.35

µmE for adipose [1/cm] 0.2079 0.1693 0.97

µmE for cort. bone [1/cm] 0.8223 0.3898 1.92

µmE for marrow [1/cm] 0.2193 0.1755 1.00

3.4.4 Chosing a threshold for tissue classification

The choice of a threshold value for the tissue classification methods, see Section 2.10, was done by taking the linear attenuation coefficients for protein at 80 kV spectrum and then use the topology to select the marrow as well. This can be done since the assumption that marrow only occurs inside bones has been made. For compensation of noise a slightly higher value than the one found in Table 3.1 were taken. As it turned out, for these experiments the choice was not very sensitive and could be done without further considerations about noise levels etc.

3.5 Phantoms

The phantom used in the evaluation process is a water phantom that has a radius of 9 cm. It contains five insertions each with a radius of 2 cm, one being a bone/marrow mixture and the other four mixtures of water, adipose and protein. The insertion containing bone tissue, R2, has a somewhat different geometry than

the others, it contains three sections, the outer section R2acontaining 100% bone,

the middle section R2bcontaining a mixture of 50% cortical bone and 50% marrow

and the inner most R2c containing 100% marrow. The marrow tissue is actually

a mixture of 50% red marrow and 50% yellow marrow. The phantom is described in detail in Fig. 3.7 and Tables 3.2-3.3.

In the later studies, Section 4.3.1, more advanced and realistic phantoms will be used. Phantom 2 is a simpler version of the well known Shepp-Logan phantom and is described in Fig. 3.8 and Tables 3.4-3.5. Phantom 3 is a phantom trying to resemble the pelvis, with two bones and some different soft tissue mixtures in between. It is described in Fig. 3.9 and Tables 3.6-3.7.

The densities of the bone tissue mixtures are calculated according to

ρ =m V = m m1 ρ1 + m2 ρ2 = _w 1 1 ρ1 + w2 ρ2 , (3.18)

where w1 and w2 is the mass fraction for cortical bone and marrow, respectively

(48)

Figure 3.7: Regions of phantom 1, sizes in Table 3.2.

Table 3.2: Sizes and positions of regions in phantom 1

Region center x [mm] center y [mm] radius [mm]

R1 0 0 90 R2a 0 60 13-20 R2b 0 60 7-13 R2c 0 60 7 R3 57.1 18.5 20 R4 35.3 -48.5 20 R5 -35.3 -48.5 20 R6 -57.1 18.5 20

Table 3.3: Mass fractions [%] of individual components and density, ρ in [g/cm3_]

of the mixtures for the regions of phantom 1

Region Cort. bone Marrow Water Protein Adipose ρ

R1 0 0 100 0 0 -R2a 100 0 0 0 0 1.92 R2b 50 50 0 0 0 1.32 R2c 0 100 0 0 0 1.00 R3 0 0 0 75 25 -R4 0 0 25 0 75 -R5 0 0 30 40 30 -R6 0 0 50 50 0

(49)

-3.5 Phantoms 37

Figure 3.8: Regions of phantom 2, sizes in Table 3.4.

Region center x [mm] center y [mm] radius [mm] angle [degrees]

R1 0 0 75/90 0 R2 0 0 65/80-70/85 0 R3 0 40 20/25 0 R4 -28.8 0 13.4/38.4 -70 R5 28.8 0 13.4/38.4 70 R6 -15 -50 10 0 R7 15 -50 10 0

Table 3.5: Mass fractions [%] of individual components and density, ρ in [g/cm3] of the mixtures for the regions of phantom 2

R1 0 0 100 0 0 -R2 100 0 0 0 0 1.92 R3 0 0 0 75 25 -R4 0 0 0 0 100 -R5 0 0 0 25 75 -R6 0 0 30 40 30 -R7 0 0 0 100 0

(50)

-Figure 3.9: Regions of phantom 3, sizes in Table 3.6.

Region center x [mm] center y [mm] radius [mm] angle [degrees]

R1 0 0 95/70 0 R2a ±60 -15 13-20 0 R2b ±60 -15 7-13 0 R2c ±60 -15 7 0 R3 25 0 15/40 -70 R4 -25 0 15/40 70 R5 0 40 20/10 0 R6 0 -50 20 0

Table 3.7: Mass fractions [%] of individual components and density, ρ in [g/cm3] of the mixtures for the regions of phantom 3

R1 0 0 100 0 0 -R2a 100 0 0 0 0 1.92 R2b 50 50 0 0 0 1.32 R2c 0 100 0 0 0 1.00 R3 0 0 30 40 30 -R4 0 0 0 100 0 -R5 0 0 0 0 100 -R6 0 0 0 75 25

An iterative reconstruction algorithm for quantitative tissue decomposition using DECT

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

An iterative reconstruction algorithm for

quantitative tissue decomposition using DECT

An iterative reconstruction algorithm for

quantitative tissue decomposition using DECT

Examensarbete utfört i Datorseende

vid Tekniska högskolan i Linköping

av

Abstract

Sammanfattning

Acknowledgments

Contents

List of Figures

Chapter 1

Introduction

1.1

Purpose and goal

1.2

Problem formulation

1.3

Approach

1.4

Related research

Chapter 2

Theory

2.1

General projection generation

2.2

Projection generation in CT

2.2.1

Monochromatic projections

2.2.2

Polychromatic projections

2.3

Backprojection

2.4

The Fourier-Slice Theorem

2.5

Filtered backprojection

2.5.1

Parallel-beam filtered backprojection

2.5.2

Fan-beam filtered backprojection

2.6

Fan-beam to parallel-beam rebinning

2.7

Forward projection

2.8

Monochromatic and polychromatic projections

2.9

Partial volume

2.10

Tissue classification

2.10.1

Simple threshold

2.10.2

Mixture of gaussians

2.11

Beam hardening

2.12

Three material decomposition

2.13

Two material plus density decomposition

Chapter 3

Methods

3.1

Simple water beam hardening correction

3.2

The iterative algorithm

3.2.1

Extensions to the 2MD plus density

3.3

Simple water beam hardening correction and

iterative algorithm

3.4

Experimental setup

3.4.1