Master of Science Thesis in Biomedical Engineering
Department of Electrical Engineering, Link¨oping University, 2021
Implementation of the Weighted Filtered
Backpro-jection Algorithm in the Dual-Energy Iterative
Al-gorithm DIRA-3D
Master of Science Thesis in Biomedical Engineering
Implementation of the Weighted Filtered Backprojection Algorithm in the Dual-Energy Iterative Algorithm DIRA-3D
Markus Tuvesson LiTH-ISY-EX--21/5358--SE
Supervisor: Alexandr Malusek
HMV, Link¨opings universitet
Examiner: Maria Magnusson
ISY, Link¨opings universitet
Department of Electrical Engineering Link¨oping University
SE-581 83 Link¨oping, Sweden February 19, 2021
Abstract
DIRA-3D is an iterative model-based reconstruction method for dual-energy helical CT whose goal is to determine the material composi-tion of the patient from accurate linear attenuacomposi-tion coefficients (LACs). Possible applications are, for example, to aid in calculations of radia-tion transport and dose calcularadia-tions in brachytherapy with low energy photons, and in proton therapy. There was a need to replace the cur-rent image reconstruction method, the PI-method, with a weighted fil-tered backprojection (wFBP) algorithm for image reconstruction, since wFBP is used for image reconstruction in Siemens’s CT-scanners. The new DIRA-3D algorithm implemented the program take for cone-beam projection generation and the FreeCT wFBP algorithm for image recon-struction. Experiments showed that the accuracies of the resulting LACs for the DIRA-3D algorithm using wFBP for image reconstruction were comparable to the one using the PI-method for image reconstruction. The relative LAC errors reached a value below 0.2 % after 10 iterations.
Acknowledgments
I want to thank my supervisor, Alexandr Malusek, and my examiner, Maria Magnusson, for their support during this project. A thanks also goes out to the IT-department at LiU for their help during the FreeCT WFBP installation. I also want to thank my mother and father for their support during my education at LiU, and my sister Josefine Tuvesson for her help with the proof-reading of this thesis report.
CONTENTS CONTENTS Contents 1 Introduction 6 1.1 Motivation . . . 6 1.2 Aim . . . 6 1.3 Problem Formulation . . . 6 1.4 Limitations . . . 7 2 Background 8 2.1 Photon Attenuation . . . 8 2.2 Helical CT-Scanner . . . 9 2.3 Dual-Energy CT . . . 11 2.4 Filtered Backprojection . . . 13 2.5 Iterative Reconstruction . . . 14 2.6 CT Artifacts . . . 14 2.6.1 Beam Hardening . . . 14
2.6.2 Partial Volume Artifacts . . . 15
2.6.3 Cone-Beam Artifacts . . . 16 2.6.4 Approximate-Reconstruction Artifacts . . . 16 3 Theory 17 3.1 DIRA-3D PI-Method . . . 17 3.2 DIRA-3D wFBP . . . 18 3.3 PI-Method . . . 19
3.4 Material Decomposition in DIRA-3D . . . 21
3.4.1 Two-Material Decomposition . . . 22
3.4.2 Three-Material Decomposition . . . 22
3.5 Forward Projection Generation in DIRA-3D . . . 22
3.6 Weighted Filtered Backprojection . . . 24
CONTENTS CONTENTS
4 Method 30
4.1 Phantom Definition . . . 30
4.2 Projection Generation . . . 31
4.3 Simulation of Measured Projections . . . 32
4.4 Forward Projections . . . 33
4.5 Reconstruction with FreeCT . . . 33
4.5.1 Row Direction Filtering . . . 34
4.6 Tissue Classification . . . 35
4.7 Error Calculation . . . 35
4.8 Computer Specifications . . . 36
4.9 Code Optimization . . . 36
5 Results 37 5.1 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Hanning Filter . . . 37
5.2 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Ram-Lak Filter . . . 40
5.3 DIRA-3D wFBP Reconstruction with Pitch 44.16 mm and Hanning Filter . . . 44
5.4 Time Results for DIRA-3D . . . 47
5.5 Dose Efficiency . . . 48
5.5.1 PI-Method Dose Efficiency . . . 48
6 Discussion 51 6.1 Results . . . 51
6.1.1 Comparison with PI-Method Results . . . 51
6.2 Method . . . 54
6.2.1 Phantom . . . 54
6.2.2 Source Criticism . . . 55
6.3.1 Parallel Computing . . . 55 6.3.2 Evaluating the DIRA-3D wFBP Algorithm with
a Higher Resolution . . . 56 6.3.3 Evaluating DIRA-3D wFBP for Different Q Values 57 6.3.4 The Long Object Problem . . . 57 6.4 Ethical and Societal Considerations . . . 57
7 Conclusions 59
References 60
A Appendix 63
1 INTRODUCTION
1 Introduction 1.1 Motivation
DIRA-3D is an iterative model-based reconstruction algorithm for dual-energy computed tomography (DECT), used to determine material com-position of an object from linear attenuation coefficients (LACs). These estimations might be used to aid calculations of radiation transport and dose calculations in for example brachytherapy, with low energy pho-tons, and in proton therapy. DIRA-3D is also able to produce mono-energetic images. For reconstruction of CT-images, DIRA-3D uses the PI-method, however since weighted filtered backprojection (wFBP) is used by Siemens in their CT-scanners, there is a need for the DIRA-3D algorithm to use wFBP for image reconstruction.
1.2 Aim
The aim of this thesis is to replace the PI-method, used for image re-construction in the DIRA-3D framework, with a wFBP algorithm and evaluate it.
1.3 Problem Formulation
This thesis main goal is to implement a wFBP algorithm for image re-construction in the DIRA-3D framework. The main changes that needs to be made to DIRA-3D are:
• Replace the current PI-method with the wFBP algorithm.
• Change how projection data is generated, since the current projec-tion data generaprojec-tion is specific to the PI-method.
This thesis will also try to answer the following questions:
• Can the wFBP algorithm, described by Stierstorfer et al. [1] and used in the FreeCT Project [2], replace the PI-method in DIRA-3D?
1.4 Limitations 1 INTRODUCTION
• The wFBP algorithm in FreeCT uses the GPU to speed up the algorithm. Can this be used effectively in our algorithm?
• How does DIRA-3D fair with the new wFBP algorithm compared to the previous PI-method algorithm? The comparison should be made in both image quality1, speed and dose efficiency2.
• Are there any other advantages or disadvantages with the new wFBP algorithm compared to the previous PI-method, for DIRA-3D? • Are there any parameters besides scanner geometry that affect the
reconstruction, and if so, in what way?
• Are there any possible improvements that can be done to the DIRA-3D algorithm?
This report aims to follow the directions given by the electrical engineer-ing department (ISY) at www.isy.liu.se/edu/xjobb/anvisnengineer-ingar_ exjobbare.html for the report layout.
1.4 Limitations
• A voxelized phantom with a low resolution was used to save com-putational time.
• A relatively small detector array and number of projection angles were used in order to save memory and computational time.
2 BACKGROUND
2 Background
2.1 Photon Attenuation
For photon beam attenuation, there are five interactions of photons with matter to take into consideration: photoelectric effect, Compton effect, Rayleigh scattering, pair production and photonuclear interactions [3, p. 125]. Of these five, Compton effect and photoelectric effect dominate in the diagnostic CT energy range of 15–150 keV [4, p. 41]. Rayleigh scattering, also known as coherent scattering, is a process where a pho-ton is redirected by an atom [3, p. 153]. The process is elastic so the energy loss of the photon is minimal and the redirection angle for the photon is small, causing the process to have only a small contribution to the overall photon attenuation in the considered energy range. Pair production and photonuclear interactions require higher photon ener-gies than those produced in diagnostic CT, so they will be disregarded in this report. Pair production have an energy threshold above 1 MeV and photonuclear interactions first occur with a photon energy of a few MeV [3, p. 146-147, 154].
Photoelectric effect is a process in which a photon is absorbed by a bound electron, causing the electron to depart from the atom with a kinetic energy
T = hv − Eb, (1)
where hv is the photon energy, Planck’s constant multiplied with the frequency of the photon, and Eb is the potential (binding) energy of the
bound electron [3, p. 138]. The chance that a photoelectric effect occurs is proportional to the macroscopic photoelectric effect cross-section Σph
[3, p. 140].
Compton effect, also known as incoherent scattering, is a process where a photon collides with an, assumed stationary and unbound, electron,
2.2 Helical CT-Scanner 2 BACKGROUND
transmitting a part of its energy to the electron, causing it to depart with a kinetic energy
T = hv − hv0, (2)
where hv0 is the photon’s energy after the collision [3, p. 126]. The chance that a Compton effect occurs is proportional to the macroscopic Compton effect cross-section ΣCo [3, p. 132-133].
For primary mono-energetic photons, the X-ray intensity I after travel-ing a distance x in a material can be calculated as
I = I0e−µx, (3)
where I0 is the initial X-ray intensity and µ is the linear attenuation
coefficient (LAC) for the material [4, p. 46]. The LAC describes the chance that a photon is either scattered or absorbed in the matter and is calculated as
µ = Σph+ ΣCo + ΣRa, (4)
where ΣRa is the macroscopic cross section for Rayleigh scattering. In
radiology, it is common to use a notation where µCo = ΣCo, µph = Σph,
and µRa = ΣRa. The mass attenuation coefficient, µm = µ/ρ, can then
be written as µm = µ ρ = µCo ρ + µph ρ + µRa ρ , (5)
where ρ is the density of the matter [3, p. 154] [5].
2.2 Helical CT-Scanner
In a helical CT-scanner, a gantry with an X-ray source and detector array, opposite one another, is rotating around an object to collect cone-beam projections for imaging of the object [4, p. 312-319]. The object can be described as a 3D function of µ(x, y, z). Projections, line inte-grals, are measured on this function. The cone-beam projections are
2.2 Helical CT-Scanner 2 BACKGROUND
is used to reconstruct the image. Putting the object on a table that moves through the gantry produces helical projections of the object, hence helical (spiral) CT. The distance covered by the table during one 360◦ gantry rotation is referred to as the pitch. The detector array con-sist of detectors organized into rows in the z-direction and into channels in the XY -plane (see figure 1).
z ζ fanangle γ Source R height
Figure 1: A cylindrical CT detector of the size 8×4 (channels × rows). In reality, the number of channels are much larger (see for example ref [6, p. 8] were the SOMATOM go.Sim have 920 channels). Image source: [7]. Used with permission from author.
The acquisition field of view (FOV) is defined by the fan-beam, see figure 2.
2.3 Dual-Energy CT 2 BACKGROUND
Figure 2: The max (acquisition) FOV of a fan-beam angle.
The reconstruction FOV is the size of the reconstructed image, which is smaller or equal to the size of the acquisition FOV.
2.3 Dual-Energy CT
In dual-energy CT (DECT), the dependence of the photoelectric effect and the Compton effect on the photon energy and the atomic number Z of the material is used to differentiate between different materials [8]. The Rayleigh effect also makes a small contribution [5]. This is possible by using two different tube voltages or by separating high and low X-ray energies [8]. It should be noted that the X-ray spectra are not monoenergetic, and the averages of photon energies are different for the two different tube voltages (see figure 10, section 4.3). The benefits of dual-energy information are that it is possible to differentiate and characterize different materials, which could lead to better pre-planning for operations and better diagnostics of diseases, injuries and afflictions.
2.3 Dual-Energy CT 2 BACKGROUND
There are several solutions on how to acquire the dual energy for mate-rial differentiation:
• Dual Spiral performs two scans subsequently with low and high tube voltage.
• TwinBeam has a split filter at the X-ray tube to generate two energy levels, separating low and high energies.
• Dual Source has two source-detectors in the same imaging plane at an angle 90◦. Both sources simultaneously operate at two different tube voltages.
• Dual Layer Detector uses a dual layer detector, for example the first layer absorbs low energies and the second layer absorbs high energies.
• Fast kV Switching alternates between low and high tube voltage at each projection.
• Slow kV Switching switches tube voltages between two gantry rotations in a single scan.
Moreover, dual-energy imaging can also be performed with photon count-ing systems with detectors able to resolve X-ray photon energies [9]. This technology allows for increased spatial resolution, reduced radia-tion exposure and the correcradia-tion of beam hardening artifacts compared to regular CT.
The above solutions have different applications and performances. The performance is evaluated on spectral separation, temporal coherence, temporal resolution, spatial resolution and dose efficiency [8]. Spectral separation is how good a modality can separate the two different energy level spectra. Temporal coherence denotes the time delay between the acquisition of the low and the high energy data. Temporal resolution
2.4 Filtered Backprojection 2 BACKGROUND
denotes the time necessary for collecting the data for image reconstruc-tion. Spatial resolution is related to the sharpness of the image and dose efficiency is about the dose needed to achieve desired results.
Disregarding photon counting systems, fast kV switching is, according to Faby et al. [8], an inferior solution in nearly all ways, while dual source is either similar or superior compared to the other solutions, if a tin-filter is used to remove low energy photons. Photon counting systems are currently available as experimental prototypes only.
2.4 Filtered Backprojection
Backprojection is a method of image reconstruction that smeare (back-project) the measured projections onto the image for each projection angle [10, p. 63-65]. For parallel projections in 2D it can be described as
fbp(x, y) =
Z
pθ(x cos θ + y sin θ)dθ, (6)
where pθ are the measured projections at projection angle θ. This does
however produce rather unsharp images and in order to fix that pθ(t)
is Fourier transformed into Pθ(ω) and is multiplied with a Ram-Lak
filter (ramp filter), or a window or function multiplied with a Ram-Lak filter. Finally, inverse Fourier transformation is used to obtain the filtered backprojection, which mathematically can be described as:
ff bp(x, y) = Z qθ(x cos θ + y sin θ)dθ, (7) where qθ(t) = Z Pθ(ω)|ω|ei2πωtdω. (8)
Pθ(ω) is the Fourier transform of a projection at a constant angle θ and
|ω| is the frequency response of the ramp filter, which means that qθ(t)
2.5 Iterative Reconstruction 2 BACKGROUND
only using simple backprojection. The filtered backprojection ff bp(x, y)
corresponds to the LAC values.
2.5 Iterative Reconstruction
In CT, the projection data is measured while the CT image (with LAC values) is unknown, which is why it needs to be reconstructed. With an estimation of the CT image, forward projection can be used to compare the estimation with the projection data [4, p. 357-358]. The difference between the forward projection and the measured projections is referred to as the error matrix, which is used to update the next image. The goal with each iteration is to reduce the error matrix. A final image will be the result after a certain number of iterations or when the error matrix is small enough.
There are several benefits of iterative reconstruction compared to regular FBP, for example it can produce images with higher signal to noise-ratio (SNR) at the same dose or produce images with similar SNR as FBP with lower dose. It can also model physical parameters that FBP cannot, for example the X-ray spectrum. Iterative reconstruction also allows for beam hardening artifacts reduction [11].
2.6 CT Artifacts
There are several types of CT artifacts. This section will describe a few that will be relevant to this thesis.
2.6.1 Beam Hardening
Beam hardening happens when a polyenergetic X-ray beam passes through an object causing the mean energy of the beam to increase, i.e. the beam becomes ”harder” [12]. This occurs due to the low energy photons being absorbed faster than the high energy photons. There are two
manifes-2.6 CT Artifacts 2 BACKGROUND
tations of the beam hardening: cupping artifacts, and streak or dark band artifacts.
Cupping artifacts appear as a darkening in the center of an object com-pared to the object’s outer regions (see figure 3); the intensity has a U-shape profile.
Dark bands and streaks can appear between two dense objects in very heterogeneous cross sections. This happens because the portion of the beam that passes through one of the objects at certain tube positions is hardened less than when it passes through both objects at other tube positions (see figure 3).
Figure 3: Beam hardening artifacts (bottom row) from different geometries (top row). Image source: [13] which is licensed under CC-BY.
2.6.2 Partial Volume Artifacts
Partial volume artifacts occur when a voxel contains more than one material (tissue), causing the attenuation value of the voxel to be the average between the materials in the voxel [12]. These artifacts can be reduced by increasing the acquisition resolution. To decrease par-tial volume artifacts below the size of the voxel size, the reconstruction resolution also needs to be increased.
2.6 CT Artifacts 2 BACKGROUND
2.6.3 Cone-Beam Artifacts
The use of a large cone beam angle may lead to undersampling in the cone-angle dimension, leading to cone-beam artifacts [4, p. 369-370]. Cone-beam artifacts tends to blur the image in an irregular way.
2.6.4 Approximate-Reconstruction Artifacts
Neither the PI-method nor wFBP are exact reconstruction methods [11, 1]; the inaccuracies may give rise to artifacts. These artifacts will be referred to as approximate-reconstruction artifacts in this report, as they were in the paper by Magnusson et al. [11].
It should be mentioned that an exact filtered backprojection method, the Katsevich reconstruction algorithm [14], does exist. Its implementation was investigated by Noo et al. [15].
3 THEORY
3 Theory
3.1 DIRA-3D PI-Method
DIRA-3D [11] is an extension of the 2D-algorithm DIRA developed by Malusek et al. [5]. It is a model-based iterative image reconstruction al-gorithm that uses DECT projections and automatic quantitative tissue classification for the generation of the reconstructed volume [5]. The polyenergetic projections acquired from the DECT-scanner are itera-tively converted into monoenergetic projections. The corrections are calculated by simulating monoenergetic and polyenergetic projections using the reconstructed volume for each of the tube voltages.
The algorithm is presented in figure 4. It performs the following steps [5, 11]:
1. Acquire measured semi-parallel polyenergetic projections, PM,Ui, for
two different tube voltages, Ui, i = 1, 2.
2. Acquire reconstructed volumes, µi = ∆µi+ µPi, at energy Ei. ∆µi is
the result from reconstructing PM,Ui − PUi with the PI-method (see
section 3.3). PUi are semi-parallel polyenergetic forward projections
and µPi is the reconstructed LAC from parallel monoenergetic
for-ward projection at energy Ei. During the first iteration µi is the
reconstruction of PM,Ui only, since PUi = PEi = 0 for the first
itera-tion.
3. Segment µi into pre-defined tissues and organs.
4. Acquire base material reconstructed volumes µC after material
de-composition, see section 3.4.
5. Calculate polyenergetic forward projections, PUi, for semi-parallel
3.2 DIRA-3D wFBP 3 THEORY
6. Calculate monoenergetic forward projections, PEi, for parallel
ge-ometry, see equation (21) in section 3.5.
Points 2-6 are repeated a predefined number of times. The final result is the two monoenergetic reconstructions µ1 and µ2, and the base material
reconstruction volumes µC. 1 2. µ 2 2. µ C 4. µ 5. P 6. PE2 U2 5. P 6. PE1 U1 ∆µ1 ∆µ P1 µ 2 µP2 3. µ 3. µ 1 2 material reconstr. volumes decom− threshold position material parallel FBP base
polyenergetic forward projection monoenergetic forward projection polyenergetic forward projection monoenergetic forward projection
1. PM,U2 measured projections 1. PM,U1 parallel FBP PI−method PI−method segmenta− tion
Figure 4: A flowchart of the DIRA-3D PI-Method algorithm. Image source: [11] which is licensed under CC BY 3.0.
The DIRA-3D PI-Method algorithm assumes that measured cone-beam helical projections are rebinned to semi-parallel projections first. For the sake of simplicity, this step was omitted in this implementation of DIRA-3D; the measured projections were generated as semi-parallel projections directly.
3.2 DIRA-3D wFBP
In DIRA-3D wFBP the PI-method was replaced by the FreeCT wFBP algorithm (see section 3.7) for reconstruction and the previous method for projection generation was replaced with the program take (see sec-tion 3.8). take generates cone-beam projecsec-tions [7, p. 21] which means
3.3 PI-Method 3 THEORY
that measured projections PM,Ui and polyenergetic forward projections
PUi now are cone-beam projections and not semi-parallel projections.
The new flowchart is shown in figure 5.
1 2. µ 2 2. µ C 4. µ 5. P 6. PE2 U2 5. P 6. PE1 U1 ∆µ1 ∆µ P1 µ 2 µP2 3. µ 3. µ 1 2 material reconstr. volumes decom− threshold position material parallel FBP base
polyenergetic forward projection monoenergetic forward projection polyenergetic forward projection monoenergetic forward projection
1. PM,U2 measured projections 1. PM,U1 parallel FBP segmenta− tion wFBP wFBP
Figure 5: A flowchart of the new DIRA-3D wFBP algorithm. Image source: [11] which is licensed under CC BY 3.0 (The image has been modified).
3.3 PI-Method
Note: The geometry and the notations are the same as for the wFBP-algorithm in section 3.6 with the sole exception for the parameter s, that is normalized in the wFBP-algorithm.
For the PI-method a virtual detector, the Tam window, is located at the isocenter, aligned with the z-axis of a CT-scanner (see figure 6). The Tam window ensures that the projection data is complete and non-redundant. [16, p. 99] [11]
3.3 PI-Method 3 THEORY
Figure 6: The rebinned geometry for the PI-method window, with the PI-detector (Tam window) aligned with the z-axis. Image source: [11] which is licensed under CC BY 3.0.
The rows of the PI-detector projected onto the cylindrical detector array can be described by the function:
q(s, γ) = s + γ P π cos γ , s ∈ −P 4, P 4 , (9)
where q is the height of the row in the specific coordinate, P is the pitch, γmax is the fan beam angle and s is the distance in the z-direction of a
specific row from the middle of the detector when γ = 0 [16, p. 101]. A simple step-wise explanation of the PI-method, which is used in DIRA-3D, is [16, p. 105]:
1. Perform rebinning pP(θ, t, s) = p(β, γ, s), where θ = β + γ and t = R sin γ (see figure 7). Cone-beam projections p(β, γ, s) are thus converted into semi-parallel projections pP(θ, t, s).
2. Discard projection data that is outside the Tam window. 3. Perform pre-weighting with the cone angle cos κ.
3.4 Material Decomposition in DIRA-3D 3 THEORY
5. Perform three-dimensional backprojection.
The PI-method is not an exact method, one reason being that the row-wise ramp filtering is one-dimensional [16, p. 105]. This causes small artifacts (see approximate reconstruction artifacts, section 2.6.4) in the reconstructed images that tend to grow in size with increased cone-beam angle, κmax [16, p. 105] [11].
3.4 Material Decomposition in DIRA-3D
The material decomposition in DIRA-3D comprises two-material and three-material decomposition [5]. Each can be written as a system of linear equations, which in the matrix form are
Aw = b, (10)
where A is the system matrix, w is the mass fraction vector and b = (0, 0, 1)T is the right-hand side vector. The elements in system matrix A are known while the elements in the mass fraction vector w are unknown. The mass fraction vector can be obtained with
w = A−1b. (11)
The system of linear equations (10) is derived from the summation rule for the mass attenuation coefficient of a mixture and the conservation of partial volumes in a mixture. The former states
µm(E) =
X
i
wiµm,i(E), (12)
where µm,i(E) is the mass attenuation coefficient at photon energy E,
and wi is the mass fraction for the ith base material,
P
iwi = 1. µm(E)
is thus the mass attenuation coefficient for the mixture of materials at photon energy E. The latter states
3.5 Forward Projection Generation in DIRA-3D 3 THEORY
where ρi is the mass density of the ith material and ρ is the mass density
of the mixture of materials.
3.4.1 Two-Material Decomposition
The system matrix, A, of the two-material decomposition can be written as µm,1(E1) µm,2(E1) −µ(E1) µm,1(E2) µm,2(E2) −µ(E2) 1 1 0 , (14)
where µ(E1) and µ(E2) are the LACs of the mixture at photon energy E1
and E2 [5]. The mass fraction vector is written as w = (w1, w2, 1/ρ)T,
where w1 and w2 are mass fractions and ρ is the mass density of the
mixture. The two-material decomposition is for instance used for the decomposition of bone tissue to compact bone and bone marrow.
3.4.2 Three-Material Decomposition
The system matrix, A, of the three-material decomposition can be writ-ten as µ(E1) ρ1 − µm,1(E1) µ(E1) ρ2 − µm,2(E1) µ(E1) ρ3 − µm,3(E1) µ(E2) ρ1 − µm,1(E2) µ(E2) ρ2 − µm,2(E2) µ(E2) ρ3 − µm,3(E2) 1 1 1 , (15)
and the mass fraction vector is written as w = (w1, w2, w3)T, where w1,
w2 and w3 are mass fractions [5]. The three-material decomposition is
for instance used for the decomposition of soft tissues to lipid, protein, and water.
3.5 Forward Projection Generation in DIRA-3D
The polyenergetic projection is calculated as P = ln I0
3.5 Forward Projection Generation in DIRA-3D 3 THEORY
where I is the beam intensity at the detector position when the attenu-ating object is present, and I0 is the same intensity when the attenuating
object is absent [5]. I0 is calculated as
I0 =
Z Emax
0
EN (E)dE, (17)
where E is the photon energy and N (E) is the energy spectrum of photons emitted from the X-ray tube. I is calculated as
I = Z Emax 0 EN (E) exp − Z L µ(x, y, z, E)dl dE, (18)
where µ(x, y, z, E) is the LAC of voxel (x, y, z) at energy E and RLdl is a line integral through the object [11].
By using material decomposition the line integralRLdl can be calculated through volume fractions of different base materials. The intensity I is then given by:
I = Z Emax 0 EN (E) exp " −X k µk(E)lk # dE, (19)
where µk is the LAC of the kth base material and lk is computed as
lk = ρ−1k
Z
L
ρ(x, y, z)wk(x, y, z)dl, (20)
where ρk is the density of the kth base material, ρ(x, y, z) is the
calcu-lated density in voxel (x, y, z) and wk(x, y, z) is the mass fraction of the
kth material in voxel (x, y, z). See two- and three-material decomposi-tion (secdecomposi-tions 3.4.1 and 3.4.2 respectively) on how the density ρ and the mass fraction wk are obtained.
The monoenergetic projection is calculated as PEi =
X
k
3.6 Weighted Filtered Backprojection 3 THEORY
equals the energy fluence weighted mean LAC for water
µw(Ei) = REmax 0 EN (E)µw(E)dE REmax 0 EN (E)dE . (22)
3.6 Weighted Filtered Backprojection
Weighted filtered backprojection (wFBP) is an iterative reconstruction algorithm for multislice spiral CT. The algorithm is the same as in the paper by Stierstorfer et al. [1] for a non-arbitrary pitch, but the notation has been changed to be the same as in the paper by Magnusson et al. [11] (see figure 7).
Figure 7 represents the geometry of a dual-source DECT-scanner. In the figure, κ is the cone-angle, γ is the fan-angle and P represents the pitch of the helix [11]. κmax is the maximum value of the cone angle
and is called cone-beam angle. γmax is the maximum value of the fan
angle and is called fan-beam angle. The s-axis is aligned with the z-axis and starts at the cross-section C between the central ray and the z-axis. The blue and red spherical markers represents the X-ray sources. The detectors are cylindrical, with axes parallel to the z-axis, and they travel in unison with their corresponding X-ray source in a helical trajectory around the z-axis. The helix path is given by:
x = R cos β, y = R sin β, z = P β
2π , (23)
3.6 Weighted Filtered Backprojection 3 THEORY
Figure 7: Dual-source helical CT geometry. Image source: [11] which is licensed under CC BY 3.0.
The wFBP algorithm consists of the following steps [1]: 1. Perform row-wise rebinning, p(θ, t, s) = p(β, γ, s) where
θ = β + γ, t = R sin γ, s ∈ [−1, 1]. (24)
Consequently, the rebinning converts the cone-beam projections into semi-parallel projections.
2. Filter along the row direction t with a kernel K (ramp filter), pconv(θ, t, s) =
Z
dt0p(θ, t0, s)K(t − t0). (25)
3. Backprojection is performed by calculating the following sum over all half-turns k and all detector rows s, for each voxel (x, y, z) and for each angle ˜θ ∈ [0, π[:
Vx,y,z(˜θ) = 1 Hx,y,z(˜θ) X k,s WQ(s)pconv(˜θ + kπ, ˆt, s) (26)
3.6 Weighted Filtered Backprojection 3 THEORY
4. The final voxel value is calculated by summing over ˜θ, Vx,y,z(θ) =
X
˜ θ
Vx,y,z(˜θ). (27)
In equation (26) the normalizing weight sum Hx,y,z is given by
Hx,y,z(˜θ) =
X
k
WQ(ˆs), (28)
and the relations between the detector coordinates ˆt, ˆs and the voxel coordinates (x, y, z) are given by
ˆ t = x sin θ − y cos θ (29) and ˆ s = z − P (˜θ + kπ − arcsin ˆ t R)/(2π)
(pR2 − ˆt2 − x cos θ − y sin θ) tan(κ
max/2)
. (30)
The weight function is given by
WQ(s) = 1 f or |s| < Q, cos2(π2|s|−Q1−Q ) f or Q ≤ |s| < 1, 0 otherwise, (31) where 0 < Q < 1.
The dose efficiency of the wFBP algorithm can be calculated as Dose efficiency = 2(1 + Q)
2
3 + 5Q . (32)
A Q value of 0.7 then gives a dose efficiency of 0.889, i.e. 11.1 % of the data is lost in the reconstruction. The dose usage is independent of the pitch. wFBP reconstruction is a non-exact method, causing low-frequency shade-like cone-beam artifacts in the reconstructed image (see approximate reconstruction artifacts, section 2.6.4). A dose efficiency of 1 is obtained for Q = 1, but has shown to give artifacts in the recon-structed image [1].
3.7 FreeCT Project 3 THEORY
3.7 FreeCT Project
The FreeCT Project is an open-source software project started at UCLA [2]. It offers a GPU-focused implementation of weighted filtered back-projection, FreeCT WFBP, for the reconstruction of 3rd generation he-lical CT projection data. FreeCT WFBP is based on algorithms de-scribed in the papers by Stierstorfer et al. [1] and Flohr et al. [17]. It also draws inspiration from Zinßer et al. [18]. It is written as a C program and uses the NVIDIA CUDA framework for the GPU execu-tion. The FreeCT wFBP algorithm reconstruction can also be executed without GPU, but according to Hoffman et al. [19] the CPU version is approximately 250–350 times slower for the scanner geometry and reconstruction parameters given.
FreeCT WFBP uses two configuration files. A geometry of the scanner is given in the scanner file, see table 1, and reconstruction parameters are given in the parameter file, see table 2.
FreeCT WFBP is executed via the command fct_wfbp prm/samples/filename.prm
in the FreeCT WFBP directory when using the GPU and with fct_wfbp --no-gpu prm/samples/filename.prm
when not using the GPU [20, p. 11]. The algorithm requires a raw data file with cone-beam projections as input and produces an output file with the reconstructed volume [20, p. 9-10].
3.7 FreeCT Project 3 THEORY
Table 1: Parameters in the scanner file [2]. Lengths are in millimeters, angles are in radians.
.
RSrcToIso Distance between X-ray source and isocenter.
RSrcToDet Distance between X-ray source and the detector.
AnodeAngle The anode angle of the X-ray source.1
FanAngleInc The increment angle between two adjacent channels.
ThetaCone The cone-angle κmax.
CentralChannel The center of the detector described in channels.
NProjTurn Number of projections per gantry rotation.
NChannels Number of channels in the detector.
ReverseRowInterleave 2
ReverseChanInterleave 2
1 The anode angle is the angle of the anode in the X-ray tube [4, p. 181-183]. 2
Used for focal flying spot, but they have no significance in this thesis.
Table 2: Parameters in the prm file. [20, p. 12] RawDataDir Directory of the raw data file.
RawDataFile Name of the raw data file.
OutputDir Output directory.
Nrows Number of rows in the detector.
CollSlicewidth Spacing between detector rows.
StarPos Location of the first slice to be reconstructed. EndPos Location of the last slice to be reconstructed.
Can be lower, higher or equal to the start position. SliceThickness Reconstructed slice thickness.
(Voxel size in the z-direction) TubeStartAngle (in whole degrees)
PitchValue Table travel per gantry rotation.
AcqFOV Acquisition field of view (see section 2.2). ReconFOV Reconstruction field of view.
ReconKernel Reconstruction kernel K (see equation (25)) selection. Readings Total number of projections contained in the raw data file.
Xorigin Spatial location of the horizontal reconstruction center in millimeters. Yorigin Spatial location of the horizontal reconstruction center in millimeters. Zffs Flag to FreeCT wFBP, to tell if Z flying focal spot was used.
Phiffs Flag to FreeCT wFBP, to tell if Phi flying focal spot was used.
Scanner Name of the scanner file.
FileType Type of data in the raw data file. FileSubType Not currently in use.
RawOffset Byte offset to the projection data in the raw data file.
Nx Width of the reconstruction volume in pixels. Should be a multiple of 32. Ny Depth of the reconstruction volume in pixels. Should be a multiple of 32.
3.8 Projection Generation Code take 3 THEORY
FreeCT WFBP also comes with a MATLAB file, ramp.m, that allows the user to specify a defined filter for the reconstruction kernel K (see equation (25)) in the wFBP algorithm.
An overview and performance testing of FreeCT WFBP is described in a paper by Hoffman et al. [19].
3.8 Projection Generation Code take
The projection generation code take is a MATLAB/C program for the simulation of X-ray projections from 3D volume data [7, p. 4]. The 2D X-ray projections can then be used as input for 3D reconstruction algorithms. take can use a polychromatic X-ray source, but it does not use energy-dependent detector sensitivity.
All parameters to the program are given in a number of text-files [7]: • take.txt is used as the initialization file.
• trj.txt describes the trajectory of the X-ray source and detector. • phm.txt describes the phantom as a mathematical object or link to
a voxel volume to use as a phantom. • det.txt describes the detector.
• spm.txt describes the energy-spectrum of the X-ray source.
• material files.txt describes the material in the phantom. This is only used for mathematical phantoms.
take can be executed as a stand-alone C program where the output data is then returned in a binary file, or as a MATLAB command where the out-data is returned as a MATLAB variable.
4 METHOD
4 Method
Experiments, with the DIRA-3D wFBP, were done for two different pitch values, P = 88.32 mm and P = 44.16 mm. For pitch 88.32 mm a Hanning and a Ram-Lak filter (see section 4.5.1) were used, while for pitch 44.16 mm only the Hanning filter was used.
air 800 400 slice 10 slice 14 slice 26 200 600 x water protein z, rotation axis lipid compact bone
0
slice 30 403 mm
muscle adipose tion areareconstruc−
Figure 8: (a) Phantom slice at z = 10. (b) The scanning geometry of the cone-beam projection generation though the phantom. The pitch is 88.32 mm and both drawings are to scale. Image source: [11] which is licensed under CC BY 3.0 (The image has been modified).
4.1 Phantom Definition
The phantom used for the experiments is the same as in the paper by Magnusson et al. [11]. It consist of six ellipsoids divided into two layers, one layer of four ellipsoids and one layer of two ellipsoids (see figures 8 and 9). The phantom had a resolution of 128 × 128 × 48 voxels and a voxel size of ∆x = ∆y = ∆z = 2.76 mm. The six ellipsoids had a diameter of 70.66 mm, 70.66 mm and 35.33 mm and they had their centers 105.98 mm from the isocenter. The phantom was loaded into MATLAB and all values above zero were set to 1 to ensure that its edges were sharp. A 3D-image of the phantom is shown in figure 9.
4.2 Projection Generation 4 METHOD
Figure 9: A 3D-image of the defined phantom before the ellipsoids were assigned any material.
Each ellipsoid in the phantom was assigned a different defined material, see table 3 and figure 8 for the material definitions used and to see what ellipsoid corresponds with what material.
Table 3: Elemental composition in mass fractions, corresponding density and tabulated LAC values µ(E1) and µ(E2), at E1= 50.0 keV and E2= 88.5 keV, for
the materials in the six ellipsoids.
Material Mass fraction Density µ(E1) µ(E2)
(%) (g/cm3) (m−1) (m−1) Adipose H 11.4, C 58.8, N 0.8, O 28.7, 0.95 20.2 16.6 tissue Na 0.1, S 0.1, Cl 0.1 Muscle H 10.2, C 14.2, N 3.4, O 71.1, 1.05 23.8 18.5 Na 0.1, P 0.2, S 0.3, Cl 0.1, K 0.4 Lipid H 11.8, C 77.3, O 10.9 0.92 19.1 16.0 Compact bone H 3.6, C 15.9, N 4.2, O 44.8, 1.92 79.2 38.7 Na 0.3, Mg 0.2, P 9.4, S 0.3, Ca 21.3 Protein H 6.6, C 53.4, N 17.0, O 22.0, 1.35 28.1 22.7 S 1.0 Water H 11.2, O 88.8 1.00 22.7 17.7 4.2 Projection Generation
Projection generation was done with the help of the program take (see section 3.8). A detector size of N × N = 192 × 64, a voxel size ∆x =
4.3 Simulation of Measured Projections 4 METHOD
angle κmax = 6.26◦ was used. See table 4 for the number of projection
angles, the number of helix turns and the pitch of the helix used in the experiments. The projection geometry is shown in figure 8 for the experiment with pitch P = 88.32 mm.
Table 4: Cone-beam projection parameters for the experiments. Projection angles Helix turns Pitch of helix
Experiments P = 88.32 800 2 88.32 mm (32 voxels)
Experiment P = 44.16 1600 4 44.16 mm (16 voxels)
4.3 Simulation of Measured Projections
The simulated measured projections of the phantom, described in sec-tion 4.1, were acquired for each ellipsoid separately to work around a limitation of take; the voxel phantom may consist of one material only (see Appendix A for the txt-files used for take). Consequently, take com-puted cone-beam line integrals, which were later converted into polyen-ergetic projections by using equations (16) and (19) in section 3.5. The energy spectra used are shown in figure 10.
0 50 100 150 energy (keV) 0 2 4 6 8 10 12 14 16
number of photons per bin
106 Energy Spectra
Figure 10: The photon energy spectra used in the DIRA-3D algorithm for the X-ray tube voltages 80 and 140 kV.
4.4 Forward Projections 4 METHOD
4.4 Forward Projections
3D forward projection generation was done with the help of the program take (see section 3.8) where take calculated the cone-beam line integrals from the base material reconstructed volumes (see figure 5). The cone-beam line integrals were later converted into polyenergetic projections by using equations (16) and (19) in section 3.5. The monoprojections are were obtained the same way as in ref [11].
4.5 Reconstruction with FreeCT
For the cone-beam reconstruction, the open-source FreeCT wFBP pro-gram was used. DIRA-3D wFBP created a binary file, with cone-beam projections from take, for the FreeCT wFBP application to read as in-put. The application was called with the terminal command
fct_wfbp prm/samples/filename.prm ,
hence the GPU was used in the experiments. The reconstructed image was stored in an img-file for DIRA-3D wFBP to use.
The FreeCT wFBP code was downloaded from github.com/FreeCT 2020-03-05, and the only things edited for the experiments was the parameter file and the scanner file (see tables 8, 9 and 10 in Appendix B).
In order to reconstruct the whole phantom, due to restrictions in FreeCT WFBP, the projection data had to be padded with zeros in the z-direction. For
example a data set of 800 projections was padded with 402 and 398 zero-filled projections at the beginning and the end of the array, respectively. The padding before the projection data had two zeros more than the padding after the projection data in order to fix a misalignment issue. The misalignment issue was a result from how the FreeCT algorithm reads binary projection data, where the last projection in the projection
4.5 Reconstruction with FreeCT 4 METHOD
StartPos used also needed to be larger than the EndPos used (see pa-rameter file in Appendix B), which in turn causes the FreeCT algorithm to read the binary file from the end to the beginning. This caused the FreeCT algorithm to read the cone-beam projections from take with a projection angle in offset, so this offset was removed by adding a zero to the beginning of the projection data.
4.5.1 Row Direction Filtering
Two filters were tested for the filtering along the row direction t in the wFBP algorithm (see K in equation (25), section 3.6). The filters were created using the FreeCT accompanying ramp.m file. For the Hanning filter, the impulse response h[k] and the frequency response H[θ] are shown in figure 11 and for the Ram-Lak filter the impulse and frequency responses are shown in figure 12. The Hanning filter, in this report, is a Hanning window multiplied with a Ram-Lak filter. The variables used in ramp.m to create the respective filters are shown in tables 5 and 6.
Figure 11: Impulse response h[k] and frequency response H[θ] for the Hanning filter, with θ being the normalized frequency.
4.6 Tissue Classification 4 METHOD
Table 5: Parameters used in ramp.m to create the Hanning filter. k [−384 −383 ... 382 383]
ds 402.7392 · sin (0.4386/2)
c 1
a 0.5
Figure 12: Impulse response h[k] and frequency response H[θ] for the Ram-Lak filter, with θ being the normalized frequency.
Table 6: Parameters used in ramp.m to create the Ram-Lak filter. k [−384 −383 ... 382 383]
ds 402.7392 · sin (0.4386/2)
c 1
a 1
4.6 Tissue Classification
The whole reconstructed volume was segmented into bone, soft tissue and air. Bone was decomposed into bone and bone marrow (see sec-tion 3.4.1). Air was decomposed into lipid and water. Soft tissue was decomposed into lipid, protein and water (see section 3.4.2).
4.7 Error Calculation
4.8 Computer Specifications 4 METHOD
where µt is the tabulated value (see table 3) and ¯µ is the average of the
calculated LAC in a spherical region of interest (ROI). The ROI was defined as a sphere with a radius of one third of the evaluated ellipsoid radius in the XY -plane, and positioned in the center of the evaluated ellipsoid.
4.8 Computer Specifications
For the experiments a computer with an 8-core 3.70 GHz Intel(R) Xeon(R) W-2145 CPU, an NVIDIA GeForce GTX 1080 with 8 GiB GPU, and 16 GiB RAM was used.
4.9 Code Optimization
MATLAB parfor -loops were replaced with for -loops in the tissue de-composition section of the DIRA-3D main file prog128.m. This was done since the parfor -loops was making the decomposition calculations slower. A parfor -loop was replaced with a for -loop for the summation of ∆µ and µP (see figure 5), since it made the calculations faster.
5 RESULTS
5 Results
The ability of the new DIRA-3D algorithm to accurately reconstruct LAC values is demonstrated in section 5.1, where the pitch was set to 88.32 mm, and in section 5.3, where the pitch was set to 44.16 mm.
5.1 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Hanning Filter
Figure 13 shows how the reconstructed LAC for the different materials (see figure 8 and table 3) in slice 10 changes with iteration number for the energies E1 = 50.0 keV and E2 = 88.5 keV. The LAC for compact
bone at low energy is the only LAC that shows a notable improvement, changing from approximately 59.80 m−1 in iteration 1 to 79.14 m−1 in iteration 10, and to 79.30 m−1 in iteration 25.
Figure 13: Reconstructed LAC images of the phantom for iterations 1, 2, 3, 10 and 25, slice 10, photon energies E1 = 50.0 keV (top row) and E2 = 88.5 keV
(bottom row), pitch 88.32 mm and filtering with the Hanning filter (see section 4.5.1).
Figure 14 shows how the reconstructed LAC changes for protein and water (see figure 8 and table 3) in slice 26 with iteration number for the two energies. There is no clear visible improvement of the LACs in this
5.1 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Hanning Filter 5 RESULTS
Figure 14: Reconstructed LAC images of the phantom shown for iterations 1, 2, 3, 10 and 25, slice 26, photon energies E1= 50.0 keV (top row) and E2= 88.5 keV
(bottom row), pitch 88.32 mm and filtering with the Hanning filter (see section 4.5.1).
An improvement of LACs is, however, shown in figure 15, where the relative errors (see section 4.7) of LACs averaged over the regions of interest (ellipsoids of the different materials) are plotted against the iteration number. For low energy the relative LAC errors of the soft tissues reached below the threshold 0.2 % after 5 iterations and after 10 iterations for compact bone. At high energy the relative LAC errors of the soft tissues reached below the threshold 0.2 % after 3 iterations and after 7 iterations for compact bone.
Figure 15: The relative errors of the LAC in the six different ellipsoids as a function of the number of iterations for pitch 88.32 mm, filtering with the Hanning filter (see section 4.5.1) and photon energies (a) E1= 50.0 keV and (b) E2= 88.5 keV.
5.1 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Hanning Filter 5 RESULTS
The suppression of beam hardening artifacts for lipid and adipose tissue is shown in figure 16 for slice 10 at iterations 1 and 10 for the two photon energies. For adipose tissue the LAC changed from 19.27 to 20.21 m−1 and for lipid it changed from 18.11 to 19.10 m−1, for low energy.
Figure 16: Profiles through the adipose tissue and lipid ellipsoids at iterations 1 (blue curve) and 10 (brown curve), in slice 10 and x-coordinate 64 for pitch 88.32 mm (true values are plotted as a black dashed curve), the Hanning filter (see section 4.5.1), and photon energies (a) E1= 50.0 keV and (b) E2= 88.5 keV.
Figures 17 and 18 show both the suppression of beam hardening artifacts and approximate-reconstruction artifacts for slices 14 and 30. The LAC was restricted to the interval -1 m−1 to 1 m−1, since the approximate-reconstruction artifacts are mainly visible in the air region [11], where the LAC is close to zero. Note that the shape of the artifacts for low and high energies is different; the reason being that the X-ray sources are positioned 90◦ with respect to each other.
5.2 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Ram-Lak Filter 5 RESULTS
Figure 17: Suppression of beam hardening and approximate-reconstruction arti-facts. Images of reconstructed LAC (in m−1) in slice 14 for iterations 1, 2, 3, 10, and 25 and photon energies E1= 50.0 keV (top row) and E2= 88.5 keV (bottom
row), pitch 88.32 mm and filtering with the Hanning filter (see section 4.5.1).. The range of LAC values was restricted to the interval [−1, 1]m−1.
Figure 18: Suppression of beam hardening and approximate-reconstruction arti-facts. Images of reconstructed LAC (in m−1) in slice 30 for iterations 1, 2, 3, 10, and 25 and photon energies E1= 50.0 keV (top row) and E2= 88.5 keV (bottom
row), pitch 88.32 mm and filtering with the Hanning filter (see section 4.5.1).. The range of LAC values was restricted to the interval [−1, 1]m−1.
5.2 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Ram-Lak Filter
The same experiment, as in previous section, was done after changing the reconstruction kernel K from the Hanning filter to the Ram-Lak
5.2 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Ram-Lak Filter 5 RESULTS
filter (see figure 12 in section 4.5.1). The difference of the LAC values (see figure 19 and 20) are minimal compared to when the Hanning filter was used and the relative LAC errors (see figure 21) reached below the threshold 0.2 % after the same number of iterations for both energy spec-tra. The curves in the profiles (see figure 22) are less smooth compared to when the Hanning filter was used.
Figure 19: Reconstructed LAC images of the phantom shown for iterations 1, 2, 3 and 10, slice 10, photon energies E1 = 50.0 keV (top row) and E2 = 88.5 keV
(bottom row), pitch 88.32 mm and filtering with the Ram-Lak filter (see section 4.5.1).
Figure 20: Reconstructed LAC images of the phantom shown for iterations 1, 2, 3, and 10, slice 26, photon energies E1 = 50.0 keV (top row) and E2 = 88.5
5.2 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Ram-Lak Filter 5 RESULTS
Figure 21: The relative errors of the LAC in the six different ellipsoids as a function of the number of iterations for pitch 88.32 mm, filtering with the Ram-Lak filter (see section 4.5.1) and photon energies (a) E1= 50.0 keV and (b) E2= 88.5 keV.
Figure 22: Profiles through the adipose tissue and lipid ellipsoids at iterations 1 (blue curve) and 10 (red curve), in slice 10 and x-coordinate 64 for pitch 88.32 mm (true values are plotted as a black dashed curve), the Ram-Lak filter (see section 4.5.1) and photon energies (a) E1= 50.0 keV and (b) E2= 88.5 keV.
There are more apparent approximate-reconstruction artifacts in the air region (see figures 23 and 24) compared to when the Hanning filter was used (see figures 17 and 18).
5.2 DIRA-3D wFBP Reconstruction with Pitch 88.32 mm and Ram-Lak Filter 5 RESULTS
Figure 23: Suppression of beam hardening and approximate-reconstruction arti-facts. Images of reconstructed LAC (in m−1) in slice 14 for iterations 1, 2, 3, and 10, photon energies E1 = 50.0 keV (top row) and E2 = 88.5 keV (bottom row),
pitch 88.32 mm and filtering with the Ram-Lak filter (see section 4.5.1). The range of LAC values was restricted to the interval [−1, 1]m−1.
Figure 24: Suppression of beam hardening and approximate-reconstruction arti-facts. Images of reconstructed LAC (in m−1) in slice 30 for iterations 1, 2, 3, and 10, photon energies E1 = 50.0 keV (top row) and E2 = 88.5 keV (bottom
row), pitch 88.32 mm and filtering with the Ram-Lak filter (see section 4.5.1). The range of LAC values was restricted to the interval [−1, 1]m−1.
5.3 DIRA-3D wFBP Reconstruction with Pitch 44.16 mm and Hanning Filter 5 RESULTS
5.3 DIRA-3D wFBP Reconstruction with Pitch 44.16 mm and Hanning Filter
The figures in this section show the same thing as the figures in the previous sections, but for the experiment with the pitch set to 44.16 mm and the reconstruction kernel K set to the Hanning filter (see section 4.5.1). The LAC results shown in figures 25 and 26 are similar to the corresponding LAC results in the previous section and the profiles in fig-ure 28 is similar to the profiles in figfig-ure 16. One notable difference from the experiment with the pitch set to 88.32 mm is that for low energy the relative LAC errors of the soft tissues reached below the threshold 0.2 % after 6 iterations and after 10 iterations for compact bone (see figure 27). At high energy the relative LAC errors of the soft tissues reached below the threshold 0.2 % after 4 iterations and after 7 itera-tions for compact bone. Another notable difference is the reduction of approximate-reconstruction artifacts in the air region (compare figures 29 and 17, and figures 30 and 18).
Figure 25: Reconstructed LAC images of the phantom for iterations 1, 2, 3, and 10, slice 10, pitch 44.16 mm, filtering with the Hanning filter (see section 4.5.1) and photon energies E1= 50.0 keV (top row) and E2= 88.5 keV (bottom row).
5.3 DIRA-3D wFBP Reconstruction with Pitch 44.16 mm and Hanning Filter 5 RESULTS
Figure 26: Reconstructed LAC images of the phantom for iterations 1, 2, 3, and 10, slice 26, pitch 44.16 mm, filtering with the Hanning filter (see section 4.5.1) and photon energies E1= 50.0 keV (top row) and E2= 88.5 keV (bottom row).
Figure 27: The relative errors of the LAC in the six different ellipsoids as a function of the number of iterations for pitch 44.16 mm, filtering with the Hanning filter (see section 4.5.1) and photon energies (a) E1= 50.0 keV and (b) E2= 88.5 keV.
5.3 DIRA-3D wFBP Reconstruction with Pitch 44.16 mm and Hanning Filter 5 RESULTS
Figure 28: Profiles through the adipose tissue and lipid ellipsoids at iterations 1 (blue curve) and 10 (red curve), in slice 10 and x-coordinate 64 for pitch 44.16 mm (true values are plotted as a black dashed curve), the Hanning filter (see section 4.5.1) and photon energies (a) E1= 50.0 keV and (b) E2= 88.5 keV.
Figure 29: Suppression of beam hardening and approximate-reconstruction arti-facts. Images of reconstructed LAC (in m−1) in slice 14 for iterations 1, 2, 3, and 10, photon energies E1 = 50.0 keV (top row) and E2 = 88.5 keV (bottom row),
pitch 44.16 mm and filtering with the Hanning filter (see section 4.5.1). The range of LAC values was restricted to the interval [−1, 1]m−1.
5.4 Time Results for DIRA-3D 5 RESULTS
Figure 30: Suppression of beam hardening and approximate-reconstruction arti-facts. Images of reconstructed LAC (in m−1) in slice 30 for iterations 1, 2, 3, and
10, photon energies E1 = 50.0 keV (top row) and E2 = 88.5 keV (bottom row),
pitch 44.16 mm and filtering with the Hanning filter (see section 4.5.1). The range of LAC values was restricted to the interval [−1, 1]m−1.
5.4 Time Results for DIRA-3D
The averaged time results for the different parts of new DIRA-3D (see section 3.2) is shown in table 7 for pitches of 88.32 and 44.16 mm, along with averaged time results for the PI-method. The parameters used for the PI-method were the same as for the experiments with pitch 88.32 mm for wFBP, except for the detector size, which was 192 × 32 instead of 192 × 64, i.e. half the size. In the DIRA-3D PI-method implementa-tion, for simplicity, the projection data were sampled directly on the PI-detector, i.e. the Tam window (see figure 6), and for the DIRA-3D wFBP implementation the projection data were sampled on the cylin-drical detector (see figure 7). If the Tam window is projected on the cylindrical detector, it gets the curved shape between the red lines in figure 31, and the minimal detector for wFPB while maintaining the same size of the detector elements and thus the same pixel size, is the
5.5 Dose Efficiency 5 RESULTS
wFBP was used, see the rectangle between the green lines.
The time results for DIRA-3D with the PI-method are from before addi-tional optimization. Note that these numbers, obviously, are dependent on the computer used to perform the iterative reconstruction algorithm.
Table 7: Averaged time results for the different parts of DIRA-3D with pitches 88.32 mm and 44.16 mm for wFBP, and for the PI-method.
Pitch: 88.32 mm 44.16 mm PI-method
Reconstruction: 8.7 s 11.1 s 29.6 s
Parallel FBP: 0.9 s 0.9 s 1.7 s
Segmentation: 0.008 s 0.007 s 0.030 s
Material decomposition: 2.7 s 2.7 s 6.5 s
Monoenergetic forward projections: 5.6 s 5.7 s 5.6 s Polyenergetic forward projections: 303.0 s 607.3 s 53.0 s
Total (one iteration): 320.9 s 627.7 s 96.4 s
5.5 Dose Efficiency
The Q-value used in the FreeCT wFBP algorithm for the weight function WQ was Q = 0.6, which is the same Q-value used in the paper by
Stierstorfer et al. [1]. Equation (32) in section 3.6 then results in a dose efficiency of 85 %.
5.5.1 PI-Method Dose Efficiency
Figure 31: The area inside the green lines is the detector used for the wFBP, the area inside the blue lines is the area of the cylindrical detector for the PI-method that is hit by the X-ray beam Abeam, and the area inside the red lines is the Tam
5.5 Dose Efficiency 5 RESULTS
For the PI-method the area of the detector array (corresponding to the outer detector in the thesis by Turbell [16, p. 100]) hit by the incoming X-ray beam, Abeam (see figure 31), can be derived from equation (9).
The maximum z-coordinate of the top row of the PI-detector projected onto the cylindrical detector array is
q(P/4, γmax) = (P/4 + γmaxP/π)/ cos(γmax)
and the minimum z-coordinate of the bottom row of the PI-detector projected onto the cylindrical detector array is
q(−P/4, −γmax) = (−P/4 − γmaxP/π)/ cos(−γmax).
The angular interval of Abeam is [−γmax, γmax] and the arc length of an
angle γ is rγ, where r is the radius of the cylindrical detector. Conse-quently Abeam is given by
Abeam = Z γmax −γmax 2 P/4 + γmaxP/π cos(γmax) r dγ = 4rP γmax 1/4 + γmax/π cos(γmax) . (34)
The area between the top and bottom row projected on the detector from the PI detector, API (see figure 31), is also derived from equation (9)
by setting s = P/4 and integrating over the arc length of the detector array, API = Z γmax −γmax 2P/4 + γP/π cos γ rdγ = P r 2 Z γmax −γmax sec γ dγ + 2P r π Z γmax −γmax γ sec γ dγ = P r 2 ln tan π 4 + γmax 2 − ln tan π 4 − γmax 2 . (35)
The function γ sec γ is odd and so its definite integral in equation (35) equals 0. The dose efficiency, ηPI, is then given by
ηPI = API A = cos γmax ln tan π4 + γmax 2 − ln tan π4 − γmax 2 (2 + γ 8/π)γ .
5.5 Dose Efficiency 5 RESULTS
Note that the dose efficiency of the PI-method is pitch independent, but dependent on the fan-angle γmax. For γmax = 0.44 rad (see section 4.2),
6 DISCUSSION
6 Discussion 6.1 Results
In the pitch 88.32 mm experiments, where Hanning and Ram-Lak filters were compared for row direction filtering in the wFBP algorithm, the corresponding LAC results in figures 13 and 19, figures 14 and 20, and the relative LAC error in figures 15 and 21 are similar. It is however shown in the profile figures (see figures 16 and 22) that the use of the Hanning filter resulted in a smoother curve. The Hanning filter also reduced approximate-reconstruction artifacts in the air region (compare figure 17 with figure 23 and figure 18 with figure 24).
Comparing experiments with pitches 88.32 mm and 44.16 mm, the shorter pitch, for high photon energy 88.5 keV, required more iterations for the relative LAC errors to get below the 0.2 % threshold. This may how-ever not be due to the change in pitch, but caused by the iteratively updated reconstructed images being different. A shorter pitch produced less approximate-reconstruction artifacts in the air region (compare low energy figures 17 and 29). Halving the pitch doubled the number of cone-beam projections so a better reconstruction result was expected.
6.1.1 Comparison with PI-Method Results
After 10 iterations, the results for DIRA-3D PI-method in the paper by Magnusson et al. [11] and the results for the DIRA-3D wFBP recon-struction with the Hanning filter and pitch 88.32 mm were similar (see figure 32).
6.1 Results 6 DISCUSSION
Figure 32: DIRA-3D results after 10 iterations for both DIRA-3D PI-method (right) and DIRA-3D wFBP algorithm (left), and for the photon energy spectras E1 = 50.0 keV (top row) and E2 = 88.5 keV (bottom row). The range of LAC
values are restricted to the interval [−1, 1]m−1. The DIRA-3D wFBP results are for pitch 88.32 mm and the Hanning filter.
Relative LAC errors for the DIRA-3D PI-method and the DIRA-3D wFBP algorithms (see figure 33) got below the 0.2 % threshold after 10 iterations for the low energy E1 = 50.0 keV and after 7 iterations for the
6.1 Results 6 DISCUSSION
Figure 33: The relative errors of the LAC in the six different ellipsoids as a function of the number of iterations for the photon energies (a) E1 = 50.0 keV
and (b) E2 = 88.5 keV for both the DIRA-3D PI-method (bottom) and for the
DIRA-3D wFBP algorithm (top). The relative LAC error results for DIRA-3D wFBP are for pitch 88.32 mm and the Hanning filter.
While the wFBP code was faster than the PI-method code, the DIRA-3D wFBP iterations were slower. This was due to the cone-beam projection generation from take taking longer time than the semi-parallel projec-tion generaprojec-tion used in DIRA-3D PI-method. The wFBP code being faster than the PI-method code is likely due to the wFBP code utilizing the GPU while the PI-method did not. Optimizing the algorithm have however resulted in the parallel FBP, segmentation and material decom-position being faster in DIRA-3D wFBP than in DIRA-3D PI-method. While the simulations in this report have been done without noise, noise does occur in CT. The level of noise is determined by the number of incident photons collected by the detector, where less photons collected results in more noise [21, p. 4]. The wFBP algorithm used in the experiments has a higher dose efficiency than the PI-method, for the
6.2 Method 6 DISCUSSION
regarded was higher for the wFBP than for the PI-method resulting in a lower dose with the same SNR. How much of a dose reduction this provides is uncertain, but if the change results in a dose reduction of more than 10-15 % the wFBP algorithm should be considered above the PI-method [22].
For the PI-method the detector height and detector width, dependent on γmax, are dependent on the pitch while for the wFBP, detector height and
detector width are not dependent on the pitch, making the PI-method a less flexible reconstruction method.
6.2 Method
A possible issue with utilizing the FreeCT wFBP software may be that it may change, and thus compatibility issues may arise. When this thesis was written DIRA-3D and FreeCT wFBP were not fully integrated into each other, but instead DIRA-3D executed the FreeCT wFBP via a terminal command. This meant that the compatibility issues that may arise should be contained to directory issues, issues with the parameter and scanner files, and issues with the raw-file (binary file with cone-beam projections) and output-file (binary file with reconstructed image). A larger detector size (192 × 64 compared to 192 × 32) was used for DIRA-3D wFBP compared to DIRA-3D PI-method. This was done since it is possible and necessary to have a larger detector in wFBP, while the PI-method is restricted by the Tam window.
6.2.1 Phantom
The phantom used for the experiments was a relative simple one with low resolution. Further validation of the DIRA-3D wFBP algorithm with more complex phantoms is needed.
6.3 Possible Extensions or Improvements for DIRA-3D 6 DISCUSSION
6.2.2 Source Criticism
DIRA-3D PI-method is based on refs [5, 16] and its performance is demonstrated in ref [11]. The FreeCT wFBP algorithm is based on refs [17, 1] while being inspired by ref [18]. Both ref [5] and ref [17] are peer-reviewed and published in Medical Physics. Ref [11] is pub-lished and peer-reviewed in Biomedical Physics & Engineering Express, ref [1] is peer-reviewed and published in Physics in Medicine and Biol-ogy and ref [16] is published as a phD thesis at Link¨oping University. Ref [18] is published in conjunction with the 12th International Meet-ing on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine.
The FreeCT wFBP algorithm is evaluated in ref [19], which is peer-reviewed and published in Medical Physics.
The evaluation of the new DIRA-3D algorithm was done with the help of information from refs [5, 11, 12]. Ref [12] is published and peer-reviewed in RadioGraphics. The dose efficiency of the PI-method was derived according to information in refs [1, 16].
6.3 Possible Extensions or Improvements for DIRA-3D
6.3.1 Parallel Computing
At the time this thesis was written the cone-beam projection generation done by take was the most time-consuming part of DIRA-3D wFBP, so utilizing parallel computing like the function parfor in MATLAB or the GPU in take may shorten its calculation time.
It was tried to speed up the polyenergetic forward projection by com-bining low and high energy calculation of doublets (volumes from two-dimensional decomposition) and the calculation of low and high energy
6.3 Possible Extensions or Improvements for DIRA-3D 6 DISCUSSION
This did however have unexpected results and without an explanation on why it did not work, the idea was abandoned. Had the idea worked, the calculation of the polyenergetic forward projection would have taken approximately half the time, given that the computer could handle at least six parallel computations with parfor. It is however possible that this idea was correct, but its implementation was wrong, and it might be worth looking into it again.
The parfor -loops replaced in the material decomposition section (see section 4.9) might be faster than the for -loops if the phantom or the resolution are changed, i.e. if material decomposition functions’ calcu-lations become more time consuming.
Another way to shorten the computation time of the polyenergetic pro-jections might be to make take utilize the GPU for its calculations. It is doubtful that both of these methods can be utilized at the same time, so if the GPU is used by take, the current calls to take by DIRA-3D might need a rework. This is due to the calls to take being done in parfor -loops.
6.3.2 Evaluating the DIRA-3D wFBP Algorithm with a Higher Resolution
Current evaluations of DIRA-3D have been done with rather low acqui-sition and reconstruction resolution. Increasing these resolutions should diminish the pronunciation of the approximate-reconstruction artifacts, visible near compact bone in figures 17, 23 and 29, by making them thinner. Increasing the resolution also have the benefit of making the images clearer.
Decreasing the computational time of the polyenergetic forward projec-tions should be a priority over increasing the resolution, since increasing the acquisition and reconstruction resolution would increase the
compu-6.4 Ethical and Societal Considerations 6 DISCUSSION
tational time exponentially.
It is also likely that the amount of RAM-memory (16 GB RAM was used for this project) needs to be increased when the resolutions are increased.
6.3.3 Evaluating DIRA-3D wFBP for Different Q Values
According to Stierstorfer et al. [1] setting the Q value in WQ to 1 results
in artifacts. Evaluating DIRA-3D wFBP with Q set to 1, with the hope that the iterative reconstruction in DIRA-3D make the artifacts mentioned by Stierstorfer disappear, might be of interest. The benefits would be a dose efficiency of 100 %.
Changing the value of Q would be done by changing it in backproject.cuh if the GPU is used, and in backproject cpu.cu if it is not and re-install FreeCT WFBP. This was not done for this project due to the difficulties with installing FreeCT WFBP on the computer used.
6.3.4 The Long Object Problem
DIRA-3D PI-method have an algorithm that can handle the long object problem (see the paper by Magnusson et al. [11] for more information). The implementation was removed since it was custom made for the PI-method and could not easily be applied for DIRA-3D wFBP. A proposal for the future is to look at the long object problem again to solve it either directly in the FreeCT wFBP algorithm or just make it a part of DIRA-3D and compatible with the FreeCT wFBP algorithm.
6.4 Ethical and Societal Considerations
DIRA-3D is a proof-of-concept code that will not be used clinically. The algorithm may, however, be implemented in clinical CT-scanners
6.4 Ethical and Societal Considerations 6 DISCUSSION
yet implemented in DIRA-3D, the code could improve the accuracy of reconstructed CT numbers. The change from the PI-method to the wFBP algorithm increased the dose efficiency.
