Automatisch kalibreren van CT-scanners tijdens
iteratieve beeldreconstructie
Jim van Nunen
Promotoren: prof. dr. Roel Van Holen, prof. Christian Vanhove
Begeleider: dr. ir. Bert Vandeghinste
Masterproef ingediend tot het behalen van de academische graad van
Master of Science in de ingenieurswetenschappen: computerwetenschappen
Vakgroep Elektronica en Informatiesystemen
Voorzitter: prof. dr. ir. Jan Van Campenhout
Faculteit Ingenieurswetenschappen en Architectuur
Academiejaar 2013-2014
Automatisch kalibreren van CT-scanners tijdens
iteratieve beeldreconstructie
Jim van Nunen
Promotoren: prof. dr. Roel Van Holen, prof. Christian Vanhove
Begeleider: dr. ir. Bert Vandeghinste
Masterproef ingediend tot het behalen van de academische graad van
Master of Science in de ingenieurswetenschappen: computerwetenschappen
Vakgroep Elektronica en Informatiesystemen
Voorzitter: prof. dr. ir. Jan Van Campenhout
Faculteit Ingenieurswetenschappen en Architectuur
Academiejaar 2013-2014
Voorwoord
Er waren allerlei redenen waarom ik voor dit thesisonderwerp koos. Enerzijds denk ik dat een
computerwetenschapper zich niet moet beperken tot zijn eigen vakgebied, maar zijn vaardigheden
ook in andere onderzoeksdomeinen moet inzetten. In dit geval de biomedische wetenschappen,
een discipline waar ik aanvankelijk niet mee vertrouwd was. Het is anderzijds een belangrijke
eigenschap van een ingenieur om zich snel vertrouwd te maken met nieuwe kennis. Dit onderzoek
zou waarschijnlijk op een andere manier tot stand zijn gekomen mocht iemand met een achtergrond in de biomedische wetenschappen zich erover gebogen hebben, maar ik acht de confrontatie
van uiteenlopende ideeën van het grootste belang in de wetenschappelijke vooruitgang. Laat ik
hier op zijn minst een onderzoek presenteren waar op verdergebouwd kan worden.
Een voorwoord als dit verlangt traditioneel een hele lijst bedankjes. Laat ik het simpel houden
en het grootste woord van dank voorbehouden aan Bert Vandeghinste, mijn thesisbegeleider, die
mij wanneer nodig de juiste inzichten heeft bijgebracht en altijd klaarstond om mijn vragen te
beantwoorden. Een tweede woord van dank gaat naar Roel Van Holen en Christian Vanhove,
mijn promotors, die dit uitdagende onderwerp ter beschikking hebben gesteld. Voor de rest acht
ik een thesis geen plaats om lief, vrienden, familie, kat of hond in de kijker te zetten. Daar
zijn betere manieren voor en — laten we eerlijk wezen — de kans dat ze dit ooit gaan lezen, is
uitermate klein.
Er is toch een kleine uitzondering die ik wil maken. Door een drukke agenda, zowel binnen als
buiten de universiteit, moest het indienen van mijn thesis noodgedwongen wachten tot augustus. Door er vaart achter te zetten kon ik toch nog vroeger dan gepland werk maken van het
weerzien met een kennis in Sint-Petersburg. Het was Joris die er ondertussen voor zorgde dat
de papieren versies van de thesis op hun bestemming aankwamen. Dank voor deze praktische
beslommeringen.
Een motivatie en een dankwoord. Ontbreekt er nog iets? Een citaat wellicht. Naar aloude
gewoonte een stukje tekst waarvan het twijfelachtig is of de vermeende auteur er ooit iets mee
te maken heeft gehad. “A professor is someone who talks in someone else’s sleep.” W.H. Auden.
Voldoet aan voornoemd criterium. Sluit samen met deze thesis een vijfjarige ingenieursopleiding
af.
Jim van Nunen, 1 augustus 2014
Toelating tot bruikleen
De auteur geeft de toelating deze masterproef voor consultatie beschikbaar te stellen en delen
van de masterproef te kopiëren voor persoonlijk gebruik.
Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit
deze masterproef.
Jim van Nunen, 1 augusuts 2014
Automatisch kalibreren van CT-scanners
tijdens iteratieve beeldreconstructie
door
Jim van Nunen
Masterproef ingediend tot het behalen van de academische graad van
Master of Science in de ingenieurswetenschappen: computerwetenschappen
Academiejaar 2013–2014
Promotoren: prof. dr. Roel Van Holen, prof. Christian Vanhove
Begeleider: dr. ir. Bert Vandeghinste
Faculteit Ingenieurswetenschappen en Architectuur
Universiteit Gent
Vakgroep Elektronica en Informatiesystemen
Voorzitter: prof. dr. ir. Jan Van Campenhout
Samenvatting
In deze thesis wordt onderzocht hoe de geometrie van CT-scanners automatisch gekalibreerd
kan worden tijdens iteratieve beeldreconstructie. Een foutieve inschatting van de geometrische
parameters kan op reconstructies leiden tot wazige randen en een afname in spatiale resolutie.
Het kalibreren wordt geabstraheerd tot het oplossen van een optimalisatieprobleem. Als objectieven komen in de projectieruimte de projectiefout aan bod en in de beeldruimte de scherpte.
Beide maten worden vergeleken en er wordt bepaald voor welke parameters ze het best geschikt
zijn. Deze studie gebeurt op tweedimensionale reconstructies. De behaalde inzichten worden
uiteindelijk toegepast op de driedimensionale reconstructie van een geplastineerde muis, met
projectiedata afkomstig van de TriFoil Triumph-II X-O CT-scanner.
Trefwoorden
computertomografie, geometrische kalibratie, iteratieve beeldreconstructie, projectiefout, scherpte
Self-calibration of CT scanners during iterative
image reconstruction
Jim van Nunen
Supervisor(s): Bert Vandeghinste, Christian Vanhove, Roel Van Holen
Abstract—Computed tomography (CT) is a nondestructive imaging technique that uses reconstruction algorithms to create threedimensional tomographic images from scanner data. The reconstruction algorithms do not
always adhere to the geometry of the scanner. Misalignments in scanner
geometry are known to lead to edge blurring and a loss of spatial resolution. Especially in the context of micro-CT, known for its high spatial
resolution, effects of misalignments are a severe threat to the overall image
quality. We present an object-independent method to calibrate CT scanners
during iterative image reconstruction. Our calibration method is based on
the minimization or maximization of objective functions in projection and
image space. We compare both the reprojection error and sharpness and
determine which of these objectives is best suited for each geometric parameter. The proposed calibration scheme is able to optimize individual
parameters independently of other parameters, resulting in a considerable
speed-up compared to other calibration methods that simultaneously take
all geometric parameters into account.
Keywords— micro computed tomography, geometric calibration, iterative image reconstruction, reprojection error, sharpness
methods have exhibited some major flaws over time. Maximizing sharpness has been proven to outperform the minimization
of entropy: entropy ignores spatial information and indicates
segmentability instead of overall sharpness of the reconstructed
image.
The Nelder-Mead algorithm is often cited as an ideal algorithm to find the correct geometric parameters [6][8]. This optimization algorithm has good convergence properties, but only
if the objective functions are free of local extrema [9]. Sawall et
al. [10] propose a genetic algorithm for misalignment estimation to escape from such local extrema. Despite its success, the
algorithm still tries to optimize all geometric parameters simultaneously, thereby making the search space unnecessarily large.
I. I NTRODUCTION
Our goal is to determine which objective is best suited for
particular parameters. Moreover, to reduce the search space of
the optimization algorithms, we will determine the parameters
that can be estimated independently of other parameters. Our
initial research was carried out over 2D reconstructions of a selfdesigned phantom. Ultimately we present a stepwise calibration
scheme for 3D reconstructions. The geometry used is shown
in Figure 1. All data (scans of a plastinated mouse [11]) were
acquired with the TriFoil Triumph-II X-O CT scanner.
C
OMPUTED TOMOGRAPHY is a nondestructive imaging
technique that uses reconstruction algorithms to create 3D
tomographic images from scanner data. The reconstruction algorithms do not always adhere to the geometry of the scanner,
resulting in edge blurring and a loss in spatial resolution. Misalignments in scanner geometry are due to mechanical effects
(e.g. gravity) or changes in atmosphere (e.g. temperature). In
the context of micro-CT, known for its high spatial resolution,
these effects are a severe threat to the overall image quality.
Successful calibration methods have been proposed, but often lack efficiency or make specific assumptions about the used
phantoms.
A first method to calibrate CT scanners involves the examination of projection images of a specifically designed geometrical object [1]. As a consequence, this step should precede the
main CT experiment. A second class of methods does not rely
on dedicated calibration phantoms and is object-independent by
placing tracking markers on the object [2]. The accuracy of such
methods depends on the correct placement of the markers [3]. A
third class of methods, and the focus of this research, is to directly use the acquired projection images. During iterative image reconstruction, the quality of the reconstruction is numerically determined by objective functions in projection or image
space.
In projection space, the calculation of the reprojection error
is most common [4][5]. However, there is an increased trend
toward the use of objective functions in image space, since optima can exist for reconstructions that do not have the lowest
reprojection error possible. The work of Kyriakou et al. [6] is
entropy-based and has yielded promising results. Kingston et
al. [7] propose sharpness as objective since the entropy-based
II. G OAL
III. PARAMETERS
We define eight parameters in our model for 3D reconstructions. ∆x, ∆y and ∆z denote the offsets of the object, which
we consider stationary. ∆u and ∆v are the detector’s offsets. φ
is the detector twist angle. SO and SD denote the distance between the source and object and the distance between the source
and detector, respectively. Table I shows the parameter set Ω
which we consider an optimum because of the sufficient sharp-
Fig. 1: Model for 3D reconstruction.
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Ω
ˆ
Ω
0.00 mm
0.00 mm
0.00 mm
0.00 mm
0.00 mm
0.00 mm
-38.39 pixels
-30.00 pixels
-155.87 pixels
-162.00 pixels
-0.28◦
0.00◦
299.64 mm
293.00 mm
131.17 mm
140.00 mm
Ωcal
0.00 mm
0.00 mm
-4.22 mm
-37.81 pixels
-151.99 pixels
-0.26◦
305.69 mm
127.61 mm
TABLE I: Parameter sets for 3D reconstruction.
(k)
(k)
SΩ = k∇ˆfΩ k2 =
XX
x
Fig. 2: 3D reconstruction of the plastimouse.
ˆ from where to start the
ness and detail, and our initial guess Ω
optimization process. Figure 2 compares both reconstructions.
IV. O BJECTIVE FUNCTIONS
The object is mathematically defined as the image distribution
f . g is the data as measured by the scanner. The relationship
Hf = g holds, in which H denotes the system matrix. Since this
matrix is often too complex to calculate, iterative reconstruction
methods estimate f by using parameterized projection operators
(k)
H. We call an estimation of f after k iterations ˆf .
A. Reprojection error
(k)
The reprojection error EΩ is defined as the error between
(k)
the measured projection g and the data Hˆf projected by the
Ω
(k)
(k)
reconstruction algorithm. Note that ˆfΩ , and thus EΩ , depends
on both Ω and the number of iterations k. Formally [5]:
(k)
EΩ =
P
X
(k)
|gi − (HˆfΩ )i |
(2)
y
The image gradient can be calculated by using the differentiation property of the discrete Fourier transform, but a simple approximation can be obtained by using gradient convolution kernels such as the Sobel mask. Sharpness is an objective
function to be maximized in the image space by an optimization algorithm. Its usefulness will be limited for the optimization of SD and SO, because sharpness is strictly not scale in(k)
variant [7]: SΩ increases with increasing magnification factor
M = SO/SD.
ˆ (transverse)
(b) Ω
(a) Ω (transverse)
(k)
|∇fˆΩ (x, y)|2
(1)
i=1
(k)
gi , i ∈ {1, . . . , P }, is a single element of g. The form of EΩ is
closely related to the cost that SIRT (Simultaneous Iterative Reconstruction Technique) tries to minimize. Therefore, we will
make use of SIRT exclusively. The reprojection error is an objective function to be minimized in the projection space by an
optimization algorithm.
B. Sharpness
Sharp images generally contain the maximum amount of
high-frequency information. Blurred edges typically suppress
these high frequencies. A possible measure for sharpness is the
L2 norm of the image gradient [7]:
V. C ALIBRATION OF 3D RECONSTRUCTIONS
We propose the optimization scheme for 3D reconstructions
in Figure 3, based on a preliminary study of 2D reconstructions.
Due to the presence of local extrema, the genetic algorithm was
chosen over the Nelder-Mead procedure.
First, the parameters that affect the displacement of the detector in its plane are optimized by using sharpness. Note that
the calculation of sharpness is based on edge detection. Typical artefacts from a misaligned detector involve double edges.
These edges can be evaluated as ‘sharper’. An upper bound on
the reprojection error is thus needed to avoid false maxima due
to double edges. Sharpness leads toward the optimal parameters
a few iterations faster than the reprojection error (k = 2 was
chosen).
Note that we only optimized ∆u and φ. Misalignments of
both ∆v and ∆z can introduce artefacts in projection space, on
the edge(s) where the object crosses the detector plane in the vdirection. These artefacts can be suppressed by either ∆v or ∆z,
but only particular values for both ∆v and ∆z will ultimately
yield the sharpest reconstruction. They are optimized together
by considering both sharpness and reprojection error by means
of a MOEA (Multiple Objective Evolutionary Algorithm) [12].
The outcome will be a trade-off between a low reprojection error
(no artefacts in the projection space) and high sharpness in the
image space. Once the directions and regions of the parameters
are known, we discard sharpness and minimize the reprojection
error. k = 4 was sufficient.
The correct values of SO and SD can only be found by optimizing them both by means of the reprojection error. SO cannot
be estimated if SD is still suboptimal, and vice versa. Moreover,
a high number of iterations is needed, e.g. k = 100 or k = 200.
The search space consists of many local extrema, which makes
it hard to estimate these parameters, even for the genetic algorithm. The genetic algorithm was stopped before the termination
criteria were met, yielding a suboptimal result Ωcal for SO and
SD, compared to the optimum Ω (see Table I). However, the reconstruction was sharp and virtually indistinguishable from the
reconstruction with Ω (see Figure 4). The estimation was good
enough due to SO and SD having less impact on artefacts and
edge blurring.
VI. C ONCLUSIONS
Most calibration methods try to optimize all geometric parameters at once, thereby making the search space of the optimization algorithms unnecessarily large. We have shown that
certain geometric parameters can be optimized independently of
other parameters, which results in a much smaller search space
at a time and a faster calibration. Not all optimization steps
require as many iterations of the iterative reconstruction algorithm, in our case SIRT. Moreover, sharpness is not an equivalent of the reprojection error, since it has local extrema by
favouring double edges. However, if we can set an upper bound
to the reprojection error (which has no local extrema for the
parameters we can optimize by maximizing sharpness), maximizing sharpness can estimate the parameters a few iterations
faster than the reprojection error. Sharpness is also proven to
Fig. 3: Stepwise optimization scheme.
(a) Ω (transverse)
(c) Ω (sagittal)
(b) Ωcal (transverse)
(d) Ωcal (sagittal)
Fig. 4: Comparison of supposed optimum and estimation.
direct the parameters in the right direction, but the final result
will still be suboptimal. We made use of this observation when
designing the MOEA. The Nelder-Mead algorithm is only useful in absence of local extrema, a situation that seldom occurs
in the context of the calibration of CT scanners, especially when
optimizing all parameters at once. The genetic algorithm was
powerful enough to escape the various local extrema.
R EFERENCES
[1] C. Mennessier, R. Clackdoyle, and F. Noo, “Direct Determination of
Geometric Alignment Parameters for Cone-Beam Scanners”, Physics in
Medicine and Biology, vol. 54, 2009, pp. 1633-1660.
[2] I. S. Kyprianou, S. Paquerault, B. D. Galas, A. Badano, S. Park, and K.
J. Myers, “Framework for Determination of Geometric Parameters of a
Cone Beam CT Scanner for Measuring the System Response Function and
Improved Object Reconstruction”, Proceedings of the IEEE International
Symposium on Biomedical Imaging, 2006, pp. 1248-1251.
[3] F. Stopp, A. J. Wieckowski, M. K¨aseberg, S. Engel, F. Fehlhaber, and E.
Keeve, “A Geometric Calibration Method for an Open Cone-Beam CT
System”, Proceedings of the 12th International Meeting on Fully ThreeDimensional Image Reconstruction in Radiology and Nuclear Medicine,
2013, pp. 106-109.
[4] J. Wicklein, H. Kunze, W. A. Kalender, and Y. Kyriakou, “Image Features for Misalignment Correction in Medical Flat-Detector CT”, Medical
Physics, vol. 39, no. 8, 2012, pp. 4918-4931.
[5] W. Wein, A. Ladikos, and A. Baumgartner, “Self-Calibration of Geometric and Radiometric Parameters for Cone-Beam Computed Tomography”,
Proceedings of the 11th International Meeting on Fully Three-Dimensional
Image Reconstruction in Radiology and Nuclear Medicine, 2011, pp. 327330.
[6] Y. Kyriakou, R. M. Lapp, L. Hillebrand, D. Ertel, and W. A. Kalender,
“Simultaneous Misalignment Correction for Approximate Circular ConeBeam Computed Tomography”, Physics in Medicine and Biology, vol. 53,
2008, pp. 6267-6289.
[7] A. Kingston, A. Sakellariou, T. Varslot, G. Myers, and A. Sheppard, “Reliable Automatic Alignment of Tomographic Projection Data by Passive
Auto-Focus”, Medical Physics, vol. 38, no. 9, 2011, pp. 4934-4945.
[8] D. Panetta, N. Belcari, A. Del Guerra, and S. Moehrs, “An OptimizationBased Method for Geometrical Calibration in Cone-Beam CT without Dedicated Phantoms”, Physics in Medicine and Biology, vol. 53, 2008, pp.
3841-3861.
[9] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence
Properties of the Nelder-Mead Simplex Method in Low Dimensions”, SIAM
Journal on Optimization, vol. 9, no. 1, 1998.
[10] S. Sawall, M. Knaup, and M. Kachelrieß, “An Adaptive Genetic Algorithm
for Misalignment Estimation (AGAME) in Circular, Sequential and Spiral
Cone-Beam Micro-CT”, Proceedings of the 11th International Meeting on
Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear
Medicine, 2011, pp. 187-190.
[11] Frank Verhaegen, Maastro Clinic, The Netherlands.
[12] E. Zitzler, M. Laumanns, and L. Thiele, “SPEA2: Improving the Strength
Pareto Evolutionary Algorithm”, TIK-Report 103, 2001.
Inhoudsopgave
1 Inleiding
1
2 Computertomografie
3
2.1
Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.2
CT-scanner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.2.1
X-stralenbuis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.2.2
Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
Beeldreconstructie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.3.1
Analytisch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3.2
Iteratief . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.4
Micro-CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.5
Kalibreren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.5.1
Probleemstelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.5.2
Literatuur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.5.3
Doelstelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.5.4
Verloop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.3
3 Iteratieve beeldreconstructie
12
3.1
Beeld- en projectieruimte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
3.2
Lineair model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
3.3
Iteratieve reconstructie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
3.3.1
Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
3.3.2
Algoritmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
SIRT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
3.4
4 Reconstructie in twee dimensies
17
4.1
Geometrische parameters
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
4.2
Object . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
4.2.1
Definitie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
4.2.2
Parameterkeuze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
4.2.3
Reconstructies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
i
INHOUDSOPGAVE
ii
5 Optimalisatietechnieken
5.1
5.2
5.3
5.4
23
Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
5.1.1
Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
5.1.2
Objectieffuncties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
5.1.3
Optimalisatiealgoritmes . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
Objectieffuncties voor het CT-probleem . . . . . . . . . . . . . . . . . . . . . . .
26
5.2.1
Entropie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
5.2.2
Projectiefout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
5.2.3
Scherpte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
Optimalisatiealgoritmes voor het CT-probleem . . . . . . . . . . . . . . . . . . .
30
5.3.1
Nelder-Mead algoritme . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
5.3.2
Genetisch algoritme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
6 Minimalisatie projectiefout
6.1
37
Invloed parameters op projectiefout . . . . . . . . . . . . . . . . . . . . . . . . . .
37
6.1.1
Horizontale uitwijking detector ∆u . . . . . . . . . . . . . . . . . . . . . .
37
6.1.2
Verticale uitwijking detector ∆v
. . . . . . . . . . . . . . . . . . . . . . .
39
6.1.3
Detectortilt θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
6.1.4
Afstand stralingsbron tot object SO . . . . . . . . . . . . . . . . . . . . .
42
6.2
Gedrag projectiefout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
6.3
Optimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
6.3.1
Directe optimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
6.3.2
Stapsgewijze optimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . .
45
6.3.3
Aanpassingen genetisch algoritme . . . . . . . . . . . . . . . . . . . . . . .
46
Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
6.4
7 Maximalisatie scherpte
7.1
47
Invloed parameters op scherpte . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
7.1.1
Horizontale uitwijking detector ∆u . . . . . . . . . . . . . . . . . . . . . .
47
7.1.2
Verticale uitwijking detector ∆v
. . . . . . . . . . . . . . . . . . . . . . .
48
7.1.3
Detectortilt θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
7.1.4
Afstand stralingsbron tot object SO . . . . . . . . . . . . . . . . . . . . .
49
7.2
Gedrag scherpte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
7.3
Optimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
7.3.1
Directe optimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
7.3.2
Stapsgewijze optimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
7.4
INHOUDSOPGAVE
iii
8 Gecombineerde objectieffuncties
52
8.1
Probleemstelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
8.2
Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
8.2.1
Lineaire combinatie van objectieffuncties . . . . . . . . . . . . . . . . . . .
53
8.2.2
Pareto-efficiënte oplossingenverzameling . . . . . . . . . . . . . . . . . . .
54
Optimalisatiealgoritme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
8.3.1
Bruikbaarheidsfunctie . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
8.3.2
Creatie nieuwe populatie . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
8.4
Testresultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
8.5
Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
8.3
9 Reconstructie in drie dimensies
58
9.1
Geometrische parameters
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
9.2
Specificaties GPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
9.3
Voorwaartse en terugwaartse projectie . . . . . . . . . . . . . . . . . . . . . . . .
60
9.3.1
Stralingsbron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
9.3.2
Object . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
9.3.3
Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
Objectieffuncties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
9.4.1
Projectiefout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
9.4.2
Scherpte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
SIRT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
9.4
9.5
10 Kalibreren
65
10.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
10.2 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
10.3 Eerste tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
10.3.1 Scherpte maximaliseren . . . . . . . . . . . . . . . . . . . . . . . . . . . .
68
10.3.2 Projectiefout minimaliseren . . . . . . . . . . . . . . . . . . . . . . . . . .
70
10.3.3 Gecombineerde objectieffuncties . . . . . . . . . . . . . . . . . . . . . . . .
71
10.4 Invloed parameters in beeld- en projectieruimte . . . . . . . . . . . . . . . . . . .
71
10.4.1 Uitwijkingen van het object . . . . . . . . . . . . . . . . . . . . . . . . . .
71
10.4.2 Detectorparameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
10.4.3 Vergrotingsfactor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
10.5 Stappenplan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
76
10.5.1 Optimaliseren van ∆u en φ . . . . . . . . . . . . . . . . . . . . . . . . . .
76
10.5.2 Optimaliseren van ∆v en ∆z . . . . . . . . . . . . . . . . . . . . . . . . .
77
10.5.3 Optimaliseren van SD en SO . . . . . . . . . . . . . . . . . . . . . . . . .
77
10.6 Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
79
10.6.1 Herhaalbaarheid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
79
INHOUDSOPGAVE
iv
10.6.2 Nauwkeurigheid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
10.6.3 Tijd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
10.6.4 Samenvatting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
11 Conclusie
82
11.1 Resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
11.2 Toekomstig werk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
A Inhoud digitale drager
84
A.1 CPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
84
A.2 GPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
B Configuratiebestanden
87
B.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
87
B.2 Formaat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
87
C Gebruik optimalisatiealgoritme
90
C.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
C.2 Generieke structuur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
D Reconstructies
92
Gebruikte afkortigen en symbolen
CPU
Central Processing Unit
CT
Computed Tomography
FBP
Filtered Backprojection
GPU
Graphics Processing Unit
MLTR
Maximum Likelihood Transmission Reconstruction
MOEA
Multiple Objective Evolutionary Algorithm
MRI
Magnetic Resonance Imaging
OS-SART
Ordered-Subset Simultaneous Algebraic Reconstruction Technique
SIRT
Simultaneous Iterative Reconstruction Technique
SPEA
Strength Pareto Evolutionary Algorithm
f
objectdistributie
g
gemeten projectiedata
H
systeemmatrix (voorwaartse projectieoperator)
H−1
inverse systeemmatrix (terugwaartse projectieoperator)
α
rotatiehoek
∆x, ∆y, ∆z
uitwijkingen object in x-, y- en z-richting
∆u, ∆v
uitwijkingen detector in u- en v-richting
θ
detectortilt
φ
detectortwist
SD
afstand stralingsbron tot detector
SO
afstand stralingsbron tot object
M
vergrotingsfactor
Ω
parameterverzameling
p
parameter
c
objectieffunctie
Tt
tableau t van optimalisatiealgoritme
At
archief t van optimalisatiealgoritme
(k)
projectiefout na k iteraties van reconstructie met Ω
(k)
scherpte na k iteraties van reconstructie met Ω
EΩ
SΩ
Hoofdstuk 1
Inleiding
Computertomografie, afgekort tot CT van het Engelse computed tomography, levert twee- of
driedimensionale afbeeldingen of tomogrammen op die de inhoud van een driedimensionaal object
weergeven. Het object in kwestie wordt hierbij niet ‘opengemaakt’, waardoor men spreekt van
een zogenaamd niet-invasieve onderzoekstechniek. Het tomogram moet digitaal, aan de hand
van een algoritme, gereconstrueerd worden uit tweedimensionale stralingsdata. Verschillende
reconstructiealgoritmes kunnen toegepast worden op dezelfde data. Ze staan dus volledig los van
het eigenlijke scanproces, waardoor aanpassingen en verbeteringen mogelijk zijn zonder dat de
(dure) apparatuur op dezelfde ingrijpende wijze veranderingen moet ondergaan.
CT-scanners maken gebruik van x-stralen (ioniserende straling). Ze vinden voornamelijk hun
toepassing in de medische beeldvorming, waar ze een contrastrijker beeld opleveren dan de traditionele röntgenopname. Het grote voordeel is dat CT dwarsdoorsnedes van het object levert,
zoals te zien op figuur 1.1. In klinische beeldvorming zijn opnames van het hoofd, de longen,
de (onder)buik en het bekken het meest gebruikelijk. CT laat onder meer toe om tumoren
nauwkeurig op te sporen.
(a) klassieke röntgenopname
(b) CT-scan [1]
Figuur 1.1: Scans van de menselijke romp met verschillende technologieën.
1
HOOFDSTUK 1. INLEIDING
2
De CT-scanner heeft ook een intrede gedaan in de preklinische beeldvorming. Hier gaat het om
micro-CT met als kenmerken een hogere spatiale resolutie door een sterke vergrotingsfactor, om
weefsels, organen en botstructuren van kleine proefdieren in vivo of ex vivo te scannen. In de
preklinische beeldvorming kunnen CT-scans bijdragen tot biomedische experimenten, maar ook
tot het uitvoerig testen van nieuwe CT-technieken alvorens ze op grote schaal in te voeren.
De kwaliteit van een CT-scan is deels afhankelijk van de parameters waarmee de algoritmes het
beeld reconstrueren. Deze parameters beschrijven de geometrie van de gebruikte CT-scanner. De
werkelijke geometrie kan afwijken van de veronderstelde door bijvoorbeeld mechanische effecten.
Onscherpte aan de randen van de reconstructie of een afname in spatiale resolutie zijn het gevolg
en verlagen de algemene kwaliteit van het beeld, zoals te zien op figuur 1.2. Een te lage kwaliteit
betekent dat mogelijk cruciale details om een bepaald (ziekte)beeld te herkennen, onzichtbaar
blijven. De arts zal niet langer in staat zijn om met grote zekerheid een bepaalde diagnose te
bevestigen of pathologieën uit te sluiten.
(a) niet-gekalibreerd
(b) gekalibreerd
Figuur 1.2: Invloed van kalibreren op reconstructiekwaliteit [3].
Het bepalen van de geometrische parameters is mogelijk door periodiek CT-scans te nemen van
een object waarvan men de geometrische eigenschappen exact kent. Men weet dan welke reconstructie men mag verwachten om zo bepaalde parameters bij te sturen in het geval van onscherpe
of door artefacten aangetaste beelden. Het is echter ook mogelijk om het reconstructiealgoritme
zelf de geometrische paramaters te laten schatten om een zo optimaal mogelijk beeld te bekomen. Een objectieffunctie in functie van de geometrische parameters moet numeriek bepalen
in hoeverre een beeld optimaal is. Het kalibreren van CT-scanners kan met andere woorden
geabstraheerd worden tot het oplossen van een optimalisatieprobleem. Dit wiskundige concept
heeft reeds zijn nut bewezen in domeinen als de economie en de mechanica, niet in het minst
sinds de komst van de computer en bijhorende rekenkracht.
Deze thesis zal tot doel hebben een systeem te ontwikkelen waarmee een CT-scanner zichzelf
kan kalibreren, of beter: de gangbare reconstructiealgoritmes zullen automatisch de verzameling
parameters bepalen die een zo goed mogelijke reconstructie opleveren.
Na een korte uiteenzetting over de werking van computertomografie in het volgende hoofdstuk,
zullen we in sectie 2.5 de doelstellingen van dit onderzoek nauwkeurig afbakenen.
Hoofdstuk 2
Computertomografie
2.1
Principe
CT levert op niet-invasieve wijze twee- of driedimensionale afbeeldingen op van de inhoud van
een driedimensionaal object. Deze afbeeldingen zijn doorsnedes van het object en maken de
techniek dus uitermate geschikt voor medische toepassingen. Het waren Allan Cormack en
Godfrey Hounsfield die in 1979 samen de Nobelprijs voor de Fysiologie of Geneeskunde in de
wacht sleepten na hun ontwikkeling van de eerste CT-scanner. Toch waren de achterliggende
principes al minstens 60 jaar vroeger bekend, meer bepaald uit het werk van de wiskundige
Johann Radon. Het is een lot waarmee veel wiskundigen zich moeten verzoenen. Inmiddels is
computertomografie uitgegroeid tot een populaire en goed presterende scantechniek, al gaat de
innovatie onverminderd verder. Figuur 2.1 toont het uitzicht van een moderne CT-scanner.
Een CT-scanner bestaat ruwweg uit drie belangrijke componenten. Het object, in medische
context een patiënt of proefdier, rust op een tafel die zich tussen een x-stralenbuis en detector
bevindt. De x-stralenbuis is een stralingsbron die tijdens de scan rond het object roteert. De
detectoren meten de intensiteit van de uitgestuurde x-stralen nadat ze doorheen het object ge-
Figuur 2.1: Typisch ontwerp van een CT-scanner in medische context [2].
3
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
4
Figuur 2.2: Onderlinge positie van stralingsbron, detector en object.
propageerd zijn. Een reconstructiealgoritme zal van deze tweedimensionale, planaire data een
driedimensionale weergave maken. Voornoemde onderdelen en principes zullen nu verder toegelicht worden, in zoverre dat nodig is voor een goed begrip van de iteratieve reconstructiealgoritmes
die aan de basis van het kalibreren zullen liggen.
2.2
2.2.1
CT-scanner
X-stralenbuis
De x-stralenbuis is de bron van ioniserende straling die tijdens acquisitie in de CT-scanner roteert
rond het object. Het object ligt stil, terwijl de x-stralenbuis en detector zich recht over elkaar
bevinden en samen roteren rond een centrale as, zoals getoond op figuur 2.2.
X-straling of röntgenstraling is elektromagnetische straling in het energiebereik tussen 100 eV en
100 keV. In de x-stralenbuis wordt een elektromagnetisch veld tussen een anode en kathode aangelegd. De kathode zal onder invloed van een verhittingsproces elektronen vrijgeven die door het
elektromagnetisch veld een versnelling ondergaan. De elektronen produceren een remstralingsbundel bij inval op het wolfraam waaruit de anode (net als de kathode) meestal vervaardigd is.
Het proces gaat overigens gepaard met hitteproductie, omdat het slechts 1% van de elektronen
omzet in x-stralen.
Recente CT-scanners maken gebruik van een kegelvormige stralenbundel (cone beam) die het
hele object kan omvatten, hoewel dit niet zal lukken in humane studies. Op deze manier kunnen
meerdere doorsnedes tegelijk gemeten worden en is de stralingsdosis voor de patiënt kleiner. De
x-stralen hebben voldoende hoge energie om atomen te ioniseren en atoombindingen te verbreken.
Het effect ervan op levende weefsels kan tot uiting komen in stralingsziekte of kanker. Minstens
even belangrijk als de stralingsdosis is de tijd waarover men aan de dosis blootgesteld wordt.
Een Belg zou per jaar vier keer meer straling ondergaan dan een Nederlander of Brit [4]. Er zou
sprake zijn van te veel nutteloze CT-scans, waarop de FOD Volksgezondheid in maart 2014 liet
weten meer MRI-scans te willen aanmoedigen [5]. MRI staat voor magnetic resonance imaging en
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
5
Figuur 2.3: Radiologisch contrast door verschil in gedetecteerde stralingsintensiteit.
maakt in tegenstelling tot CT geen gebruik van ioniserende straling, maar van magneetvelden en
radiogolven. De techniek is meer geschikt voor zachte weefsels dan voor botstructuren, waardoor
beide technologieën niet onderling inwisselbaar zijn. Toepassingen ontwikkelen uit theoretisch
wetenschappelijk onderzoek kan nooit zonder aan het maatschappelijk debat deel te nemen.
2.2.2
Detector
De fysica achter de productie van x-stralen is verder niet belangrijk voor het beoogde kalibreren:
er zullen immers geen wijzigingen aangebracht worden aan de sturing of hardware van de CTscanner. Inzicht in de interactie van de x-stralen met het object is daarentegen wel van cruciaal
belang om te begrijpen welke data de detectoren meten en opslaan. Deze data zullen samen met
de geometrische parameters de invoer van het reconstructiealgoritme vormen.
Het medium waardoor x-stralen propageren, zal een deel van deze straling absorberen. Voor
mono-energetische x-stralen geldt voor absorptie de wet van Lambert-Beer:
I = I0 e−
RL
0
µ(s) ds
(2.1)
Vergelijking 2.1 legt het verband tussen de intensiteit I0 van de uitgestuurde x-stralenbundel en
de intensiteit I van de bundel waargenomen door de detector, na transmissie door het object.
L is de dikte van dit object met s een bepaalde plaats in het medium. De absorptiecoëfficiënt
µ is materiaalafhankelijk. Gezien het object uit meerdere materialen (bijvoorbeeld weefsels of
botstructuren) kan bestaan, is µ afhankelijk van de plaats s in het object. De absorptiecoëfficiënt
is de probabiliteit dat een foton verdwijnt per eenheid van laagdikte en is verder recht evenredig
met de massadichtheid. Figuur 2.3 toont het principe van radiologisch contrast door detectie
van verschillende stralingsintensiteiten naargelang de aard van het gepenetreerde medium. Het
contrast op de gereconstrueerde afbeelding ontstaat uit dergelijke verschillen in waargenomen
intensiteit.
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
6
Figuur 2.4: Venster bepaalt welke CT-nummers op de het grijswaardebereik afgebeeld worden.
Er moet opgemerkt worden dat vergelijking 2.1 enkel geldig is voor mono-energetische x-stralen.
In werkelijkheid is het spectrum van de uitgestuurde x-stralen polyenergetisch, waardoor aanpassingen aan de uitdrukking van de absorptiecoëfficiënt µ nodig zijn. Het heeft voor het verder
verloop van dit onderzoek geen zin om deze aanpassingen in detail te bespreken, maar het volstaat op te merken dat de absorptiecoëfficiënt zal afnemen bij verhoogde energie. Fotonen met
hoge energie zullen het medium gemakkelijker penetreren met minder ruis tot gevolg, maar tegen
de prijs van een lager contrast.
Net omdat bij conventionele radiografie en -scopie een afbeelding van de tweedimensionale beeldreceptor meteen het eindresultaat is, kunnen enkel structuren die sterk verschillen in absorptiecoëfficiënt bij gegeven x-stralenkwaliteit goed gevisualiseerd worden. Bovendien scheiden dergelijke systemen niet de ruimtelijke structuren in de richting van de x-stralenbundel, waardoor
diepte ontbreekt. CT komt aan beide problemen tegemoet door vanuit verschillende hoeken
tweedimensionale data te registreren en ze bij reconstructie terug samen te voegen tot een driedimensionaal geheel.
2.3
Beeldreconstructie
Een structuur met gemiddelde absorptiecoëfficiënt µ stemt overeen met een bepaald CT-nummer
of Hounsfieldeenheid HUµ volgens onderstaande formule:
HUµ = 1000 ×
µ − µwater
µwater
(2.2)
Het resultaat van een CT-scan is echter een tomogram met grijswaarden. De Hounsfieldeenheden, die waarden van -1000 tot 3000 kunnen aannemen, moeten afgebeeld worden op een 8 bit
grijswaardebereik (8 bit: 28 = 256 verschillende waarden). Een vensterbreedte moet bepalen
welke CT-nummers overeenkomen met het grijswaardebereik op het weergavescherm. Een klein
venster levert hoog contrast op: er zijn veel grijswaarden te kiezen voor een klein bereik. Een
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
7
groot venster levert laag contrast op: hetzelfde aantal grijswaarden moet nu meer CT-nummers
visualiseren. Figuur 2.4 illustreert het vensterprincipe voor 256 grijswaarden. Merk op dat het
hier niet om het intrinsieke contrast gaat, maar om de loutere weergave op een scherm.
De reconstructie kan analytisch of iteratief tot stand komen.
2.3.1
Analytisch
Beeldreconstructie werd aanvankelijk beschouwd als het inverteren van de discrete Radontransformatie, een probleem dat filtered backprojection (FBP) exact oplost [6]. De Radontransformatie
zelf stelt het projectieproces voor. De projectie van het object levert een zogenaamd sinogram
op, een twee- of driedimensionale dataverzameling die verschillende projecties van het gescande
object bevat, namelijk een projectie per hoek vanwaar de x-stralen het object gepenetreerd hebben. FBP veronderstelt dat de detectordata exacte lijnintegralen van de objectdistributie zijn.
Onder meer door de beperkte bruikbaarheid van vergelijking 2.1 als model om het werkelijke
stralingsgedrag te beschrijven, is analytische beeldreconstructie slechts een vereenvoudigd model
dat tot ongewenste beeldartefacten kan leiden.
2.3.2
Iteratief
Een andere klasse van reconstructiealgoritmes zijn de iteratieve technieken. Zij gaan uit van een
rijker (al dan niet lineair) model dat de voorwaartse projectie voorstelt. Het model gaat beter
om met de mechanismen die onscherpte en absorptie beïnvloeden in het reconstructieproces [6].
De generieke structuur van zo’n model zal aan bod komen in sectie 3.2. Er zal blijken dat het
niet praktisch is om het bekomen stelsel van vergelijkingen op te lossen, waardoor algoritmes
zijn ontworpen die in plaats daarvan iteratief de afbeelding moeten reconstrueren. Tijdens elke
iteratie tracht men dichter bij de afbeelding te komen die de werkelijke inverse van het probleem
zou opleveren. De beeldkwaliteit is beter dan bij de analytische methodes, maar het iteratieve
proces vraagt wel meer rekentijd.
Een iteratief reconstructiealgoritme zal een numerieke ‘kost’, bepaald uit de detector- of beelddata, trachten te minimaliseren. De link met kalibreren is dus vlug gemaakt: daar willen we aan
de hand van een optimalisatiealgoritme een bepaalde kost in functie van de geometrische parameters minimaliseren. Gezien het belang van iteratieve reconstructiealgoritmes in dit onderzoek,
wijden we Hoofdstuk 3 volledig aan dat onderwerp.
2.4
Micro-CT
De preklinische beeldvorming maakt vrijwel uitsluitend gebruik van micro-CT. De gebruikte
scanners zijn kleiner en de beelden vertonen een hogere spatiale resolutie door een sterkere
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
8
Figuur 2.5: De TriFoil Triumph-II X-O CT-scanner gebruikt in dit onderzoek.
vergrotingsfactor. De onderzoeksgroep MEDISIP (Medical Imaging and Signal Processing) van
de vakgroep ELIS (Elektronica en Informatiesystemen) aan de UGent kan beroep doen op de
diensten van Infinity, een beeldlab voor kleine proefdieren.
De ultieme tests van dit onderzoek zullen uitgevoerd worden op data bekomen van de TriFoil
Triumph-II X-O CT-scanner, afgebeeld op figuur 2.5. Secties 10.1 en 10.2 lichten in meer detail
het gebruikte object toe, samen met andere eigenschappen van deze scanner.
2.5
2.5.1
Kalibreren
Probleemstelling
De iteratieve reconstructiealgoritmes gaan uit van een geparameteriseerd model van de CTscanner (de precieze parameters komen aan bod in tabellen 4.2 en 9.2, respectievelijk voor reconstructie in twee dimensies en in drie dimensies). De kwaliteit van het tomogram is afhankelijk
van de kwaliteit van het model en de parameters waarmee de iteratieve algoritmes het beeld
reconstrueren. De projecties die de CT-scanner oplevert, voldoen niet altijd aan de geometrische veronderstellingen van het model. Dit resulteert in wazige, onscherpe tomogrammen met
mogelijks een afname in spatiale resolutie.
De werkelijke geometrie kan afwijken van de veronderstelde door bijvoorbeeld mechanische effecten. Zo speelt de invloed van de zwaartekracht op de stralingsbron een rol, alsook atmosferische
invloeden zoals temperatuur of druk. Er is nood aan het kalibreren van CT-scanners. Indien het
geparameteriseerde model de CT-scanner afdoende beschrijft, en we kunnen parameters vinden
die met dit model een goede reconstructie opleveren, dan kennen we met grote zekerheid de
werkelijke geometrische parameters van de CT-scanner.
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
2.5.2
9
Literatuur
Het kalibreren van CT-scanners kan ruwweg op drie manieren. De eerste manier bestaat erin een
object te scannen waarvan men de geometrische eigenschappen exact kent. Door het tomogram
te vergelijken met het ideale tomogram van dat object, kan men de geometrische parameters
bijsturen [7]. Men leidt de parameters af uit de soorten artefacten die optreden. De tweede
manier brengt markeringen aan op het te scannen object. Door de markeringen nauwkeurig te
volgen op de beeld- en projectiedata kan men ook een inschatting van de geometrische parameters
maken [8]. Deze methode komt voornamelijk voor in de context van elektronentomografie [9].
De derde methode laat het iteratief reconstructiealgoritme zelf de parameters schatten, zonder
voorkennis over het object of markeringen erop aangebracht. In de literatuur komen twee pistes voor om CT-scanners te kalibreren tijdens iteratieve beeldreconstructie: optimalisatie in de
beeldruimte [10] en optimalisatie in de projectieruimte [11]. Optimalisatie in de beeldruimte
betekent dat het tomogram wordt geëvalueerd op basis van zijn visuele eigenschappen om een
uitspraak te doen over de accuraatheid van de gebruikte geometrische parameters. Optimalisatie
in de projectieruimte betekent dat de projectiedata geëvalueerd wordt, door een vergelijking te
maken tussen de data bekomen van de CT-scanner en de projectiedata die het reconstructiemodel
oplevert. Beide pistes hebben gemeen dat men een objectieffunctie definieert die geminimaliseerd
of gemaximaliseerd moet worden. Het kaliberen van CT-scanners herleidt zich tot het oplossen
van een optimalisatieprobleem.
De meeste methodes om een objectieffunctie in de beeldruimte te minimaliseren of maximaliseren, bouwen verder op het werk van Kyriakou et al. [12]. Zij stellen voor om de entropie van
het bekomen tomogram te minimaliseren. Het minimaliseren van de entropie betekent het bevoordelen van hoog contrast. Dat een hoog contrast vanzelfsprekend een optimale reconstructie
zou betekenen, wordt in vraag gesteld door Kingston et al. [13]. Zij argumenteren dat entropie
spatiale informatie van het beeld negeert en geen scherpte in de hand werkt, maar eerder het
kunnen onderscheiden van segmenten in het tomogram. Zelf stellen zij scherpte als te maximaliseren objectieffunctie in de beeldruimte voor, gedefinieerd als de norm van de beeldgradiënt. Een
belangrijke kanttekening is dat scherpte strikt genomen niet schaalinvariant is: scherpte wijzigt
bij een variërende vergrotingsfactor. De vergrotingsfactor constant houden, is een aanname die
men in de praktijk niet altijd kan maken. We zullen ons moeten afvragen of parameters geoptimaliseerd kunnen worden bij een constante vergrotingsfactor die zelf niet noodzakelijk optimaal
is.
In de projectieruimte is het minimaliseren van de projectiefout gangbaar. Dit is een maat voor het
verschil tussen de scannerdata en de projecties die het scannermodel en reconstructiealgoritme
opleveren. Toonaangevend zijn onder meer Wicklein et al. [14] en Wein et al. [15]. Het is echter
onduidelijk of een optimaal tomogram altijd bekomen wordt bij de minimale projectiefout.
Debbeler et al. [16] geven de voorkeur aan een optimalisatie in de projectieruimte omwille van de
snelheid: het berekenen van het tomogram kan achterwege gelaten worden. De objectieffunctie
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
10
die zij voorstellen brengt de aanwezigheid van overtollig opgemeten x-stralen in rekening. Hiermee kunnen zij inderdaad tot nauwkeurige kalibratie komen, maar er stellen zich twee problemen.
Ten eerste bezitten niet alle projectiedata overtollig opgemeten x-stralen. Hun voorgestelde oplossing brengt dan weer niet alle geometrische parameters in rekening, waardoor het kalibreren
enkel lukt onder geïdealiseerde voorwaarden met zekerheid over de parameters die men niet kalibreert. Ten tweede blijken de objectieffuncties lokale extrema te bezitten, waardoor vraagtekens
bij de betrouwbaarheid geplaatst moeten worden.
Ook Patel et al. [17] werken in de projectieruimte, en behalen veelbelovende resultaten met betrekking tot het inschatten van rotatie- en translatieparameters van het object. Ook bewegingen
van het object tijdens het scanproces kunnen zij corrigeren. Ondanks het succes beperkt dit de
bruikbaarheid van hun voorgestelde methode aanzienlijk, gezien andere parameters (bijvoorbeeld
met betrekking tot de detector en stralingsbron) buiten beschouwing gelaten worden. Zij worden
gekend verondersteld.
Bronnikov [18][19] garandeert met een nieuwe, algemene theorie en bijhorend algoritme stabiele
convergentie naar de optimale parameters van het systeem. Deze methode brengt alle parameters
in rekening, maar maakt gebruik van de systeemmatrix van het CT-model die vaak onbekend of
analytisch moeilijk te hanteren is.
Stopp et al. [20] werken een kalibratieproces uit door markeringen op het object aan te brengen.
In de praktijk is het niet altijd mogelijk dit nauwkeurig te doen. Fouten op de behaalde resultaten
zijn inderdaad voornamelijk te wijten aan onnauwkeurig aangebrachte markeringen. Verder komt
het werk voornamelijk tegemoet aan CT-systemen waarbij de stralingsbron geen planair traject
aflegt.
Als optimalisatiealgoritme duikt vaak het Nelder-Mead algoritme op [11] [12], al lijkt een genetisch algoritme uit recent onderzoek van Sawall et al. [3] op het eerste zicht veelbelovender.
Opvallend is dat zelden parameters apart geoptimaliseerd worden, onafhankelijk van welke methode men gebruikt.
2.5.3
Doelstelling
Uit dit beknopt literatuuroverzicht is gebleken dat er nog geen algemene methode bestaat om alle
geometrische parameters te optimaliseren — de methode van Bronnikov [18] veronderstelt kennis
van de systeemmatrix. Scherpte als objectieffunctie is nog onvoldoende gekarakteriseerd en werd
niet vergeleken met resultaten behaald uit het gebruik van de projectiefout. Een combinatie van
objectieffuncties zal sowieso nodig zijn, willen we alle parameters optimaliseren.
Concreet trachten we in dit onderzoek na te gaan of er een algemeen toepasbare methode bestaat, die zo weinig mogelijk veronderstellingen over de geometrische parameters en het object
maakt, zonder gebruik te maken van de systeemmatrix. We willen alle geometrische parameters
HOOFDSTUK 2. COMPUTERTOMOGRAFIE
11
van het model kunnen kalibreren en nagaan of we bepaalde parameters onafhankelijk van de andere kunnen optimaliseren. Er moet nagegaan worden welke objectieffuncties het best geschikt
zijn voor welke parameters, door de invloed van de parameters op beeld- en projectiedata na te
trekken. Ook de tijd speelt een rol: het moet duidelijk worden na hoeveel iteraties van het reconstructiealgoritme we de objectieffuncties betrouwbaar kunnen evalueren, iets wat in de literatuur
vaak in het midden gelaten wordt. Het is verder interessant te weten welk optimalisatiealgoritme
aangewezen is om specifiek dit CT-probleem aan te pakken.
2.5.4
Verloop
Het verloop van dit onderzoek is ruwweg als volgt.
• Hoofdstuk 3 verklaart de grondslagen van iteratieve beeldreconstructie. Deze voorkennis
is nodig, omdat het tijdens deze reconstructie is dat het kalibreren zal plaatsvinden.
• Hoofdstuk 4 legt de basis voor een voorstudie van kalibreren in twee dimensies.
• Hoofdstuk 5 bespreekt en vergelijkt mogelijke optimalisatiealgoritmes om het kalibreren
uit te voeren. Ook de objectieffuncties worden er gedefiniëerd.
• Hoofdstukken 6 en 7 behandelen respectievelijk de projectiefout en scherpte. Prestaties worden geëvalueerd en eerste resultaten voorgelegd. Hoofdstuk 8 gaat na of het
nut heeft om deze objectieffuncties te combineren. Deze tweedimensionale voorstudie moet
voldoende inzicht in het kalibreren geleverd hebben om vlot over te gaan naar een implementatie op de GPU, om ‘echte’ data van kleine proefdieren in drie dimensies te reconstrueren.
• Hoofdstuk 9 beschrijft welke aanpassingen nodig zijn voor een uitbreiding naar drie dimensies, met een implementatie op de GPU. Hoofdstuk 10 behandelt in detail het kalibreren van de CT-scanner.
• Hoofdstuk 11 geeft een conclusie bij het volledige onderzoek.
Implementatiedetails, bijvoorbeeld met betrekking tot specifieke C++ syntax, worden zoveel
mogelijk vermeden in de uiteenzetting die volgt. Toch maakte het programmeren een aanzienlijk
deel van dit onderzoek uit, zodat de code wel geannoteerd op een digitale drager aan deze thesis
is toegevoegd. Bijlage A beschrijft de inhoud van deze cd-rom.
Hoofdstuk 3
Iteratieve beeldreconstructie
In dit hoofdstuk lichten we de principes van iteratieve beeldreconstructie toe, omdat het tijdens
dit proces is dat we de beeld- en/of detectordata zullen evalueren om de geometrische parameters
van de CT-scanner te optimaliseren.
3.1
Beeld- en projectieruimte
Met elke tweedimensionale pixel (picture element) van het tomogram komen in het object één
of meerdere driedimensionale voxels (volume elements) overeen. De dikte van de doorsnede
bepaalt de derde dimensie van de voxel. Deze dikte is bij CT dun (tussen 1 mm en 10 mm)
en uniform over de doorsnede. Het beeld wordt weergegeven in grijswaarden, die zoals getoond
op figuur 2.4 uit de Hounsfieldeenheden worden afgeleid. De Hounsfieldeenheden worden op
hun beurt afgeleid uit de gemiddelde absorptie-eigenschappen door de detectoren geregistreerd,
zoals vergelijkingen 2.1 en 2.2 reeds toonden. Het tomogram wordt middels een algoritme uit de
detectordata gereconstrueerd. Er bestaat met andere woorden een afbeelding tussen de projectieen beeldruimte.
3.2
Lineair model
Het te scannen object stellen we wiskundig voor met de objectdistributie f . Zeer algemeen
wordt een object ook wel ‘fantoom’ genoemd, wegens het nog onbekend zijn van de inwendige
structuren die we net door het object te scannen, willen visualiseren. Het resultaat van de
scan zijn de projectiewaarden g die door de detector tijdens de acquisitie gemeten werden. De
systeemmatrix H geeft het verband tussen f en g [6]:
g = Hf
12
(3.1)
HOOFDSTUK 3. ITERATIEVE BEELDRECONSTRUCTIE
13
Figuur 3.1: Verband tussen f en g door systeemmatrix H.
f is een kolomvector met dimensie N ×1, g heeft dimensie P ×1 en bijgevolg is H een P ×N matrix.
N is het aantal voxels van de te reconstrueren afbeelding (f kan hier, zonder de algemeenheid te
schaden, ook tweedimensionaal zijn, bijvoorbeeld in het geval van één centrale doorsnede). P is
het aantal detectorcellen. Een enkele voxel (pixel) zullen we verder fj , j ∈ {1, . . . , N } noemen.
Naar analogie wordt een enkele detectorcel omschreven als gi , i ∈ {1, . . . , P }.
Merk op dat zowel f als g hier discreet zijn. Voor het object betekent dit het invoeren van een
basisfunctie φj (x), met x een plaats in de continue beeldruimte, om een individuele pixelwaarde
te bekomen [6]:
Z
fj =
f (x)φj (x) dx
(3.2)
De meest gebruikte basisfuncties φj (x) zijn constant in een welbepaalde kubus (voxel) uit een
verzameling van kubussen waarin het object is onderverdeeld.
Voorts is de projectiedata g strikt genomen een toevalsvariabele door de emissie van fotonen
waarvan men weet dat ze een Poissondistributie volgen. In dat opzicht moeten we vergelijking
3.1 als volgt herschrijven:
E[g] = Hf
(3.3)
Men heeft het hier over de verwachtingswaarde of gemiddelde waarde van de gemeten data g.
De systeemmatrix H is een zuiver wiskundig concept om het verband tussen f en g gestalte te
geven. Het gebruik van H is in de praktijk moeilijk, wegens vaak te groot om uit te rekenen.
Ook het bepalen van H−1 zou dan erg rekenintensief zijn. Toch is net deze inverse nodig om uit
vergelijking 3.1 f te bepalen, het doel van beeldreconstructie:
f = H−1 g
(3.4)
Het (praktisch) niet kunnen inverteren van H geeft aanleiding tot iteratieve reconstructietechnieken.
HOOFDSTUK 3. ITERATIEVE BEELDRECONSTRUCTIE
3.3
14
Iteratieve reconstructie
Een iteratieve reconstructie bestaat steeds uit twee componenten. De eerste component is het
criterium dat zal bepalen welke afbeelding ˆf de beste schatting is van de echte afbeelding f , die
het inverteren in vergelijking 3.4 zou opleveren. De tweede component is het algoritme dat zo’n
schatting ˆf iteratief tot stand zal brengen. Bij elke iteratie k van dit algoritme zal de schatting
beter worden en zal de beeldkwaliteit verbeteren:
lim ˆf (k) = f
k→+∞
3.3.1
(3.5)
Criteria
De criteria kunnen onderverdeeld worden in twee categorieën: een klasse van criteria die de
probabiliteit van f en/of g mee in rekening brengt en een klasse die geen veronderstellingen over
het object f doet. Er volgen nu twee voorbeelden die voor het verder verloop van dit onderzoek
belangrijk zullen zijn.
Criterium van de grootste aannemelijkheid
Het Poissonkarakter van g dat aanleiding gaf tot vergelijking 3.3, leidt formeel tot deze probabiliteitswet voor g:
P
Y
g¯igi e−¯gi
(3.6)
p(g; f ) =
gi !
i=1
g¯i is element i van E[g] = Hf :
g¯i =
N
X
hij fj
(3.7)
j=1
Hierin stelt hij het element op positie (i, j) van de systeemmatrix H voor. Men hanteert vervolgens de meest aannemelijke schatter om ˆf uit de aannemelijkheidsfunctie 3.6 te bepalen [6]:
ˆf = arg max p(g; f )
f
(3.8)
Men kiest met andere woorden een afbeelding waarvoor de aannemelijkheidsfunctie p(g, f ) een zo
groot mogelijke waarde oplevert. Deze techniek wordt onder meer gebruikt in MLTR (maximum
likelihood transmission reconstruction).
Criterium van de kleinste kwadraten
Een ander criterium bestaat erin ˆf te kiezen zodat de projectie ervan de kleinste Euclidische
afstand tot de projectiedata g vertoont. Deze schatter wordt formeel als volgt uitgedrukt, volgens
HOOFDSTUK 3. ITERATIEVE BEELDRECONSTRUCTIE
15
Figuur 3.2: De generieke structuur van een iteratief reconstructiealgoritme.
de methode van de kleinste kwadraten [6]:
ˆf = arg min kg − Hf k2 = arg min
2
f
f
2
P N
X
X
gi −
hij fj
i=1
(3.9)
j=1
Vergelijking 3.9 ligt onder meer aan de basis van SIRT (simultaneous iterative reconstruction
technique), in meer detail besproken in sectie 3.4.
3.3.2
Algoritmes
De generieke structuur van een iteratief reconstructiealgoritme is weergegeven op figuur 3.2. Het
algoritme maakt gebruik van een computermodel van de CT-scanner dat virtuele projectieoperators gebruikt in plaats van de systeemmatrix H. Een schatting ˆf (k) wordt geprojecteerd op
ˆ (k) = Hˆf (k) bekomen wordt. Deze projectie wordt
een virtuele detector zodat een projectie g
vergeleken met de gemeten projectiedata g. Deze vergelijking houdt het bepalen van een fout
in. Het is deze fout die (wederom virtueel) terugwaarts geprojecteerd wordt om opnieuw een
afbeelding te bekomen. De oude schatting ˆf (k) wordt met deze nieuwe afbeelding aangepast tot
schatting ˆf (k+1) . De reconstructie stopt na een vooraf bepaald aantal iteraties of wanneer de
afbeelding onvoldoende verbetert per toenemende iteratie k. Een mogelijke metriek hiervoor is
het bepalen van de L2 -norm van het verschil tussen ˆf (k+1) en ˆf (k) . Het doel van een iteratief
reconstructiealgoritme is ervoor zorgen dat ˆf (k) (zo snel mogelijk) naar f convergeert, gelet op
vergelijking 3.5.
3.4
SIRT
De simultaneous iterative reconstruction technique, kortweg SIRT, maakt gebruik van het criterium van de kleinste kwadraten uit vergelijking 3.9. Algoritme 3.1 verduidelijkt de werking van
HOOFDSTUK 3. ITERATIEVE BEELDRECONSTRUCTIE
16
Algoritme 3.1: SIRT
1
2
foreach iteratie k do
foreach projectiehoek α do
3
(k)
gα,err ← gα −VoorwaartseProjectie(ˆf (k) ,α);
4
gα,err ← Normaliseer(gα,err );
5
ferr ← ferr +TerugwaartseProjectie(gα,err ,α);
(k)
(k)
(k)
(k)
(k)
6
end
7
8
ferr ← Normaliseer(ferr );
(k)
ˆf (k) =ˆf (k) + ferr
;
9
k ← k + 1;
10
(k)
(k)
end
SIRT [21]. De voorwaartse en terugwaartse projectie worden bij gebrek aan H uitgevoerd door
geprogrammeerde projectieoperators.
In sectie 5.2.2 zullen we de projectiefout als te minimaliseren objectieffunctie definiëren met oog
op het kalibreren van CT-scanners. Deze metriek is direct verwant met de fout die in SIRT
terugwaarts wordt geprojecteerd. Daarom beslissen we om in dit onderzoek steeds SIRT te gebruiken, mede door de beschikbaarheid van een robuuste SIRT-implementatie op zowel de CPU
als GPU. Bovendien zal in sectie 9.5 blijken dat SIRT zich gemakkelijker leent tot een parallelle
implementatie op de GPU dan het intrinsiek snellere OS-SART (ordered-subset simultaneous algebraic reconstruction technique). MLTR toepassen op de scannerdata bleek onmogelijk, wegens
te grote dimensies van de detector en dus beperkingen op het geheugen van de GPU. MLTR zal
hier dus buiten beschouwing gelaten worden. Bovendien zou het aanpassingen vergen om ook in
MLTR met de projectiefout uit sectie 5.2.2 te werken, gezien MLTR het criterium van de grootste aannemelijkheid gebruikt en niet het criterium van de kleinste kwadraten dat rechtstreeks
verband houdt met de projectiefout uit sectie 5.2.2.
Hoofdstuk 4
Reconstructie in twee dimensies
We leggen nu de basis van tweedimensionale reconstructies. We gebruiken zoals gezegd SIRT,
waarbij de voorwaartse en terugwaartse projectie afstandsgebaseerd (distance-driven) zijn [22].
Om de complexiteit binnen de perken te houden, gaat het om de reconstructie van één centrale doorsnede van een object. Dit reduceert het probleem tijdelijk tot een tweedimensionale
reconstructie die binnen redelijke tijd op een CPU uitgevoerd kan worden. Het maakt het definiëren van objectieffuncties en het testen van optimalisatiealgoritmes aanvankelijk vlotter dan een
rechtstreekse implementatie op een GPU, dat meteen meer programmeerwerk zou vergen. We
gaan pas over naar een implementatie op de GPU nadat in deze voorstudie de meest efficiënte
technieken om te kalibreren op punt werden gesteld.
4.1
Geometrische parameters
Figuur 4.1 toont het gebruikte coördinatensysteem van het model voor reconstructie in twee
dimensies. Merk op dat de detector zijn eigen (u, v)-coördinatensysteem heeft, terwijl de detector
wel samen met de stralingsbron in het getoonde (x, y)-coördinatensysteem van het object roteert.
Figuur 4.1: Coördinatensysteem voor reconstructie in twee dimensies.
17
HOOFDSTUK 4. RECONSTRUCTIE IN TWEE DIMENSIES
18
Beschrijving
Eenheid
Totale rotatie
Aantal projecties
Dimensie object in x-richting
Dimensie object in y-richting
Dimensie objectpixel (isotroop)
Dimensie detector in u-richting
Dimensie detectorpixel
◦
—
pixels
pixels
mm
pixels
mm
Tabel 4.1: Definitie van de vaste geometrische parameters in twee dimensies.
Symbool
Beschrijving
Eenheid
∆u
∆v
θ
SD
SO
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (tilt)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
mm
mm
◦
mm
mm
Tabel 4.2: Definitie van de te optimaliseren geometrische parameters in twee dimensies.
Tabel 4.1 geeft een overzicht van de geometrische parameters die we in twee dimensies onderscheiden. Het zijn de ‘vaste’ parameters die we niet zullen kalibreren. Tabel 4.2 bevat de parameters
die we zullen optimaliseren naargelang het gedrag van de objectieffuncties. De uitwijkingen van
de detector in het (u, v)-coördinatensysteem worden gevisualiseerd op figuur 4.2, waar de pijlen
steeds de positieve zin aanduiden. We kiezen ervoor om de afstand tussen object en detector
enkel te beïnvloeden door de uitwijking ∆v van de detector aan te passen. SD blijft dan constant.
Dit brengt het aantal te optimaliseren parameters in deze voorstudie op vier.
We behandelen nu kort hoe we in de implementaties ForwardProject2D_geometry.hh en
SIRT_geometry.hh de coördinaten van de verschoven detector in het (x, y)-coördinatensysteem
van het object bepalen. Deze coördinatentransformatie hangt af van ∆u, ∆v, θ, SD en SO, maar
ook van de projectiehoek α, waarvan de zin en richting op figuur 4.1 getoond werden.
Veronderstel dat we de detector op de x-as leggen, met zijn middelpunt op de oorsprong. Noem
xgi dan de x-coördinaat van detectorcel gi . Tijdens de rotatie van stralingsbron en detector rond
het object, heeft deze detectorcel in het (x, y)-coördinatensysteem de coördinaten (x0gi , yg0 i ). Dit
coördinatenpaar bepalen we als volgt:
x0gi = cos(α) xgi cos(θ) + ∆u + sin(α) − xgi sin(θ) − (SD − SO) + ∆v
yg0 i = cos(α) − xgi sin(θ) − (SD − SO) + ∆v − sin(α) xgi cos(θ) + ∆u
(4.1)
(4.2)
Merk op dat we de parameters tijdens elke iteratie k gelijk houden voor elke projectiehoek
α. In principe geldt deze aanname enkel in het geval van een perfect planair, circulair traject
HOOFDSTUK 4. RECONSTRUCTIE IN TWEE DIMENSIES
(a) ∆u
(b) ∆v
19
(c) θ
Figuur 4.2: Uitwijkingen van de detector in het (u, v)-coördinatensysteem.
van de stralingsbron rond het object. Toch maken we deze veronderstelling omdat anders de
complexiteit van het optimalisatiealgoritme onhandelbaar groot zou worden.
4.2
4.2.1
Object
Definitie
Het object met distributie f , getoond op figuur 4.3, is in deze voorstudie een zelfontworpen
fantoom van 128 × 128 pixels met enkele eenvoudige geometrische patronen: cirkel, driehoek
en rechthoeken. De scherpe randen moeten ons in staat stellen gemakkelijk de scherpte van
reconstructies te beoordelen, die berekend zal worden aan de hand van de randen. De lijnen met
variërende dikte en grijswaarden staan in verschillende hoeken ten opzichte van elkaar opdat bepaalde artefacten zouden optreden bij het slecht inschatten van de geometrische parameters. De
exacte kennis van f en de geparameteriseerde projectieoperator H stelt ons in staat de behaalde
resultaten ondubbelzinnig te evalueren. Het is vanzelfsprekend — en net het uitgangspunt van
dit onderzoek — dat f en H in werkelijkheid onbekend zijn. Deze voorstudie zal dus voorkennis aanwenden om de objectieffuncties en optimalisatiealgoritmes zo robuust mogelijk te maken,
zodat ze later toepasbaar zijn op echte data met onbekende f en H.
Figuur 4.3: Tweedimensionaal fantoom.
HOOFDSTUK 4. RECONSTRUCTIE IN TWEE DIMENSIES
Symbool
∆u
∆v
θ
SD
SO
20
Beschrijving
Waarde
Eenheid
Totale rotatie
Aantal projecties
Dimensie object in x-richting
Dimensie object in y-richting
Dimensie objectpixel (isotroop)
Dimensie detector in u-richting
Dimensie detectorpixel
360
360
128
128
0,15625
1184
0,1
◦
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (tilt)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
-3,0
6,0
-9,0
290,0
151,0
—
pixels
pixels
mm
pixels
mm
mm
mm
◦
mm
mm
Tabel 4.3: Geometrische parameters bij voorwaartse projectie van het fantoom.
Symbool
Beschrijving
Waarde
Eenheid
∆u
∆v
θ
SD
SO
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (tilt)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
0,0
290,0
145,0
mm
mm
◦
mm
mm
Tabel 4.4: Initiële geometrische parameters bij reconstructie.
4.2.2
Parameterkeuze
Tabel 4.3 vat de parameters samen van het fantoom en van de geometrie van de (fictieve) CTscanner. Het zijn deze parameters die we zullen gebruiken bij de voorwaartse projectie van
het fantoom, om de detectordata g te simuleren. Tabel 4.4 beschrijft voor de te optimaliseren
parameters de waarden waarvan men denkt dat ze de werkelijke geometrie beschrijven. Het
ultieme doel van het kalibrereren is om vanuit de waarden in tabel 4.4 de waarden uit tabel 4.3
te bereiken. In werkelijkheid zal een betere initiële schatting mogelijk zijn, maar het gebruik van
een dergelijke slechte schatting laat ons toe te evalueren hoe goed de optimalisatiealgoritmes de
parameterruimte exploreren.
4.2.3
Reconstructies
Figuur 4.4 toont voor verschillende iteraties k van SIRT het bekomen tomogram. Het betreft
de perfecte reconstructie met parameters uit tabel 4.3. Een kwantitatief beter resultaat is met
andere woorden niet haalbaar. Kwalitatief kunnen restauratietechnieken de reconstructie even-
HOOFDSTUK 4. RECONSTRUCTIE IN TWEE DIMENSIES
(a) k = 1
(b) k = 10
(c) k = 100
21
(d) k = 500
Figuur 4.4: Reconstructie met correcte parameters (tabel 4.3).
(a) k = 1
(b) k = 10
(c) k = 100
(d) k = 500
Figuur 4.5: Reconstructie met initiële parameters (tabel 4.4).
tueel nog verbeteren, maar zij worden in dit onderzoek buiten beschouwing gelaten. Figuur 4.5
toont voor dezelfde iteraties k de afbeelding die bekomen wordt door te reconstrueren met de
parameters uit tabel 4.4.
De impact van de foutief ingeschatte parameters is duidelijk groot. Het valt verder op dat het
onderscheid van fijn detail, zoals de zichtbaarheid van de randen en de ruimte tussen de balkjes,
sterk verbetert naarmate het aantal iteraties k toeneemt. Dit geldt zowel voor de perfecte reconstructie als voor de foute schatting. Dat zelfs de perfecte reconstructie na 500 iteraties nog steeds
kleine randartefacten vertoont, is deels te wijten aan de aard van het reconstructiealgoritme, in
dit geval SIRT.
Figuren 4.6 tot en met 4.9 tonen het effect van telkens één foute geometrische parameter. De
uitwijking ∆u van de detector heeft duidelijk de grootste invloed. Een foute inschatting ervan
levert dubbele randen en een lager intrinsiek contrast op. ∆v en SO beïnvloeden in eerste plaats
de vergroting, maar ook de scherpte: de foute inschatting van SO toont wazige randen van de
fijne balkjes rechts. De detectorhoek θ tast de scherpte van de randen aan, naast het introduceren
van dubbele randen, zij het minder uitgesproken dan in het geval van ∆u. Extreme waarden van
θ beïnvloeden ook de vergroting. De vier voorbeelden tonen dat de horizontale artefacten in de
linkerbalk er telkens meer uitgesproken zijn dan in de perfecte reconstructie op figuur 4.4.
HOOFDSTUK 4. RECONSTRUCTIE IN TWEE DIMENSIES
(a) k = 1
(b) k = 10
(c) k = 100
22
(d) k = 500
Figuur 4.6: ∆u = 0, 0 mm in de aanwezigheid van correcte parameters.
(a) k = 1
(b) k = 10
(c) k = 100
(d) k = 500
Figuur 4.7: ∆v = 0, 0 mm in de aanwezigheid van correcte parameters.
(a) k = 1
(b) k = 10
(c) k = 100
(d) k = 500
Figuur 4.8: θ = 0, 0◦ in de aanwezigheid van correcte parameters.
(a) k = 1
(b) k = 10
(c) k = 100
(d) k = 500
Figuur 4.9: SO =145,0 mm in de aanwezigheid van correcte parameters.
Hoofdstuk 5
Optimalisatietechnieken
Dit hoofdstuk geeft een inleiding tot algoritmische optimalisatietechnieken. We trachten het
kalibreren van CT-scanners volledig te abstraheren tot een dergelijk optimalisatieprobleem, met
definitie en afweging van geschikte objectieffuncties en algoritmes.
5.1
5.1.1
Principe
Parameters
In tabel 4.2 werden reeds de te optimaliseren parameters beschreven. Naar analogie zal tabel 9.2
de parameters geven in het geval van driedimensionale reconstructies. Met Ω = {p1 , p2 , . . . , pn } ∈
Rn duiden we algemeen een verzameling van n parameters aan. De parameters bepalen de kwaliteit van een reconstructie, omdat reconstructiealgoritmes werken met geparameteriseerde projectieoperators H en H−1 . Wanneer we zeggen dat we CT-scanners willen kalibreren, bedoelen we
dat we tijdens iteratieve reconstructie de verzameling Ω willen bepalen om een zo goed mogelijke
reconstructie te bekomen. Gedurende k iteraties zullen we met dezelfde, door een optimalisatiealgoritme bepaalde parameterverzameling Ω werken. Een kwalitatief goede reconstructie is een
scherpe afbeelding met grote spatiale resolutie die dus voldoende detail vertoont.
Een goede reconstructie kan visueel waargenomen worden, maar om de parameters automatisch
te schatten, is er uiteraard een manier nodig om numeriek te bepalen wat een ‘goede reconstructie’
is. Er is nood aan de definitie van een zogenaamde objectieffunctie.
5.1.2
Objectieffuncties
Een objectieffunctie c : Ω ∈ Rn → R kent aan een parameterverzameling Ω een waarde toe
die men wil maximaliseren of minimaliseren, naargelang de betekenis van c. In het geval van
minimalisatie heeft men het ook wel eens over een ‘kostfunctie’ in plaats van een objectieffunctie.
23
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
(a) pmin ≤ p ≤ pmax
24
(b) pmin ≤ p ≤ pmax , a(p) ≥ 0
Figuur 5.1: Geïllustreerde terminologie in het vereenvoudigde geval Ω = {p}.
De parameters in Ω worden de decisievariabelen pi genoemd, met i ∈ {1, . . . , n}. Vergelijkingen
5.1 tot en met 5.3 zijn bijkomende, beperkende voorwaarden die men typisch oplegt in de context
van optimalisatieproblemen.
ai (Ω) = 0, i ∈ {1, . . . , m1 }
(5.1)
ai (Ω) ≤ 0, i ∈ {m1 + 1, . . . , m1 + m2 }
(5.2)
ai (Ω) ≥ 0, i ∈ {m1 + m2 + 1, . . . , m1 + m2 + m3 }
(5.3)
m = m1 + m2 + m3 is het aantal beperkende voorwaarden, met ai : Ω ∈ Rn → R gegeven functies
die ook een bepaalde kost definiëren, zonder dat die geminimaliseerd of gemaximaliseerd moet
worden. Bijkomend kan men afdwingen dat de parameters pi zich in een bepaald bereik moeten
bevinden. Het aantal geldige oplossingen in de oplossingenruimte wordt zo beperkt. Figuur 5.1
verduidelijkt het optimalisatieprincipe en de gebruikte terminologie.
Wanneer c en ai lineair zijn in de parameters van Ω, wordt het geschetste probleem opgelost
door zogenaamd lineair programmeren. In de context van CT is dit onmogelijk: H is praktisch
vaak niet te berekenen, wat ook het gebruik van H−1 bemoeilijkt. Objectieffuncties zullen
sowieso gedefinieerd moeten worden op de geschatte detectordata (minstens H benodigd) of op
de gereconstrueerde afbeeldingen zelf (H en H−1 benodigd). Sectie 5.2 zal specifiek voor het
CT-probleem objectieffuncties definiëren.
5.1.3
Optimalisatiealgoritmes
Met een optimalisatiealgoritme tracht men het globaal optimum stapsgewijs (iteratief) te bereiken. Een aantal algoritmes slagen erin te convergeren naar een dergelijke oplossing, maar dit
hangt af van de aard van het probleem. Door de onbruikbaarheid van de systeemmatrix H zijn
we beperkt tot heuristische algoritmes, waarbij het bereiken van het globaal optimum en convergentie nooit gegarandeerd zijn. We onderscheiden verschillende types heuristische algoritmes.
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
25
Generiek versus probleemafhankelijk
Generieke heuristieken zijn zonder veel aanpassingen inzetbaar voor een breed toepassingsgebied
van optimalisatieproblemen, onafhankelijk van de aard van het probleem. Dit in tegenstelling tot
probleemafhankelijke heuristieken die volledig ontworpen worden voor een specifiek optimalisatieprobleem. De intelligentie van dergelijke algoritmes beroept zich vaak op een soort historische
voorkennis van eerdere oplossingen van het probleem.
In dit onderzoek geven we de voorkeur aan generieke heuristieken, omdat hun implementaties
en prestaties over het algemeen goed gedocumenteerd zijn in de literatuur. Uiteraard zullen
specifieke aanpassingen nodig zijn naargelang de aard van de CT-context en de problemen die
zich erin stellen.
Zoekalgoritme versus constructief algoritme
Zoekalgoritmes starten vanuit een bepaalde oplossing die willekeurig is of reeds voldoende goed
bevonden werd. Het algoritme suggereert tal van andere oplossingen die uit vorige oplossingen ontstaan, door bijvoorbeeld uit goed presterende oplossingen parameters te combineren of
bepaalde parameters in een weloverwogen richting te duwen. Er ontstaat een traject door de
oplossingenruimte. Constructieve algoritmes daarentegen nemen een beslissing op basis van alle
parameters van de huidige oplossing en maken in grote mate gebruik van voorkennis over hun
gedrag en wederzijdse beïnvloeding.
Een constructief algoritme zou veel inzicht vragen in hoe de geometrische parameters van de CTscanner zich tot elkaar verhouden: hoe beïnvloeden ze elkaar? Welke parameters kunnen elkaar
opheffen? Het zou een sterk wiskundige studie vergen, die we omzeilen door aanvankelijk een
zoekalgoritme te gebruiken en later via reverse engineering alsnog empirisch verbanden tussen
de parameters vast te stellen. Dit komt uitvoerig in Hoofdstuk 6 en Hoofdstuk 7 aan bod.
Stochastisch versus deterministisch algoritme
Een laatste onderscheid in heuristische optimalisatiealgoritmes is het al dan niet gebruiken van
toevalsvariabelen die de waarden van de te optimaliseren parameters beïnvloeden. Deterministische algoritmen gebruiken vaste regels om nieuwe oplossingen te genereren, wederom gebaseerd
op voorkennis over het gedrag van de parameters. Een deterministisch algoritme zal daarom
altijd hetzelfde optimum bereiken wanneer het start vanuit dezelfde initiële oplossing.
Het werken met stochastische algoritmen zal ons in dit onderzoek toelaten om te kijken of het
optimalisatiealgoritme ondanks (dankzij?) de stochastische aard naar dezelfde oplossing streeft.
Bovendien laat het het zoeken van oplossingen toe die anders — afhankelijk van de kwaliteit van
de voorkennis om deterministische beslissingen te nemen — misschien onverkend zouden blijven.
Het is een manier om een lokaal van een globaal optimum te onderscheiden.
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
5.2
26
Objectieffuncties voor het CT-probleem
We kunnen objectieffuncties voor het CT-probleem op twee manieren definiëren. Ofwel berekent
men een numerieke waarde uit de geschatte detectordata. In het geval van één iteratie k vereist
dit enkel een voorwaartse projectie, voor k > 1 gaat het om k voorwaartse en k − 1 terugwaartse
projecties. Ofwel berekent men een numerieke waarde uit de beelddata. Dit vereist evenveel
voorwaartse als terugwaartse projecties.
Een ideale objectieffunctie heeft drie belangrijke eigenschappen [13]:
1. De objectieffunctie heeft een globaal minimum of maximum bij gebruik van de werkelijke
geometrische parameters.
2. Lokale minima of maxima ontbreken zoveel als mogelijk.
3. De objectieffunctie is betrouwbaar te evalueren bij onvolledige reconstructies (gering aantal
iteraties k, beperkt aantal projectiehoeken,. . . )
De laatste eigenschap is er vooral op gericht om het kalibreren binnen redelijke tijd uit te voeren.
5.2.1
Entropie
Het minimaliseren van de entropie E heeft als doel het contrast van de gereconstrueerde afbeelding te vergroten, in de veronderstelling dat het beste contrast optreedt bij de parameters
die zich het dichtste bij de werkelijke parameters bevinden. Entropie als objectieffunctie wordt
geëvalueerd op de beelddata. Formeel [12]:
E=−
B
X
h(b) · ln(h(b))
(5.4)
b=1
H(b)
h(b) = PB
b=1 H(b)
(5.5)
H is het histogram van de gereconstrueerde afbeelding. Een histogram verdeelt het grijswaardebereik van de gereconstrueerde afbeelding in B bins van gelijke grootte. Voor elke bin
b ∈ {1, . . . , B} geeft H(b) het aantal pixels (voxels) in de afbeelding binnen het grijswaardebereik
dat b voorstelt.
Entropie moet geïnterpreteerd worden als een maat voor ‘onzekerheid’. Een lage entropie betekent weinig onzekerheid: de kans om een pixel binnen een bepaald grijswaardenbereik aan te
treffen, is groot. Op het histogram kenmerkt zich deze situatie door een beperkt aantal pieken
bij bepaalde grijswaarden of bins. Een hoge entropie betekent dat er veel pixels zijn met verschillende grijswaarden. Men kan in dit geval op het histogram geen nauwe regio aanduiden waarin
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
(a) onscherp
27
(b) scherper
Figuur 5.2: Onscherpe randovergangen bevatten een groter aantal verschillende grijswaarden.
(a) onscherp
(b) scherper
Figuur 5.3: Histogrammen in overeenstemming met figuur 5.2.
de meeste pixels zich bevinden. Figuren 5.2 en 5.3 verduidelijken dit. De scherpe afbeelding
heeft een lagere entropie dan de onscherpe.
Het voordeel van het minimaliseren van entropie is de relatieve eenvoud van het proces. Een
histogram opstellen kan over het algemeen snel, zonder veel rekenwerk. Er zijn echter twee
belangrijke nadelen aan entropie als objectieffunctie verbonden [13]:
1. Entropie negeert spatiale informatie.
2. Entropie minimaliseren neigt naar het bevorderen van segmenten in het beeld in plaats van
algemene scherpte.
Deze beperkingen dwingen ons verder op zoek te gaan naar andere objectieffuncties, die in de
literatuur aanvankelijk minder voorkwamen dan entropie (zie ook sectie 2.5.2), maar tegemoetkomen aan voornoemde beperkingen.
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
(a) E (100) = 1, 00
28
(b) E (100) = 6, 67
(c) E (100) = 852, 02
Figuur 5.4: Genormaliseerde projectiefout voor verschillende reconstructies.
5.2.2
Projectiefout
(k)
De projectiefout EΩ is gedefinieerd op de detectordata en ziet er formeel als volgt uit [15]:
(k)
EΩ
=
P
X
(k)
|gi − (HˆfΩ )i |
(5.6)
i=1
ˆf (k) is de afbeelding die ontstaat na k iteraties van het iteratief reconstructiealgoritme met geΩ
bruik van een bepaalde parametervector Ω. Figuur 5.4 toont de projectiefout voor verschillende
reconstructies (de waarden werden genormaliseerd op basis van de waarde bij perfecte reconstructie).
Het gebruik van de projectiefout als objectieffunctie is gesteund op het feit dat (een klasse van)
iteratieve reconstructiealgoritmes deze projectiefout wil minimaliseren, volgens vergelijking 3.9.
Hoofdstuk 6 bespreekt het minimaliseren van de projectiefout in detail.
5.2.3
Scherpte
Net zoals entropie definiëren we scherpte op de beelddata. Scherpe afbeeldingen bevatten over
het algemeen de maximale hoeveelheid hoogfrequente informatie mogelijk. Slechte scherpstelling of het wazig maken van randen zullen deze hoogfrequente informatie typisch onderdrukken.
Scherpte als objectieffunctie slaagt erin hoogfrequente data kwantitatief uit te drukken en zou
over het algemeen een globaal maximum bezitten [13], hoewel net zoals bij de projectiefout zal
blijken dat dit maximum afhangt van het aantal iteraties k van het iteratief reconstructiealgoritme.
(k)
(k)
(k)
De scherpte SΩ van ˆfΩ is formeel gedefinieerd als de L2 -norm van de gradiënt van ˆfΩ [13]:
(k)
(k)
SΩ = k∇ˆfΩ k2 =
XX
x
y
(k)
|∇fˆΩ (x, y)|2
(5.7)
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
(a) reconstructie
(b) convolutie G
29
(c) convolutie GT
Figuur 5.5: Randdetectie door middel van het Sobelmasker G.
(a) reconstructie
(b) convolutie G
(c) convolutie GT
Figuur 5.6: Randdetectie door middel van het Sobelmasker G.
De gegeven formule geldt in twee dimensies. In geval van driedimensionale reconstructie kan de
scherpte per doorsnede gesommeerd worden om de scherpte van de volledige afbeelding uit te
drukken.
(k)
|∇fˆΩ |2 kan op verschillende manier bepaald worden, waaronder de afleidingseigenschap van de
eendimensionale, discrete Fouriertransformatie. Een eenvoudige benadering bestaat er echter in
het Sobelmasker G te gebruiken:
(k)
(k)
(k)
|∇fˆΩ |2 = (G ∗ ˆfΩ )2 + (GT ∗ ˆfΩ )2


−1 0 1


G = −2 0 2
−1 0 1
(5.8)
(5.9)
De operator ∗ is de convolutieoperator. Het gebruik van het Sobelmasker G betekent dat het
bepalen van de scherpte gebaseerd is op randdetectie. G detecteert horizontale randen, GT
verticale. Figuur 5.6 toont de invloed van het Sobelmasker op een gereconstrueerde afbeelding.
Figuur 5.7 toont de genormaliseerde scherpte voor verschillende reconstructies (normalisatie met
de scherpte van de perfecte reconstructie). Het Sobelmasker is ook in drie dimensies te gebruiken,
waardoor het sommeren van de scherpte per doorsnede niet meer nodig is.
Het grootste nadeel is dat de scherpte strikt genomen niet schaalinvariant is: de scherpte neemt
toe met de vergrotingsfactor M = SO/SD [13]. De veronderstelling dat M constant moet zijn, is
in de praktijk vaak niet mogelijk. Bovendien, mochten we toch met een centraal object werken,
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
30
dan kan dit object invariant zijn voor sommige parameters, waardoor de conclusies getrokken uit
het kalibreren mogelijk niet toepasbaar zijn op andere reconstructies. Hoofdstuk 7 behandelt de
maximalisatie van scherpte in detail, waaronder dit probleem.
(a) S (100) = 1, 00
(b) S (100) = 0, 94
(c) S (100) = 0, 83
Figuur 5.7: Genormaliseerde scherpte voor verschillende reconstructies.
5.3
Optimalisatiealgoritmes voor het CT-probleem
In sectie 5.1.3 beargumenteerden we reeds de keuze voor een stochastisch zoekalgoritme. We
bespreken nu twee algoritmes die aan dit criterium voldoen en in de literatuur voorkomen als
algoritmes om CT-scanners te kalibreren [3][11][13]. Het zijn algoritmes die steeds één objectieffunctie in rekening brengen. Hoofdstuk 8 zal de combinatie van objectieffuncties bespreken.
5.3.1
Nelder-Mead algoritme
De Nelder-Mead procedure is een heuristische, directe zoekmethode en staat beschreven in algoritme 5.1 (algoritme 5.2 beschrijft een hulpfunctie) [23]. De procedure geldt voor het minimaliseren van een objectieffunctie c in functie van Ω. In het geval van maximalisatie kan men
eenvoudigweg werken met 1/c(Ω) zonder het algoritme aan te passen. Een implementatie in
C++ is te vinden in de klasse Simplex.hh.
Werking
Stel dat de te optimaliseren parameterverzameling Ω bestaat uit n parameters pi . Het NelderMead algoritme zal dan een polytoop of simplex construeren met n+1 hoekpunten in n dimensies
(een polytoop is een uitbreiding van het veelvlak in meer dimensies). Elk hoekpunt komt overeen
met een andere oplossing Ωj , j ∈ {1, . . . , n + 1}. Het Nelder-Mead algoritme zal per iteratie de
hoekpunten van het polytoop verplaatsen in de hoop oplossingen te bekomen die de objectieffunctie nog verder minimaliseren. Het algoritme zal elke operatie uitvoeren op alle parameters
van een oplossing Ωj , en alle andere oplossingen mee in rekening brengen. Om het onderscheid
te maken met een iteratie k uit de iteratieve reconstructiealgoritmes, noemen we een iteratie
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
(a) reflectie
(b) expansie
(c) contractie
31
(d) contractie
Figuur 5.8: De vier basisoperaties van het Nelder-Mead algoritme.
van een optimalisatiealgoritme een tableau t. Initialisatie gebeurt doorgaans met n + 1 initiële
schattingen.
Het heeft weinig zin elke stap van algoritme 5.1 in woorden te beschrijven. We houden het bij het
illustreren van enkele belangrijke operaties. Bij aanvang van elk tableau worden de oplossingen
gerangschikt van minimale naar maximale kost, waarbij Ωn+1 de oplossing met maximale (en dus
slechtste) kost is. Lijn 4 van het algoritme toont de berekening van het zwaartepunt Ω van de n
beste oplossingen. Het is nu de bedoeling Ωn+1 te vervangen. Lijn 5 berekent het reflectiepunt Ωr
(figuur 5.8a). Lijn 10 toont de berekening van het expansiepunt Ωe (figuur 5.8b). Lijn 18 berekent
het uitwendig contractiepunt Ωc (figuur 5.8c). Lijn 25 berekent het inwendig contractiepunt Ωcc
(figuur 5.8d). Naargelang de kost verbonden met Ωr , Ωe , Ωc en Ωcc , zal Ωn+1 door een van deze
oplossingen vervangen worden.
Het Nelder-Mead algoritme gebruikt voor reflectie, expansie, contractie en inkrimping respectievelijk de parameters ρ, χ, γ en σ. Deze parameters moeten aan volgende voorwaarden voldoen:
ρ>0
(5.10)
χ>1
(5.11)
χ>ρ
(5.12)
0<γ<1
(5.13)
0<σ<1
(5.14)
In optimalisatieproblemen waar aanvankelijk weinig bekend is over het gedrag van de parameters
en objectieffuncties, kiest men meestal ρ = 1, χ = 2, γ = 1/2 en σ = 1/2 [24]. Het zijn deze
waarden die ook wij zullen hanteren.
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
32
Algoritme 5.1: Nelder-Mead procedure
1
2
3
4
Initialiseer Simplex = {Ω1 , . . . , Ωn+1 } met initiële oplossingen Ωj
while stopcriterium niet vervuld do
Rangschik Simplex = {Ω1 , . . . , Ωn+1 } zodat c(Ω1 ) ≤ c(Ω2 ) ≤ . . . ≤ c(Ωn+1 )
n
P
pi
Ω←
n
i=1
5
6
7
8
9
10
11
12
13
14
15
Ωr ← Ω + ρ(Ω − Ωn+1 )
if c(Ω1 ) ≤ c(Ωr ) < c(Ωn ) then
Simplex ← {Simplex\{Ωn+1 }} ∪ {Ωr }
else
if c(Ωr ) < c(Ω1 ) then
Ωe ← Ω + χ(Ωr − Ω)
if c(Ωe ) < c(Ωr ) then
Simplex ← {Simplex\{Ωn+1 }} ∪ {Ωe }
else
Simplex ← {Simplex\{Ωn+1 }} ∪ {Ωr }
end
else
if c(Ωr ) < c(Ωn+1 ) then
Ωc ← Ω + γ(Ωr − Ω)
if c(Ωc ) < c(Ωr ) then
Simplex ← {Simplex\{Ωn+1 }} ∪ {Ωc }
else
Simplex ← Shrink(Simplex,σ)
end
16
17
18
19
20
21
22
23
30
else
Ωcc ← Ω − γ(Ω − Ωn+1 )
if c(Ωcc ) < c(Ωn+1 ) then
Simplex ← {Simplex\{Ωn+1 }} ∪ {Ωcc }
else
Simplex ← Shrink(Simplex,σ)
end
31
end
24
25
26
27
28
29
end
32
33
34
end
end
Stopcriterium
We stoppen het optimalisatiealgoritme wanneer de kost van de optimale parameterverzameling
Ωt in tableau t onvoldoende verschilt van het minimum c(Ωt−1 ) uit het vorige tableau:
c(Ωt−1 ) − c(Ωt )
<
c(Ωt )
(5.15)
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
33
Algoritme 5.2: Nelder-Mead procedure: Shrink(Simplex, σ)
1
2
3
4
5
6
Simplex0 ← {Ω1 }
for j = 2 to n + 1 do
Ω0j ← Ω1 + σ(Ωj − Ω1 )
Simplex0 ← Simplex0 ∪ {Ω0j }
end
return(Simplex0 )
We zullen goede resultaten bekomen met = 0.0001. Om lokale minima te vermijden, stoppen
we in de praktijk het algoritme pas wanneer voor strikt opeenvolgende tableaus t een aantal keer
aan vergelijking 5.15 voldaan is.
5.3.2
Genetisch algoritme
Genetische algoritmes leveren goede benaderingen van globale extrema, zelfs in hoogdimensionale oplossingenruimtes [25]. Deze klasse van algoritmes bootsen het natuurlijke proces van
evolutie na. Een stappenplan is te vinden in algoritme 5.3 [3]. In wat volgt, beschrijven we het
optimalisatiealgoritme om een objectieffunctie te maximaliseren. Minimalisatie is mogelijk door
te werken met 1/c(Ω) in plaats van c(Ω). Een implementatie in C++ is te vinden in de klasse
Population.hh, met Genetics.hh als ouderklasse (zie Bijlage C).
Initiële populatie
Elke oplossing Ωj wordt een individu van de populatie Tt genoemd: j ∈ {1, . . . , |Tt |}. De grootte
|Tt | van de populatie wordt vooraf bepaald en blijft gelijk voor alle t (| · | staat hier symbool
voor de kardinaliteit).
Concreet geven we voor t = 0 een intiële schatting Ω1 op. Ω2 tot en met Ω|T0 | stellen we gelijk
aan Ω1 , maar we veranderen hun individuele parameters pi door er een willekeurig getal bij
op te tellen. Dit willekeurig getal is afkomstig uit een normale verdeling met gemiddelde 0 en
standaardafwijking σi . Bij gebrek aan kennis over σi , blijkt σi = 1 een bruikbare waarde.
Voortplanting
Om een nieuwe populatie Tt+1 samen te stellen, combineren we willekeurig individuele parameters pi van een gekozen ‘ouderpaar’. Dit ouderpaar bestaat uit twee individuen die in populatie
Tt goed gepresteerd hebben in termen van het maximaliseren van c. Deze voortplantingsstap
wordt geïllustreerd op figuur 5.9.
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
34
Algoritme 5.3: Genetisch algoritme
1
2
3
4
5
6
Creër een initiële populatie T0 . Zet t = 0.
Bereken voor elk individu Ωj de objectiefwaarde c(Ωj ).
Indien voldaan aan stopcriterium: individu met grootste c(Ωj ) is optimum. Stop.
Zet t = t + 1 en creër een nieuwe populatie Tt (voortplanting met elitarisme).
Pas mutatie toe op de nieuwe populatie Tt .
Bereken voor elk individu Ωj de objectiefwaarde c(Ωj ). Ga naar stap 3.
Figuur 5.9: Principe van voortplanting: elke kleur stelt een individuele parameter voor.
Om ervoor te zorgen dat niet alleen de beste individuen het voortplantingsproces domineren,
kennen we aan elk individu een kans toe dat het geselecteerd wordt als ouder. De kans dat Ωj
als ouder geselecteerd wordt, stellen we gelijk aan sj [26]:
sj =
c(Ωj )
|T
Pt |
(5.16)
c(Ωl )
l=1
Elk individu kan dus als ouder geselecteerd worden, maar deze kans neemt wel toe naarmate de
objectieffunctie voor het individu een hogere waarde teruggaf.
Mutatie
De individuen die ontstaan uit de voortplantingsstap, ondergaan mutatie. Het muteren gebeurt
eveneens door een willekeurig getal op te tellen bij de individuele parameters pi . Dit willekeurig
getal is afkomstig uit een normale verdeling met gemiddelde 0 en standaardafwijking σi . Bij
gebrek aan kennis over σi , blijkt σi = 1 een bruikbare waarde. Mutatie van een parameter
vindt slechts met een bepaalde kans plaats. Deze kans wordt best groter gemaakt naarmate
men met een ruwere initiële schatting werkt: 1/n of 1/(2n) zijn dan bruikbaar. Naarmate de
individuen naar een optimum streven, wordt de kans om een parameter te muteren best verkleind,
bijvoorbeeld naar 1/|Tt |. Deze laatste ingreep zal een invloed hebben op de nauwkeurigheid van
de finale resultaten.
De keuze om met een normale verdeling te werken [3], is ingegeven door de kennis van σi die
men heeft naarmate men het kalibreren herhaalt en van de bekomen resultaten gemiddelde en
standaardafwijking bijhoudt. Strikt genomen zou mutatie echter moeten plaatsvinden op bitre-
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
35
presentaties van de individuele parameters [25]. De methode met σi neemt echter de moeilijkheid
van het bepalen van een goede bitrepresentatie weg.
Ten slotte bewaren we uit elke populatie een aantal van de beste individuen uit de vorige populatie, zonder hun individuele parameters te muteren. Op die manier zorgen we ervoor dat er
geen goede oplossingen verloren gaan als de voortplanting niets oplevert. Dit principe noemen
we elitarisme.
In sectie 6.3.3 zal het genetisch algoritme verder aangepast worden naargelang de specifieke
problemen die zich zullen stellen bij het kalibreren van CT-scanners. De populatiegrootte zetten
we op 6 tot 8 individuen per te optimaliseren parameter [3].
Stopcriterium
We hanteren hetzelfde stopcriterium als bij de Nelder-Mead procedure: als de waarde van c voor
opeenvolgende tableaus t nog maar nauwelijks toeneemt, stoppen we het genetisch algoritme en
wordt het best presterende individu uit het laatste tableau als ‘winnaar’ gekozen.
5.4
Conclusie
Dat scherpte beter zou presteren dan entropie [13], zet ons ertoe aan enkel de projectiefout en
scherpte in volgende hoofdstukken als objectieffuncties te bespreken. Welke parameters laten
zich het best met welke objectieffunctie optimaliseren, en na hoeveel iteraties van het iteratief
reconstructiealgoritme kunnen we betrouwbaar de waarde van de objectieffuncties evalueren? Is
het echt nodig alle parameters samen te optimaliseren, zoals bijvoorbeeld in het onderzoek van
Sawall et al. [3] en Wein et al. [15]?
Wat de optimalisatiealgoritmes betreft, gaat er al een lichte voorkeur uit naar het genetisch
algoritme. Het is immers een bekend probleem dat de Nelder-Mead procedure, in tegenstelling
tot het genetisch algoritme, niet uit lokale extrema kan ontsnappen [24]. Lokale extrema vormen
in ons onderzoek een reëel gevaar [16] en het is met andere woorden te riskant om blindelings
het Nelder-Mead algoritme toe te passen op het kaliberen van CT-scanners, zoals tot nu toe het
geval was in eerdere studies [11][13]. In de volgende hoofdstukken zal blijken dat de projectiefout
en scherpte enkel onder bepaalde voorwaarden geen lokale extrema bezitten.
Toch zullen deze vaststellingen het Nelder-Mead algoritme niet geheel overbodig maken. De
Nelder-Mead procedure kan sneller dan het genetisch algoritme steile hellingen van de objectieffunctie nemen, in het geval er nauwelijks lokale extrema zijn [23].
De keuze voor een genetisch algoritme in het onderzoek van Sawall et al. [3] blijkt dus over
het algemeen gegrond, maar zij optimaliseren alle geometrische parameters tegelijkertijd. Dit
verhoogt de complexiteit van het probleem aanzienlijk. Wij willen onder meer nagaan of we
HOOFDSTUK 5. OPTIMALISATIETECHNIEKEN
36
tot een optimum kunnen komen door parameters afzonderlijk te optimaliseren. Hoe dichter we
bij het optimum komen, hoe bruikbaarder Nelder-Mead alsnog zal worden. De twee algoritmes
sluiten elkaar dus niet uit, maar moeten elkaar aanvullen.
Hoofdstuk 6
Minimalisatie projectiefout
We bespreken hier het kalibreren van CT-scanners aan de hand van het minimaliseren van de
(k)
(k)
projectiefout EΩ . De objectieffunctie EΩ kwam reeds formeel aan bod in sectie 5.2.2. Het
gebruikte optimalisatiealgoritme zal in de meeste gevallen het genetisch algoritme zijn, dat in dit
hoofdstuk verdere aanpassingen zal ondergaan om tegemoet te komen aan specifieke problemen
die zich stellen in de context van CT.
6.1
Invloed parameters op projectiefout
We willen nagaan wat de invloed van de geometrische parameters op de projectiefout is. Wanneer
we het hebben over parameterverzameling Ω, dan bedoelen we de werkelijke parameters van het
ˆ staat voor onze initiële schatting uit tabel 4.5.
systeem zoals getoond in tabel 4.3. Ω
We maken nog steeds gebruik van SIRT. Het aantal iteraties van dit reconstructiealgoritme
duiden we aan met k. Een iteratie van het optimalisatiealgoritme noemen we een tableau t.
6.1.1
Horizontale uitwijking detector ∆u
Uit de studie van Smekal et al. [27] weten we dat de positie van de detector de grootste invloed
op de reconstructiekwaliteit heeft. Dit wordt in ons onderzoek visueel gestaafd door figuur 4.6.
We verwachten deze nefaste invloed uiteraard ook in de projectieruimte waar te nemen, in de
vorm van een grote projectiefout.
We zetten een experiment op waarin we reconstructies uitvoeren voor verschillende waarden van
∆u, met alle andere parameters vastgehouden op hun correcte waarden uit Ω. Figuur 6.1 toont
(k)
de behaalde resultaten. Merk op dat we de reciproke projectiefout 1/EΩ visualiseren, die we
dus willen maximaliseren. Ook de scherpte wordt al getoond, maar een bespreking ervan volgt
37
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) k = 2
(c) k = 20
38
(b) k = 10
(d) k = 100
Figuur 6.1: Invloed van ∆u op objectieffuncties in aanwezigheid van correcte parameters (Ω).
pas in Hoofdstuk 7. De waarden werden tevens genormaliseerd met de waarden bij een correcte
parameter, in dit geval bij ∆u = −3 mm.
Het valt op dat de reciproke projectiefout meteen (i.e. na weinig iteraties) een maximum oplevert
bij de beoogde waarde ∆u = −3 mm. Er zijn geen lokale extrema aanwezig. Naarmate het aantal
iteraties toeneemt, wordt het verschil tussen de projectiefout bij ∆u = −3 mm en de projectiefout
bij waarden in zijn directe omgeving, groter.
Bij kennis van de parameters ∆v, θ en SO, kunnen we ∆u uitermate snel optimaliseren op basis
van de projectiefout. Gezien het hier om één te optimaliseren parameter gaat, is het Nelder-Mead
algoritme aangewezen (ook al omdat er slechts één extremum is).
Interessanter is natuurlijk in hoeverre ∆u zich laat optimaliseren in de aanwezigheid van andere
parameters die (nog) niet optimaal zijn. Als we ook dan ∆u onafhankelijk van de rest kunnen
optimaliseren, daalt de complexiteit van de zoekruimte aanzienlijk. We herhalen hetzelfde experiment met variabele ∆u, maar we gebruiken voor de andere parameters de waarden uit de
ˆ Het resultaat is samengevat in figuur 6.2.
initiële schatting Ω.
Het resultaat is gelijkaardig aan dat uit figuur 6.1. Dit betekent twee dingen. Ten eerste is
de invloed van ∆u op de projectiefout erg groot. Ten tweede laat ∆u zich onafhankelijk van
de andere parameters optimaliseren. Nagenoeg perfecte alignatie van de geprojecteerde data
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) k = 2
(c) k = 20
39
(b) k = 10
(d) k = 100
ˆ
Figuur 6.2: Invloed van ∆u op objectieffuncties in aanwezigheid van initiële parameters (Ω).
met de meting g is dus mogelijk in de aanwezigheid van andere parameters die ook een fout in
de projectieruimte introduceren. Verschuivingen van de detector in zijn vlak blijken vroeg een
richtinggevende invloed op de projectiefout te hebben. Een snelle optimalisatie is mogelijk.
6.1.2
Verticale uitwijking detector ∆v
De verticale uitwijking ∆v is geen verschuiving van de detector in zijn vlak, en we verwachten dus
dat de projectiefout minder bruikbaar zal zijn om deze parameter te optimaliseren. Wederom
testen we voor verschillende waarden van ∆v de invloed op de projectiefout, in aanwezigheid van
de optimale parameters uit Ω. Het resultaat is te zien op figuur 6.3.
De reciproke projectiefout vertoont pas een maximum bij ∆v = 6 mm na k = 20. Dit maximum
bevindt zich evenwel in de nabijheid van lokale maxima. Het is pas voor k = 100 dat het
maximum bij de beoogde waarde ∆v = 6 mm zich duidelijk aftekent. Er zijn nog steeds lokale
maxima, maar die bevinden zich bij een veel kleinere reciproke projectiefout. De projectiefout
is in dit geval een goede schatter voor ∆v, al moet de projectiefout na een voldoende aantal
iteraties geëvalueerd worden.
ˆ laat ∆v zich niet schatten op basis van de projectieIn aanwezigheid van de initiële parameters Ω,
fout, zoals blijkt uit figuur 6.4. Omdat ∆u een grote invloed had op de projectiefout, hebben we
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) k = 2
(c) k = 20
40
(b) k = 10
(d) k = 100
Figuur 6.3: Invloed van ∆v op objectieffuncties in aanwezigheid van correcte parameters (Ω).
(a) k = 2
(b) k = 100
ˆ
Figuur 6.4: Invloed van ∆v op objectieffuncties in aanwezigheid van initiële parameters (Ω).
het experiment ook herhaald met ∆u = −3 mm. Het resultaat is te zien op figuur 6.5: er treedt
geen verbetering op. De kennis van ∆u alleen is dus onvoldoende om ∆v in de aanwezigheid van
foute parameters te schatten op basis van de projectiefout: er tekent zich geen extremum af in
de buurt van het beoogde ∆v = 6 mm.
We vermoeden reeds dat parameters die de vergrotingsfactor M = SO/SD beïnvloeden (∆v en
SO) zich enkel samen zullen laten optimaliseren met betrouwbare kennis van de andere parameters.
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) k = 2
41
(b) k = 100
ˆ met ∆u = −3 mm.
Figuur 6.5: Invloed van ∆v in aanwezigheid van initiële parameters (Ω),
(a) k = 2
(c) k = 20
(b) k = 10
(d) k = 100
Figuur 6.6: Invloed van θ op objectieffuncties in aanwezigheid van correcte parameters (Ω).
6.1.3
Detectortilt θ
Figuur 6.6 toont voor aanpassingen van θ bij optimale parameters Ω de resulterende projectiefout.
Er tekent zich een maximum af voor de reciproke projectiefout in de buurt van het beoogde
optimum θ = −9◦ , maar dit maximum komt pas voor k = 20 op de juiste plaats te liggen.
In aanwezigheid van de initiële parameters, zie figuur 6.7, tekent er zich geen maximum af bij
θ = −9◦ , ook niet na een groot aantal iteraties k. Het is pas wanneer we ∆u aanpassen tot zijn
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) k = 2
42
(b) k = 100
ˆ
Figuur 6.7: Invloed van θ op objectieffuncties in aanwezigheid van initiële parameters (Ω).
(a) k = 2
(b) k = 100
ˆ met ∆u = −3 mm.
Figuur 6.8: Invloed van θ in aanwezigheid van initiële parameters (Ω),
optimale waarde, θ een maximum krijgt bij θ = −9◦ (zie figuur 6.8). Toch moeten we ook hier
wachten tot een hondertal iteraties om de projectiefout betrouwbaar te evalueren. Het toont wel
het belang aan van ∆u als cruciale parameter.
6.1.4
Afstand stralingsbron tot object SO
Figuur 6.9 toont voor aanpassingen van SO bij optimale parameters Ω de resulterende projectiefout. Er tekent zich pas na ongeveer k = 100 een maximum af voor de reciproke projectiefout
bij het beoogde optimum SO = 151 mm.
In aanwezigheid van de initiële parameters, zie figuur 6.10, tekent er zich geen maximum af bij
SO = 151 mm (integendeel), ook niet voor een groot aantal iteraties k. Ook het vasthouden van
∆u = −3 mm helpt niet (zie figuur 6.11). We zien hier hetzelfde fenomeen opduiken als in de
situatie met ∆v, de andere parameter die de vergrotingsfactor M rechtstreeks aantast.
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) k = 2
43
(b) k = 100
Figuur 6.9: Invloed van SO op objectieffuncties in aanwezigheid van correcte parameters (Ω).
(a) k = 2
(b) k = 100
ˆ
Figuur 6.10: Invloed van SO op objectieffuncties in aanwezigheid van initiële parameters (Ω).
(a) k = 2
(b) k = 100
ˆ met ∆u = −3 mm.
Figuur 6.11: Invloed van SO in aanwezigheid van initiële parameters (Ω),
6.2
Gedrag projectiefout
Figuur 6.12 toont hoe de projectiefout evolueert naarmate het aantal iteraties toeneemt (de
waarden werden genormaliseerd met de projectiefout bij optimale parameters voor k = 2). Ten
eerste merken we hoe de projectiefout een monotoon dalende functie is. Dit geldt voor om het
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
44
Figuur 6.12: Globaal gedrag van de projectiefout.
even welke reconstructie. Ten tweede stagneert de projectiefout voor reconstructies met foute
parameters snel. Ter vergelijking: na 100 iteraties bedraagt de projectiefout van de reconstructie
met optimale parameters nog maar 2,4% van de projectiefout na 2 iteraties; voor de reconstructie met initiële parameters is dat 79,4%. De projectiefout van de reconstructie met optimale
parameters bedraagt na 800 iteraties 1,6% van de projectiefout bij 100 iteraties; bij reconstructie
met initiële parameters is dat 99,8%.
Bovendien kunnen er snijpunten optreden tussen de curves van de projectiefout voor verschillende
reconstructies. Deze snijpunten doen zich onder meer voor wanneer de projectiefout voor een
correct ingeschatte parameter pas minimaal wordt na voldoende iteraties, zie bijvoorbeeld figuur
6.8 voor de invloed van θ op de projectiefout. We merken wel op dat bij een voldoende groot
aantal iteraties de optimale reconstructie altijd de minimale projectiefout zal hebben.
6.3
6.3.1
Optimalisatie
Directe optimalisatie
Uit voorgaande observaties wordt duidelijk dat het geen zin heeft alle parameters tegelijkertijd
te optimaliseren door de projectiefout na een beperkt aantal iteraties (e.g. k = 20) te evalueren.
De parameters ∆v en SO laten zich (samen) pas na veel iteraties optimaliseren, en enkel in de
aanwezigheid van correcte ∆u en θ. Deze laatste parameters kunnen al na veel minder iteraties
correct geschat worden, onafhankelijk van ∆v en SO. Alle parameters tegelijkertijd optimaliseren,
zoals gebruikelijk in de literatuur [3][15], vereist dus een te groot aantal iteraties van het iteratief
reconstructiealgoritme en gebruikt een te grote zoekruimte.
Toepassing van het Nelder-Mead algoritme op alle parameters, met evaluatie van de projectiefout
na 200 iteraties, leverde niet het optimum op. De resultaten waren ∆u = −2, 95 mm, ∆v =
3, 73 mm, θ = −5, 74◦ en SO = 149, 62 mm. Het algoritme is niet uit een lokaal minimum
kunnen ontsnappen, want het optimum Ω had een nog veel kleinere projectiefout dan deze
parameterverzameling. Het Nelder-Mead algoritme was bovendien niet in staat de parameter ∆u
te verfijnen tot −3 mm, wat in aanwezigheid van de andere parameters een kleinere projectiefout
had opgeleverd.
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
(a) Nelder-Mead
(b) genetisch
45
(c) combinatie
Figuur 6.13: (a) en (b): directe optimalisatie, (c): stapsgewijze optimalisatie
Het genetisch algoritme lost deze tekortkomingen op, weliswaar tegen de prijs van meer reconstructies. We behaalden het resultaat ∆u = −2, 981 mm, ∆v = 5, 882 mm, θ = −9, 784◦ en
SO = 152, 014 mm. Figuur 6.13 vergelijkt deze parameters met het resultaat van het NelderMead algoritme. Het is duidelijk dat de parameterschatting van het genetisch algoritme niet
alleen dichter bij het optimum ligt (in termen van afstand), maar dat ook de reconstructie visueel beter is. Het resultaat zou nog verfijnd kunnen worden, door het genetisch algoritme met
een kleinere standaardafwijking op de parameters te laten werken tijdens de mutatiestap. Het
is tevens mogelijk het Nelder-Mead algoritme op deze schatting toe te passen. Het meest nabijgelegen minimum, waarin de procedure hoogstwaarschijnlijk verzeild zal raken, zal immers het
globale minimum zijn.
6.3.2
Stapsgewijze optimalisatie
Het grote aantal iteraties k ≥ 200 en het aanzienlijke aantal reconstructies per tableau van
het genetisch algoritme (typisch 6 tot 8 reconstructies per parameter [3]), valt makkelijk te
omzeilen door de resultaten uit sectie 6.1 te gebruiken. We kunnen op basis hiervan een methode
voorstellen die parameters afzonderlijk optimaliseert en het aantal iteraties beperkt houdt. Het
stappenplan ziet er als volgt uit:
1. Optimaliseer ∆u met het Nelder-Mead algoritme. Zet k = 2.
2. Optimaliseer θ met het Nelder-Mead algoritme. Zet k = 100.
3. Optimaliseer ∆v en SO met het genetisch algoritme. Zet |Tt | = 12 en k = 200.
Het voordeel van het Nelder-Mead algoritme in stap 1 en 2 is de goede convergentie in deze
omstandigheden. De Nelder-Mead procedure is immers in staat om bij functies met één extremum
zeer snel steile hellingen te nemen [24]. Uit figuur 6.2 weten we dat dit voor ∆u het geval is; bij
θ is controle vereist of het algoritme niet in een lokaal extremum verzeild raakt. In de praktijk
bleek dat niet het geval.
HOOFDSTUK 6. MINIMALISATIE PROJECTIEFOUT
46
De keuze voor het genetisch algoritme in stap 3 heeft als doel uit lokale extrema te onstnappen,
die bij het aanpassen van ∆v en SO veelvuldig opduiken.
6.3.3
Aanpassingen genetisch algoritme
Keren we terug naar het experiment waarbij we alle parameters tegelijkertijd optimaliseerden met
k = 200, dan viel het op dat ∆u zeer snel naar zijn optimum −3 mm evolueerde. Deze waarde
over alle volgende tableaus vastgehouden. In de implementatie van het genetisch algoritme
zorgden we ervoor dat in dat geval, i.e. wanneer een parameter uit de best presterende oplossing
lang genoeg constant blijft, de populatiegrootte |Tt | dynamisch kleiner wordt en de bewuste
parameter als gekalibreerd beschouwd wordt. Dit verkleint de zoekruimte en dus de complexiteit
van het probleem.
Een maatregel om uit lokale extrema te ontsnappen, bestaat erin voor elke parameter uit de
best presterende oplossingen bij te houden in welke richting ze evolueren. Er zullen bijkomende
individuen gecreëerd worden die deze evolutie verderzetten, maar ook individuen die er bewust
van afwijken om sneller nieuwe richtingen in de zoekruimte te onderzoeken.
6.4
Conclusie
Het minimaliseren van de projectiefout blijkt een robuuste methode te zijn om de correcte geometrische parameters in te schatten. Zolang we het aantal iteraties k voldoende groot maken,
zal het genetisch algoritme inderdaad streven naar het beoogde optimum. In de praktijk is een
grote waarde voor k, zeker tijdens het kalibreren, vaak inefficiënt omwille van de verhoogde
reconstructietijd.
Uit een parameteranalyse is echter gebleken dat bepaalde parameters zich onafhankelijk van de
andere laten optimaliseren. Binnen redelijke tijd kan op deze manier (stapsgewijs) een goede
schatting gevonden worden, door enerzijds de zoekruimte te verkleinen en anderzijds het aantal
iteraties van het iteratief reconstructiealgoritme te verminderen. Een samenvatting is te zien op
figuur 6.14.
Figuur 6.14: Stappenplan voor het kalibreren.
Hoofdstuk 7
Maximalisatie scherpte
Het minimaliseren van de projectiefout gebeurt in de projectieruimte. In de beeldruimte heeft
(k)
de scherpte SΩ , gedefinieerd als de L2 -norm van de beeldgradiënt, veelbelovende resultaten
opgeleverd [13]. We willen de prestaties van deze objectieffunctie vergelijken met die van de
projectiefout. Scherpte kwam reeds formeel aan bod in sectie 5.2.3.
7.1
Invloed parameters op scherpte
Net als in Hoofdstuk 6, dat de minimalisatie van de projectiefout behandelde, gebruiken we hier
ˆ voor onze initiële schatting (tabel
de notatie Ω voor de correcte parameters (tabel 4.3) en Ω
4.5). Voor reconstructie gebruiken we eveneens SIRT. Op de grafieken in sectie 6.1 werd reeds
de scherpte getoond. We zullen in dit hoofdstuk dan ook naar die figuren verwijzen.
7.1.1
Horizontale uitwijking detector ∆u
Op figuur 6.1 zagen we dat de scherpte een maximum vertoont bij de optimale waarde voor ∆u,
in de aanwezigheid van andere parameters die reeds optimaal zijn. Dit maximum is echter lokaal.
Globale maxima doen zich voor bij waarden van ∆u die een erg grote projectiefout introduceren.
Figuur 7.1 toont een voorbeeld van een reconstructie volgens Ω, maar met ∆u = 1 mm. Deze
reconstructie heeft volgens de grafieken op figuur 6.1 een grotere numerieke scherpte dan de
optimale reconstructie. De berekening van de scherpte verloopt via randdetectie: de dubbele
randen zijn voldoende afgetekend om als ‘scherper’ beschouwd te worden.
ˆ duikt
Op figuur 6.2, waar ∆u variabel wordt gemaakt in aanwezigheid van parameters uit Ω,
een gelijkaardig fenomeen op. We zien dat na 100 iteraties de maximale scherpte elders komt te
liggen, op ∆u = −1 mm. Ook hier spelen de randen een rol. Figuur 7.2 toont reconstructies
ˆ
volgens Ω met bijhorende randdetecties. Figuur 7.3 doet hetzelfde voor reconstructies volgens Ω,
47
HOOFDSTUK 7. MAXIMALISATIE SCHERPTE
(a) reconstructie
48
(b) randdetectie
(c) randdetectie
Figuur 7.1: Scherper geëvalueerde reconstructie dan het optimum.
(a) k = 20
(b) randdetectie
(c) k = 100
(d) randdetectie
Figuur 7.2: Reconstructie volgens optimum Ω.
(a) k = 20
(b) randdetectie
(c) k = 100
(d) randdetectie
ˆ met ∆u = −1 mm.
Figuur 7.3: Reconstructie volgens schatting Ω,
met ∆u = −1 mm. De dubbele randen zijn na 100 iteraties voldoende doortekend om scherper
geëvalueerd te worden dan het optimum Ω. De doortekening was na k = 20 nog onvoldoende.
7.1.2
Verticale uitwijking detector ∆v
We weten dat de numerieke scherpte toeneemt met een toenemende vergrotingsfactor M =
SO/SD [13]. Gezien we de afstand SD beïnvloeden door ∆v aan te passen, is het inderdaad niet
verwonderlijk dat we op figuren 6.3, 6.4 en 6.5 een stijgende trend in scherpte waarnemen voor
toenemende ∆v. Op figuur 6.3 tekent er zich na voldoende iteraties wel een lokaal maximum af
bij het optimum. Figuur 7.4 verduidelijkt dit effect. Laten we ∆v toenemen van zijn optimale
waarde ∆v = 6 mm naar ∆v = 9 mm, dan nemen we inderdaad vergroting waar, maar de wazige
HOOFDSTUK 7. MAXIMALISATIE SCHERPTE
(a) ∆v = 6 mm
49
(b) ∆v = 9 mm
(c) ∆v = 12 mm
Figuur 7.4: Reconstructie volgens optimum Ω (k = 100).
(a) θ = −9◦
(b) θ = −45◦
Figuur 7.5: Reconstructie volgens optimum Ω (k = 100).
randen compenseren het effect dat vergroting op de numerieke scherpte heeft. Pas bij ∆v = 12
mm is de vergroting voldoende en zijn de randen scherp genoeg om van een toename in numerieke
scherpte te spreken. Erg lokaal kan scherpte dus richtinggevend zijn om ∆v te bepalen, maar
een algemene regel valt niet af te leiden.
7.1.3
Detectortilt θ
Uit figuur 6.6 leiden we af dat de scherpte vroeger dan de projectiefout (k = 10) de optimale
waarde van de detectortilt θ aanduidt, maar deze vaststelling geldt enkel wanneer alle andere
parameters al correct geschat zijn. Scherpte is niet richtinggevend voor een correcte schatting
van θ wanneer de andere parameters nog niet hun optimale waarden bezitten, zoals blijkt uit
figuren 6.6 en 6.7.
We mogen niet uit het oog verliezen dat bij extreme waarden van θ ook de vergrotingsfactor
aangetast wordt, zij het in de aanwezigheid van artefacten (dubbele randen). Figuur 7.5 toont
een reconstructies volgens Ω, maar met θ = −45◦ .
7.1.4
Afstand stralingsbron tot object SO
Naar analogie met de bespreking van ∆v moeten we ook hier besluiten dat SO niet geschat kan
worden aan de hand van de numerieke scherpte. De vergrotingsfactor M = SO/SD speelt immers
HOOFDSTUK 7. MAXIMALISATIE SCHERPTE
50
Figuur 7.6: Globaal gedrag van de scherpte.
een allesoverheersende rol.
7.2
Gedrag scherpte
In tegenstelling tot de projectiefout, treedt er in het geval van scherpte geen snelle stagnatie op
bij reconstructies met foute parameters. Belangrijker is dat suboptimale reconstructies scherper
geëvalueerd kunnen worden dan het optimum. Dit resultaat geldt zelfs wanneer we de vergrotingsfactor M buiten beschouwing laten. Zo toonde figuur 7.1 een reconstructie volgens Ω, maar
met ∆u = −1 mm. De geïntroduceerde artefacten (dubbele randen) zorgen ervoor dat deze
reconstructie scherper beoordeeld wordt dan de optimale.
Scherpte heeft geen globaal maximum bij de werkelijke parameters (zie figuur 6.1) en vertoont
lokale maxima (zie figuur 6.2). Hierdoor schaadt scherpte twee belangrijke voorwaarden voor
een goede objectieffunctie.
7.3
7.3.1
Optimalisatie
Directe optimalisatie
Alle parameters optimaliseren aan de hand van scherpte is onmogelijk, door de rol van de vergrotingsfactor M en de onmogelijkheid om met scherpte θ te bepalen zonder kennis van ∆v en
SO. ∆v en SO laten zich dan weer niet schatten zonder kennis van θ, waardoor enkel nog ∆u
overblijft om met scherpte te optimaliseren.
7.3.2
Stapsgewijze optimalisatie
Het stappenplan verschilt nauwelijks met het plan dat we voor de projectiefout uiteenzetten:
1. Optimaliseer ∆u aan de hand van scherpte, met het Nelder-Mead algoritme. Zet k = 2.
HOOFDSTUK 7. MAXIMALISATIE SCHERPTE
51
2. Optimaliseer θ aan de hand van de projectiefout, met het Nelder-Mead algoritme. Zet
k = 100.
3. Optimaliseer ∆v en SO aan de hand van de projectiefout, met het genetisch algoritme. Zet
|Tt | = 12 en k = 200.
Belangrijk is dat in stap 1 een bovengrens op de projectiefout gezet moet worden, om de scherpte
niet naar de foute maxima te sturen. Er zit weinig intuïtie achter het kiezen van zo’n bovengrens,
ˆ gekozen worden.
maar in de praktijk kan de bovengrens van de initiële schatting Ω
Op die manier zijn goede schattingen voor ∆u te bekomen, die gelijk zijn aan de resultaten die
we met de projectiefout behaalden. In het geval van driedimensionale reconstructies (zie sectie
10.5) zal scherpte het bijkomende voordeel hebben dat het een paar iteraties vroeger dan de
projectiefout de detectorparameters in de juiste richting duwt.
7.4
Conclusie
Scherpte heeft weinig voordelen vergeleken met het gebruik van de projectiefout. Toch zal in het
geval van driedimensionale reconstructies blijken dat scherpte iets sneller dan de projectiefout
een uitspraak over de detectorparameters kan doen, met name de parameters die verschuivingen
van de detector in zijn vlak veroorzaken. Voor tweedimensionale reconstructies verschilt het
stappenplan nauwelijks van het vorige, zie figuur 7.7.
Het optimaliseren van alle parameters zal nooit met scherpte kunnen, omwille van de vergrotingsfactor M .
Figuur 7.7: Stappenplan voor het kalibreren.
Hoofdstuk 8
Gecombineerde objectieffuncties
Dit hoofdstuk behandelt hoe we tegelijkertijd de projectiefout en scherpte als objectieffuncties
kunnen gebruiken binnen de klasse van de genetische algoritmes. De voorgestelde methode
moet ons in staat stellen om reeds na een beperkt aantal iteraties een betrouwbare uitspraak te
doen over de richting die de geometrische parameters moeten uitgaan, indien we alle parameters
tegelijkertijd willen optimaliseren.
8.1
Probleemstelling
Het gebruik van de projectiefout als objectieffunctie bleek robuust, maar heeft een groot nadeel:
het aantal te gebruiken iteraties per reconstructie moet groot genoeg zijn, willen we alle parameters tegelijkertijd optimaliseren. In de praktijk betekent dit dat het kalibreren op basis van de
projectiefout ontzettend lang kan duren, afhankelijk van het gebruikte reconstructiealgoritme en
ˆ Scherpte als objectieffunctie beïnvloedde dan weer de vergrotingsfactor,
de initiële schatting Ω.
waardoor het nooit alleen gebruikt kan worden om alle parameters te optimaliseren.
Toch kunnen beide objectieffuncties elkaar bijsturen. Figuur 6.1 liet zien hoe de maximale
scherpte tegenover een zeer grote projectiefout staat. Wanneer we optimaliseerden op basis van
scherpte, moesten we daarom een bovengrens op de projectiefout zetten. In het ideale geval zou
dit vanzelf gaan. Op figuur 6.9 werd duidelijk dat scherpte erg lokaal toch richtinggevend kan
zijn om het optimum te bepalen van parameters die de vergroting beïnvloeden.
We willen nagaan hoe we beide objectieffuncties samen kunnen evalueren om alle parameters
tegelijkertijd al na een beperkt aantal iteraties richting een optimum te duwen. Tabel 8.1 verduidelijkt een mogelijk scenario.
52
HOOFDSTUK 8. GECOMBINEERDE OBJECTIEFFUNCTIES
ˆ1
Ω
ˆ2
Ω
∆u
∆v
θ
-2,97 mm
-2,97 mm
5,93 mm
6,81 mm
-7,50
-6,66
◦
◦
53
SO
1/E (10)
S (10)
1/E (28)
S (28)
159,52 mm
162,17 mm
0,079
0,089
1939,37
1931,72
0,40
0,39
2995,35
2987,88
Tabel 8.1: Situatie waarin scherpte de projectiefout kan bijsturen.
(a) Ω
ˆ1
(b) Ω
ˆ2
(c) Ω
Figuur 8.1: Reconstructies met k = 100.
8.2
Principe
ˆ 1 , die nochtans het dichtst in de buurt van het optimum ligt (tabel 4.3),
Uit de tabel blijkt dat Ω
ˆ 2 , voor k = 10. Voor k ≤ 10 kan er op basis van de
toch een grotere projectiefout oplevert dan Ω
projectiefout dus geen zinnig vergelijk gemaakt worden. Het blijkt echter dat de scherpte voor
ˆ 2 voor k = 10 al aanzienlijk beter was dan de scherpte voor Ω
ˆ 1 . Scherpte kan
reconstructie met Ω
ˆ 1 op basis van de projectiefout weerleggen
dus richtinggevend zijn en het tijdelijke voordeel van Ω
of op zijn minst in vraag stellen.
Men kan zich afvragen of het dan niet beter is om (in dit geval voor k = 10) meteen de scherpte te
optimaliseren. Zoals al mocht blijken uit figuur 6.9 werkt scherpte de vergroting in de hand, die
we dan weer zouden kunnen tegenwerken door de projectiefout mee in beschouwing te nemen. We
willen beide objectieffuncties dus samen evalueren om dit soort situaties aan te pakken, wanneer
we k niet ontzettend groot willen maken om tijd te besparen. Er is echter geen wetmatigheid
voor welke k de scherpte of projectiefout sturend zijn voor elkaar. Handmatige (visuele) controle
is dus vereist.
8.2.1
Lineaire combinatie van objectieffuncties
Een eerste methode bestaat erin de afzonderlijke objectieffuncties te combineren tot een enkele
objectieffunctie. De nieuwe objectieffunctie kan bijvoorbeeld een lineaire combinatie van de
afzonderlijke objectieffuncties zijn. Formeel:
c(λ, Ω) = λc1 (Ω) + (1 − λ)c2 (Ω), 0 ≤ λ ≤ 1
(8.1)
De parameterverzameling Ω die c(λ, Ω) optimaliseert, zal afhankelijk zijn van de keuze van λ.
De keuze van λ bepaalt in grote mate het succes van de objectieffunctie c. In de literatuur staan
HOOFDSTUK 8. GECOMBINEERDE OBJECTIEFFUNCTIES
54
Figuur 8.2: Pareto-efficiënte oplossingenverzameling.
methodes beschreven om de optimale λ te bepalen bij gegeven c1 en c2 [28]. Deze methodes steunen op minimax (ook toepasbaar bij niet-lineaire combinaties van objectieffuncties) en vereisen
het kunnen afleiden van c1 en c2 . Zoals al eerder aan bod kwam, is het door de aard van het
(k)
(k)
CT-probleem niet mogelijk om EΩ en SΩ analytisch af te leiden. Er is met andere woorden
(k)
(k)
in dit geval geen trefzekere methode voorhanden om EΩ en SΩ betrouwbaar tot een enkele
objectieffunctie te herleiden.
8.2.2
Pareto-efficiënte oplossingenverzameling
Een tweede methode, in de context van genetische algoritmes bekend onder de noemer MOEA
(multiple objective evolutionary algorithm), zoekt naar oplossingen Ω die bij gegeven c1 (Ω) de
hoogste waarde voor c2 (Ω) opleveren, en vice versa (in de veronderstelling dat zowel c1 als c2
gemaximaliseerd moeten worden). Een MOEA moet tegemoetkomen aan optimalisatieproblemen
met twee of meer conflicterende objectieffuncties in een grote en complexe zoekruimte [29]. Het
mag al duidelijk geworden zijn dat het CT-probleem aan deze voorwaarden voldoet.
Concreet is de uitkomst van een MOEA een (benaderde) verzameling van Pareto-efficiënte oplossingen Ω. Die verzameling bevat alle oplossingen Ω die niet door andere oplossingen gedomineerd
worden. Een oplossing Ω1 domineert Ω2 (Ω1 Ω2 ) als aan volgende voorwaarde voldaan is [29]:
Ω1 Ω2 ⇔(c1 (Ω1 ) ≥ c1 (Ω2 ) ∧ c2 (Ω1 ) > c2 (Ω2 ))
∨(c1 (Ω1 ) > c1 (Ω2 ) ∧ c2 (Ω1 ) ≥ c2 (Ω2 ))
(8.2)
Een MOEA kan dus meerdere optima opleveren die elk een verschillende afweging tussen de objectieffuncties representeren. Dit betekent dat we alsnog manueel moeten nagaan welke oplossing
uit de verzameling van Pareto-efficiënte oplossingen de beste is, maar dit zal niet opwegen tegen
het voordeel dat we halen uit het kunnen beperken van het aantal iteraties om met een MOEA
richting het optimum te streven. Samengevat toont figuur 8.2 hoe in de objectiefruimte de
Pareto-efficiënte oplossingen overeenstemmen met parameterverzamelingen in de decisieruimte.
Een MOEA evalueert de oplossingen in de objectiefruimte om de parameterverzamelingen in de
decisieruimte aan te passen. Merk op dat beide objectieffuncties als evenwaardig worden be-
HOOFDSTUK 8. GECOMBINEERDE OBJECTIEFFUNCTIES
55
Algoritme 8.1: SPEA2 procedure
1
2
3
4
5
6
Creëer een initiële populatie T0 met extern archief A0 = ∅. Zet t = 0.
Bepaal met functie F de bruikbaarheid van elk individu in Tt en At .
Kopieer alle niet-gedomineerde individuen in Tt en At naar At+1 .
Indien voldaan aan stopcriterium: At+1 is de verzameling optima. Stop.
Creëer een nieuwe populatie Tt+1 door recombinatie en mutatie van parameters uit At+1 .
Zet t = t + 1 en ga naar stap 2.
schouwd. Bovendien is het mogelijk om ook meer dan twee objectieffuncties te gebruiken, al zal
dat niet in dit onderzoek aan bod komen.
8.3
Optimalisatiealgoritme
Om evolutionair naar de verzameling van Pareto-efficiënte oplossingen te streven, kunnen we
gebruikmaken van het SPEA (strength pareto evolutionary algorithm). Het is een techniek die
in de literatuur als beter wordt beschouwd dan andere MOEA’s [30]. Meer bepaald SPEA2
blijkt goed te presteren; in tegenstelling tot voorganger SPEA convergeert SPEA2 beter naar
een optimum [31]. Het heeft geen belang de verbeteringen van SPEA2 ten opzichte van SPEA
te bespreken. We doen meteen kort de werking van SPEA2 uit de doeken in algoritme 8.1.
Merk op dat de populatie- en archiefgrootte, respectievelijk |Tt | en |At |, constant blijven voor
alle tableaus t. De lus tussen stap 2 en stap 6 is gelijkaardig aan de lus uit het gewone genetisch
algoritme. Het zijn echter de bruikbaarheidsfunctie F in stap 2 en het creëren van een nieuwe
populatie Tt+1 in stap 5 die eigen zijn aan een MOEA in het algemeen en aan SPEA2 in het
bijzonder. Deze stappen moeten we dus definiëren om inzicht te verwerven in de precieze werking
van SPEA2.
8.3.1
Bruikbaarheidsfunctie
De bruikbaarheidsfunctie F : Ω ∈ Rn → R kent aan elk individu Ω in Tt en At een numerieke
waarde toe die aangeeft hoe goed het individu presteert en dus bruikbaar is om tot de verzameling
van Pareto-efficiënte oplossingen te komen. Formeel is F (Ω) de som van de ruwe bruikbaarheid
R(Ω) en de densiteit D(Ω) [31]:
F (Ω) = R(Ω) + D(Ω)
(8.3)
De ruwe bruikbaarheid R(Ω) geeft aan hoeveel oplossingen gedomineerd worden door de oplossingen die Ω domineren [31]:
R(Ω) =
X
Ω0 ∈Tt ∪At ,Ω0 Ω
|{Ω00 |Ω00 ∈ Tt ∪ At ∧ Ω0 Ω00 }|
(8.4)
HOOFDSTUK 8. GECOMBINEERDE OBJECTIEFFUNCTIES
56
| · | duidt hier op de kardinaliteit van een verzameling. R(Ω) = 0 duidt op een niet-gedomineerd
individu. Een grotere waarde voor R(Ω) wijst erop dat Ω door veel andere individuen wordt
gedomineerd die op hun beurt ook veel individuen domineren.
Het is meteen duidelijk dat bij het toekennen van de bruikbaarheid aan een individu de volledige
populatie mee in rekening wordt genomen. De bruikbaarheidsfunctie vervangt het gebruik van
een objectieffunctie. De individuele objectieffuncties zitten in de bruikbaarheidsfunctie vervat
door het gebruik van de operator (vergelijking 8.2) in het bepalen van de ruwe bruikbaarheid
R(Ω) (vergelijking 8.4). Mocht de bruikbaarheidsfunctie F enkel afhangen van R, dan zal F
slecht presteren in het geval er in het archief enkel individuen opduiken die niet gedomineerd
worden, maar zelf ook geen andere individuen domineren. Er is nood aan een extra factor: de
densiteit D. Om de diversiteit van de Pareto-efficiënte oplossingen te bewaren, moet voorrang
gegeven worden aan oplossingen die zich niet in de nabijheid van andere oplossingen bevinden
in de objectiefruimte. Een efficiënte metriek voor D(Ω) bestaat erin de de inverse te nemen van
de afstand δΩ tot de naaste buur van Ω. Hoe groter de densiteit D(Ω) van een individu Ω, hoe
kleiner de kans dat Ω in aanmerking komt om nog langer tot de verzameling van Pareto-efficiënte
oplossingen te behoren. Formeel [31]:
D(Ω) =
1
δΩ + 2
(8.5)
Het toevoegen van 2 in de noemer moet verhinderen dat D(Ω) groter of gelijk wordt aan 1.
R en D bepalen nu volledig F , zoals gedefinieerd in vergelijking 8.3. De complexiteit van het
bepalen van F (Ω) wordt gedomineerd door het bepalen van D(Ω) dat complexiteit O(|Tt ∪
At |2 log |Tt ∪ At |) heeft. De berekening van R bezit complexiteit O(|Tt ∪ At |2 ). Gelet op de
beperkte populatie- en archiefgrootte die we in dit onderzoek zullen gebruiken, en de lange reconstructietijd per te onderzoeken individu, zorgen deze complexiteiten nauwelijks voor beperkingen
in tijd.
8.3.2
Creatie nieuwe populatie
Uit Tt ∪ At worden de individuen geselecteerd waarvoor F (Ω) < 1: dit zijn de individuen die
Pareto-efficiënt zijn met betrekking tot de twee kostfuncties.
Een volledige implementatie is beschikbaar in MOEA.hh, met Genetics.hh als ouderklasse.
8.4
Testresultaten
Vooreerst optimaliseren we met het genetisch algoritme de parameters ∆u, ∆v, θ en SO tegelijkertijd door na 20 iteraties de projectiefout te evalueren. Vervolgens doen we hetzelfde op basis
van de scherpte. In een derde experiment evalueren we zowel de projectiefout als de scherpte in
een MOEA. We zetten |Tt | = 25 en |At | = 12.
HOOFDSTUK 8. GECOMBINEERDE OBJECTIEFFUNCTIES
(a) scherpte
(b) projectiefout
(c) combinatie
57
(d) optimum
Figuur 8.3: Een MOEA levert voor k = 20 de beste schatting.
Geen enkele optimalisatieprocedure levert het optimum met k = 20, maar het MOEA streeft wel
het beste naar de optimale parameters en levert visueel de beste reconstructie. De resultaten
staan samengevat op figuur 8.3.
8.5
Conclusie
Het gebruik van een MOEA is nuttig om bij een beperkt aantal iteraties een idee van het optimum te krijgen. Het optimum wordt nooit bereikt, maar de benadering is goed genoeg om een
bruikbaar beeld op te leveren of om als initiële schatting te gebruiken voor een optimalisatieproces met meer iteraties. In geval van driedimensionale reconstructies zullen we zien hoe een
MOEA kan bijdragen tot het vinden van een afweging tussen maximale scherpte en minimale
projectiefout.
Hoofdstuk 9
Reconstructie in drie dimensies
Driedimensionale reconstructies worden typisch uitgevoerd met de rekenkracht van een GPU. Het
is de bedoeling dat we de algoritmes ontworpen in voorgaande hoofdstukken ultiem testen op een
GPU, met gebruik van ‘echte’ data van kleine proefdieren. Door de relatief lage complexiteit heeft
het heeft weinig zin om de logica achter de optimalisatiealgoritmes naar de GPU te porteren.
De reconstructiealgoritmes daarentegen kunnen we wel nog verder optimaliseren op de GPU om
zoveel mogelijk afbeeldingen tegelijkertijd te reconstrueren. Dit zal het kalibreren aanzienlijk
versnellen. In dit hoofdstuk komen de aanpassingen aan bod om het aantal parallelle taken op
de GPU te verhogen, maar we beginnen met het definiëren van de geometrische parameters die
zich in drie dimensies manifesteren.
9.1
Geometrische parameters
Figuur 9.1 illustreert de gehanteerde coördinatensystemen voor reconstructie in drie dimensies.
Figuur 9.1: Coördinatensysteem voor reconstructie in drie dimensies.
58
HOOFDSTUK 9. RECONSTRUCTIE IN DRIE DIMENSIES
59
Beschrijving
Eenheid
Totale rotatie
Aantal projecties
Dimensie object in x-richting
Dimensie object in y-richting
Dimensie object in z-richting
Dimensie objectvoxel (isotroop)
Dimensie detector in u-richting
Dimensie detector in v-richting
Dimensie detectorpixel (isotroop)
Dimensie brandpunt stralingsbron
◦
—
voxels
voxels
voxels
mm
pixels
pixels
mm
mm
Tabel 9.1: Definitie van de vaste geometrische parameters in drie dimensies.
Symbool
Beschrijving
Eenheid
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
mm
mm
mm
pixels
pixels
◦
mm
mm
Tabel 9.2: Definitie van de te optimaliseren geometrische parameters in drie dimensies.
Tabel 9.1 geeft een overzicht van de geometrische parameters die we in drie dimensies onderscheiden, in het door ons gebruikte model. Het zijn de ‘vaste’ parameters die we niet zullen
kalibreren. Tabel 9.2 bevat de parameters die we zullen optimaliseren naargelang het gedrag van
de objectieffuncties.
Het reconstructiealgoritme haalt de geometrische parameters en de parameters van het genetisch
algoritme uit een tekstgebaseerd configuratiebestand. Bijlage B bevat een voorbeeld van het
gehanteerde formaat. De initiële waarden van de te optimaliseren parameters, alsook de gebruikte
data, zullen we beargumenteren in secties 10.1 en 10.2.
9.2
Specificaties GPU
De gebruikte GPU is een NVIDIA Tesla M2070, aanwezig in een machine met 16 CPU’s van het
type Intel Xeon (kloksnelheid 2,40GHz, 4 kernen per processor). Tabel 9.3 vat de belangrijkste
eigenschappen van de grafische kaart samen. Het zijn deze eigenschappen die het verder parallelli-
HOOFDSTUK 9. RECONSTRUCTIE IN DRIE DIMENSIES
60
GPU
CUDA versie
Hoofdgeheugen
Aantal processors
Aantal kernen per processor
Kloksnelheid processor
Maximale dimensie texture (2D)
Maximale dimensie texture (3D)
Maximaal aantal threads per processor
Maximaal aantal threads per block
Maximale dimensie block
Maximale dimensie grid
CUDA 5.0
5375 MB
14
32
1147 MHz
65536 × 65536
2048 × 2048 × 2048
1536
1024
1024 × 1024 × 64
65536 × 65536 × 65536
Tabel 9.3: Specificaties van de GPU gebruikt tijdens beeldreconstructie (NVIDIA Tesla M2070).
seren van de reconstructiealgoritmes zullen beïnvloeden. NVIDIA heeft het programmeermodel
CUDA ontwikkeld om de virtuele instructieset en het geheugen van de GPU toegankelijk te
maken voor de ontwikkelaar. Met CUDA C/C++ alloceren we vanuit C++ code in de CPUomgeving van de server geheugen op de GPU. Vervolgens stuurt de CPU instructies naar de
GPU om taken te laten uitvoeren op de data in het grafisch geheugen. Alle processorkernen van
de GPU werken parallel ten opzichte van elkaar. Resultaten komen via het grafisch geheugen
terug in het hoofdgeheugen van de CPU-omgeving terecht [33].
Elke processorkern kan een aantal threads parallel uitvoeren. Threads worden vaak samengenomen in zogenaamde blocks. Blocks kunnen we desgewenst in een grid groeperen. Omdat de
termen threads, blocks en grids zo eigen zijn aan de woordenschat van CUDA, zullen we niet
met de Nederlandstalige equivalenten werken, ook al zouden respectievelijk ‘draden’, ‘blokken’
en ‘rasters’ geschikt zijn.
9.3
Voorwaartse en terugwaartse projectie
Een robuuste implementatie van voorwaartse en terugwaartse projectie was reeds voorhanden
in de vorm van de klassen XOCT.hh (CPU) en XOCT.cu (GPU). Deze code zal dus niet meer
in detail besproken worden. Maar het parallellisme dat de GPU kan bieden, werd er ingezet
om iteratief tot één reconstructie te komen. Omdat kalibreren het vergelijken van verschillende
reconstructies inhoudt, moeten we onderzoeken hoe we het parallellisme van de GPU ook kunnen
aanwenden om meerdere reconstructies tegelijkertijd uit te voeren. De verwachte snelheidswinst
is eerder beperkt, omdat één reconstructie al tegen de limieten van geheugen en aantal threads
werkt. Toch moet het zeker mogelijk zijn een snelheid te bekomen die hoger ligt dan het serieel
uitvoeren van reconstructies. We bespreken nu kort de stappen die hiertoe benodigd zijn.
HOOFDSTUK 9. RECONSTRUCTIE IN DRIE DIMENSIES
61
Het uitgangspunt is tegelijkertijd de data in het geheugen houden die de stralingsbron, objecten
en detectors modelleren voor verschillende reconstructieparameters. De aanpassingen zitten vervat in de nieuwe klasse XOCT_calibration.hh met gebruik van kernels (functies uitgevoerd
op de GPU) uit XOCT_calibration.cu. Per iteratie van het reconstructiealgoritme voeren
we de projecties parallel uit.
9.3.1
Stralingsbron
Het model van de stralingsbron bevat enkel parameters en geen object- of detectordata. Gezien
er dus geen data permanent in het grafisch geheugen geladen moeten worden, zijn geen andere
maatregelen nodig dan het tegelijkertijd aanmaken van verschillende stralingsbronnen: één voor
elke verzameling van parameters. In de C++ code betekent dit het bijhouden van een pointer
naar een rij van Tube.hh pointers.
9.3.2
Object
Het model van een te reconstrueren afbeelding bevat niet alleen de eigenschappen van het object,
maar ook de ruwe data die zich in het grafisch geheugen bevinden. Naar analogie met de
stralingsbron trachtten we aanvankelijk een pointer naar een rij van Image_GPU.hh pointers
bij te houden, maar deze piste werd uiteindelijk niet gevolgd. Er stelden zich immers twee
problemen.
Ten eerste laat CUDA (nog?) niet toe een rij van textures aan te maken. Een texture is een
intelligent cachesysteem (ROM) dat onder alle processoren van de GPU gedeeld wordt. Tijdens
voorwaartse projectie wordt een afbeelding om snelheidsredenen in een dergelijke texture geladen.
In het geval van een pointer naar Image_GPU.hh pointers zouden de individuele afbeeldingen
alsnog geconcateneerd moeten worden om ze samen in één texture te laden (de omgekeerde
operatie is ook nodig). Dit brengt overbodig rekenwerk met zich mee.
Ten tweede bleek een pointer naar Image_GPU.hh pointers niet ideaal om tijdens terugwaartse
projectie de individuele afbeeldingen op de GPU aan te roepen. Het zou ook hier extra geheugenoperaties vereisen om dit probleem op te lossen.
De oplossing ligt echter voor de hand: alle afbeeldingen die we parallel wensen te behandelen,
concateneren we bij aanvang al in een ‘superafbeelding’. Deze concatenatie voeren we door in de
z-richting. De aangepaste versie van Image_GPU.hh laat toe een ‘superafbeelding’ aan te maken
op basis van de dimensie van een afzonderlijke afbeelding en het aantal afzonderlijke afbeeldingen
dat we wensen te concateneren. Deze aangepaste versie herleidt zich volkomen tot de originele
versie wanneer we aangeven dat de ‘superafbeelding’ slechts één afbeelding mag bevatten.
Gezien alle afzonderlijke afbeeldingen dezelfde dimensies hebben, kunnen we mits zorgvuldige
indexering op betrouwbare wijze een individuele voxel selecteren.
HOOFDSTUK 9. RECONSTRUCTIE IN DRIE DIMENSIES
62
Figuur 9.2: Verdeling van threads over elke voxel van het object.
Figuur 9.3: Verdeling van threads over elke voxel van geconcateneerde objecten.
Terugwaartse projectie
Bij terugwaartse projectie roepen we parallel voor elke voxel een functie op die de terugwaartse
projectie modelleert. Figuur 9.4 toont hoe de threads oorspronkelijk over het grid van blocks
verdeeld waren. We beschouwen een tweedimensionaal grid waarvan de x- en y-dimensie overeenkomen met de x- en y-dimensie van de afbeelding. Elk block bevat evenveel threads als er
voxels zijn in de z-richting van de afbeelding. Elke thread behandelt dus een afzonderlijke voxel.
Figuur 9.5 toont de doorgevoerde aanpassingen voor parallelle reconstructie. Hetzelfde principe
geldt, al wordt het grid nu in de y-richting uitgebreid volgens het aantal parallelle reconstructies.
9.3.3
Detector
Naar analogie met de objecten, concateneren we ook de detectordata tot één geheel. Dit doen
we in de y-richting. Ook hier geldt dat alle afzonderlijke detectors dezelfde dimensies hebben,
zodat we betrouwbaar individuele pixels kunnen selecteren. De nieuwe versie van de klasse
Detector.hh herleidt zich in het geval van een enkele detector tot de oude versie, zodat achterwaartse compatibiliteit ook hier gegarandeerd is.
HOOFDSTUK 9. RECONSTRUCTIE IN DRIE DIMENSIES
63
Figuur 9.4: Verdeling van threads over elke pixel van de detector.
Figuur 9.5: Verdeling van threads over elke pixel van geconcateneerde detectors.
Voorwaartse projectie
Bij voorwaartse projectie roepen we parallel voor elke detectorpixel een functie op die de voorwaartse projectie modelleert. Figuur 9.4 toont hoe de threads aanvankelijk over het grid van
blocks verdeeld waren. Figuur 9.5 toont de nieuwe situatie. We beschouwen een tweedimensionaal grid waarvan de x-dimensie overeenkomt met de x-dimensie van een grid in het geval van
een afzonderlijke detector. De y-dimensie van het grid is de y-dimensie van een enkel detectorgrid
vermenigvuldigd met het aantal tegelijkertijd te behandelen detectors. Zowel de x- als y-dimensie
zijn veelvouden van 16, omdat elk block 16 × 16 threads telt. Elke thread behandelt met andere
woorden een detectorpixel. Pixels die hierdoor niet meer tot de detector zouden behoren, worden
automatisch genegeerd. Het bleek niet mogelijk het aantal threads 16 × 16 = 256 te verhogen,
de specificaties van tabel 9.3 indachtig.
9.4
Objectieffuncties
De berekening van de objectieffuncties werd ook naar de GPU geporteerd.
HOOFDSTUK 9. RECONSTRUCTIE IN DRIE DIMENSIES
9.4.1
64
Projectiefout
De projectiefout doet beroep op de functie GetTotalAbsoluteValue() in de C++ klasse
BinaryData_GPU.hh met gebruik van de gelijknamige kernel uit BinaryData_GPU.cu.
(k)
Concreet komt de berekening van EΩ op deze code neer:
BinaryData_GPU<f l o a t > e s t i m a t e ;
P r o j e c t i o n s S e t <f l o a t > measured ;
ProjectionsSet_GPU<f l o a t > e r r o r ;
// b e r e k e n i n g van de g e s c h a t t e p r o j e c t i e en o p h a l e n van gemeten p r o j e c t i e
TransferFromCPUtoGPU ( measured , e r r o r ) ;
e r r o r . Subtract ( estimate ) ;
float p r o j e c t i o n E r r o r = e r r o r . GetTotalAbsoluteValue ( ) ;
9.4.2
Scherpte
De functie Sharpness() in de C++ klasse Image_GPU.hh gebruikt de gelijknamige kernel
uit Image_GPU.cu. Het gebruik is eenvoudig:
Image_GPU<f l o a t > image ;
// b e r e k e n i n g van de a f b e e l d i n g
f l o a t s h a r p n e s s = image . S h a r p n e s s ( ) ;
9.5
SIRT
MLTR gebruikt niet rechtstreeks de projectiefout uit vergelijking 5.6, maar zou omwille van
zijn snelheid en precisie geschikt zijn om de scherpte van de reconstructies te bepalen. Toch
bleek het door beperkingen op het geheugen van de GPU niet mogelijk om de door ons gescande
data te reconstrueren (zie ook sectie 10.1). De detectordata is te groot. Het grootste geteste
reconstructievolume met MLTR.cu heeft 256×256×256 objectvoxels, 1184×1120 detectorpixels,
256 projecties en 8 subsets.
OS-SART kon niet verder geparallelliseerd worden, eveneens omwille van beperkingen op het
GPU-geheugen. SIRT daarentegen liet nog ruimte over voor parellellisatie, wat resulteerde in
een speedup van 1,76 per reconstructie.
Het bescheiden succes van deze speedup — gelet op het feit dat we nog steeds van één GPU
gebruikmaken — heeft ons ertoe aangezet verder te gaan met SIRT, ook al omdat we dit algoritme
reeds in de voorstudie goed leren kennen hebben. Convergentie is gegarandeerd.
De geparallelliseerde implementatie is beschikbaar in SIRT_calibration.cu. Er kunnen
maximaal 8 reconstructies tegelijkertijd uitgevoerd worden. De grootte van een populatie —
indien we met een genetisch algoritme werken — dient dus best een veelvoud van 8 te zijn.
Hoofdstuk 10
Kalibreren
Na de voorstudie in twee dimensies en het aanpassen van de GPU-code, pakken we in dit hoofdstuk het kalibreren van CT-scanners aan met driedimensionale reconstructies. We halen tevens
de fenomenen en problemen aan die zich specifiek in de context van driedimensionale beeldreconstructie stellen.
10.1
Data
De gebruikte scanner is de TriFoil Triumph-II X-O CT en kwam al kort aan bod in sectie 2.4
omtrent micro-CT.
Het gebruikte fantoom is de zogenaamde ‘plastimouse’ [34]. Dit is een geplastineerde muis met
behoud van zijn inwendige structuren. De keuze voor een proefdier boven een zuiver geometrisch
object is ingegeven door de grilligere structuur die een proefdier inwendig typisch vertoont. Dit
zal, naar we verwachten, een bijkomende uitdaging zijn voor het optimalisatiealgoritme en de
objectieffuncties. Scherpte werd in de literatuur tenslotte voorgesteld als een vervanging voor
entropie, omdat deze laatste objectieffunctie enkel goed zou presteren bij zuiver geometrische
objecten [13]. De keuze voor een dood dier boven een levend exemplaar was tweeledig. Ten
eerste bleek het op het moment van scannen niet mogelijk om het beoogde, levende proefdier
te verdoven. Ten tweede heeft de plastimouse als voordeel dat ze niet beweegt tijdens het
scanproces. Onscherpte waargenomen op de scans kan dus over het algemeen enkel het gevolg
zijn van slechte reconstructieparameters, net hetgeen we in dit onderzoek willen aanpakken. De
plastimouse zal ons in staat stellen de scherpte beter te evalueren, omdat bewegingsonscherpte
vrijwel uitgesloten is.
Figuur 10.1 verduidelijkt de terminologie van de anatomische vlakverdeling, die later in dit
hoofdstuk gebruikt zal worden om de reconstructies te duiden.
65
HOOFDSTUK 10. KALIBREREN
66
Figuur 10.1: Anatomische vlakverdeling.
10.2
Parameters
ˆ om zijn zoektocht
Het optimalisatiealgoritme heeft nood aan een initiële parameterverzameling Ω
ˆ initialiseren op nul betekent
naar betere parameters te starten. Alle individuele parameters in Ω
niet alleen dat het optimalisatiealgoritme erg lang zal moeten zoeken, het is ook een onrealistische
aanname. Men kan immers al een (erg ruwe) schatting van bepaalde parameters maken, al is het
maar met een eenvoudige meetlat. We achten het interessanter om te kijken of we een dergelijke,
ruwe schatting snel beter kunnen maken dan om ‘uit het niets’ de juiste parameters tevoorschijn
te halen. Voortwerken op een ruwe schatting is een situatie die zich in de praktijk het meeste
zal voordoen.
Om de kwaliteit van het kalibreren te evalueren, is nood aan een goede reconstructie met ideale
Ω. Dit optimum is onbekend, maar we kunnen terugvallen op een schatting van Ω die reeds een
scherp en gedetailleerd tomogram opleverde na het scannen. De doelstelling is dan om met het
kalibreren zo dicht mogelijk in de buurt van de reconstructie met deze parameters te komen, of
zelfs beter te doen, indien mogelijk.
Tabel 10.1 geeft de vaste parameters voor de reconstructie, gevolgd door het waarschijnlijke
ˆ Merk op dat er geen eenduidigheid is
optimum Ω. Tabel 10.2 geeft onze initiële schatting Ω.
wat betreft het aantal beduidende cijfers. In het optimalisatiealgoritme zullen hier ook geen
veronderstellingen over gemaakt worden. Figuren 10.2 en 10.3 tonen voor respectievelijk Ω
ˆ het tomogram, na 100 iteraties (SIRT). Het is duidelijk dat Ω
ˆ onscherpte en vergroting
en Ω
ˆ wel degelijk nuttig maakt.
introduceert, wat het kalibreren vanuit Ω
In dit hoofdstuk zullen we steeds SIRT gebruiken. Tenzij anders vermeld, zijn alle getoonde
afbeeldingen tot stand gekomen na k = 100 iteraties van dit reconstructiealgoritme.
HOOFDSTUK 10. KALIBREREN
Symbool
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
67
Beschrijving
Waarde
Eenheid
Totale rotatie
Aantal projecties
Dimensie object in x-richting
Dimensie object in y-richting
Dimensie object in z-richting
Dimensie objectvoxel (isotroop)
Dimensie detector in u-richting
Dimensie detector in v-richting
Dimensie detectorpixel (isotroop)
Dimensie brandpunt stralingsbron
360
512
256
256
256
0,20
2368
2240
0,05
0,05
◦
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
0,0
-38,39
-155,87
-0,276579
299,64
131,17
—
voxels
voxels
voxels
mm
pixels
pixels
mm
mm
mm
mm
mm
pixels
pixels
◦
mm
mm
Tabel 10.1: De vaste parameters gevolgd door het waarschijnlijke optimum Ω.
Symbool
Beschrijving
Waarde
Eenheid
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
0,0
-30,0
-162,0
0,0
293,0
140,0
mm
mm
mm
pixels
pixels
◦
mm
mm
ˆ
Tabel 10.2: Initiële geometrische parameters Ω.
10.3
Eerste tests
Tijdens deze eerste tests, die geen finale resultaten zullen opleveren, willen we kort nagaan in
hoeverre onze voorstudie met betrekking tot tweedimensionale reconstructies standhoudt in drie
dimensies. De experimenten zullen richtinggevend zijn voor hoe we in sectie 10.5 tot de finale
resultaten zullen komen.
HOOFDSTUK 10. KALIBREREN
(a) transversaal
68
(b) coronaal
(c) sagittaal
Figuur 10.2: Reconstructie volgens Ω (tabel 10.1).
(a) transversaal
(b) coronaal
(c) sagittaal
ˆ (tabel 10.2).
Figuur 10.3: Reconstructie volgens Ω
10.3.1
Scherpte maximaliseren
Uit sectie 8.4 is gebleken dat scherpte richtinggevend kan zijn voor het vinden van de correcte
parameters. Toch presteert scherpte ondermaats in de afwezigheid van betrouwbare parameters die de positie van de detector beschrijven (sectie 7.1). De positie van de detector heeft de
meeste invloed op de reconstructiekwaliteit [27]. Omdat het berekenen van de scherpte randgebaseerd is, kan scherpte zelfs dubbele randen bevorderen, zoals we zagen in sectie 7.1.1. Ook de
vergrotingsfactor beïnvloedt de scherpte.
Desalniettemin gaan we in deze nieuwe situatie na of we de geometrische parameters, op de
uitwijkingen van het object ∆x, ∆y en ∆z na, aan de hand van scherpte kunnen optimaliseren.
We laten deze uitwijkingen buiten beschouwing, omdat we uit de kennis van Ω weten dat er reeds
een optimum bestaat met deze waarden op 0. We evalueren de scherpte na 20 iteraties van het
reconstructiealgoritme.
Een tussentijds optimum, dat bekomen werd nog voor het optimalisatiealgoritme zijn stopvoorwaarde bereikt had, toonde al dat de parameters duidelijk in de verkeerde richting evolueerden.
HOOFDSTUK 10. KALIBREREN
(a) transversaal
69
(b) coronaal
(c) sagittaal
Figuur 10.4: Reconstructie die scherper geëvalueerd wordt dan reconstructie met Ω.
Figuur 10.5: Scherpteverloop voor de suboptimale reconstructie uit figuur 10.4.
Het teleurstellende resultaat is te zien op figuur 10.4. We zien de verwachte fenomenen opnieuw
opduiken. Ten eerste is de vergrotingsfactor toegenomen, waarvan we inderdaad al wisten dat
het de numerieke scherpte positief beïnvloedt [13]. Ten tweede is de positie van de detector
zo slecht geschat, dat de geïntroduceerde artefacten (dubbele randen) foutief beschouwd worden als een teken van toegenomen scherpte. Figuur 10.5 toont wel dat het voordeel van deze
reconstructie ten opzichte van Ω steeds kleiner wordt naarmate het aantal iteraties van het reconstructiealgoritme toeneemt. Na 79 iteraties zou Ω scherper geëvalueerd worden. Maar meer
dan 100 iteraties wachten tot scherpte een betrouwbaar resultaat oplevert, is lang. Het verschil
in numerieke scherpte tussen een optimaal en slecht tomogram zal bovendien zo klein zijn, dat
lokale maxima niet uit te sluiten zijn.
Een te sterk toegenomen vergrotingsfactor compenseert het negatieve effect dat slechte detectorparameters op de scherpte hebben. Althans, dat is het resultaat als we scherpte evalueren
na te weinig iteraties van het iteratief reconstructiealgoritme: scherpte convergeert slecht, zoals
HOOFDSTUK 10. KALIBREREN
(a) transversaal
70
(b) coronaal
(c) sagittaal
Figuur 10.6: Intermediair resultaat tijdens het minimaliseren van de projectiefout.
(a) transversaal
(b) coronaal
(c) sagittaal
Figuur 10.7: Intermediair resultaat voor het combineren van scherpte en projectiefout.
getoond op figuur 10.5. Scherpte toepassen op de detectorparameters alleen zal nog mogelijk zijn,
als we met andere woorden de vergroting constant houden. Willen we scherpte toch toepassen
op het optimaliseren van SD en SO, bijvoorbeeld samen met de projectiefout in een MOEA, dan
moeten we de slechte convergentie-eigenschappen van scherpte in het achterhoofd houden.
10.3.2
Projectiefout minimaliseren
Ook in deze test optimaliseren we alle geometrische parameters, op ∆x, ∆y en ∆z na. We
evalueren de projectiefout na 20 iteraties.
∆u streeft als enige parameter het snelste naar zijn optimale waarde. De uitwijking ∆v daarentegen en de afstanden SD en SO wijken sterk af van hun optimale waarden. Vooral in het geval
van ∆v viel de afwijking met het optimum uit Ω op. Nochtans hadden we, uit de voorstudie, het
vermoeden dat de projectiefout minimaliseren een gunstig effect zou hebben op de optimalisatie
van alle verschuivingen en rotaties van de detector in zijn vlak. In sectie 10.4.1 zullen we de
oorzaak toelichten.
HOOFDSTUK 10. KALIBREREN
71
Figuur 10.8: Geprojecteerde data afkomstig van de scanner.
Figuur 10.6 toont een tussentijds optimum tijdens het minimaliseren van de projectiefout. Door
de grootte van de zoekruimte is dit intermediair resultaat niet perfect, en bovendien weten we
uit sectie 6.1 dat het optimaliseren van SD en SO sowieso meer dan 20 iteraties vergt.
10.3.3
Gecombineerde objectieffuncties
Het gebruik van een MOEA op ∆u, ∆v, φ, SD en SO lijkt vooralsnog veelbelovender. Het
algoritme is in staat een afweging te vinden tussen een resultaat met grote scherpte (maar met
mogelijks artefacten, dus een grote projectiefout) en een resultaat met kleine projectiefout, maar
daarom nog niet de beste scherpte. Figuur 10.7 toont een tussentijds optimum door na 20
iteraties scherpte en projectiefout samen te evalueren. Scherpte blijkt dus inderdaad nuttig in
het bijsturen van parameters, maar de effecten van trage convergentie uit figuur 10.5 spelen
nog een te grote rol om met gecombineerde kostfuncties na weinig iteraties tot een optimum te
komen.
Uit deze drie experimenten leiden we af dat ook in drie dimensies een stapsgewijze optimalisatie
aangewezen is, om het aantal iteraties van het iteratieve reconstructiealgoritme beperkt te houden
en per optimalisatiestap de zoekruimte niet te complex te maken. Dat deze methode zal werken,
tonen we aan in sectie 10.5.
10.4
10.4.1
Invloed parameters in beeld- en projectieruimte
Uitwijkingen van het object
In sectie 10.3.2 merkten we op dat ∆v sterk afweek van zijn optimale waarde, terwijl we uit
de voorstudie weten dat de projectiefout in staat is om bij weinig iteraties van het iteratief
reconstructiealgoritme de detectorpositie al nauwkeurig te schatten, onafhankelijk van de andere
parameters. De oorzaak ligt hier in het buiten beschouwing laten van ∆z.
HOOFDSTUK 10. KALIBREREN
72
Op figuur 10.8 tonen we een voorbeeld van de geprojecteerde data, afkomstig van de scanner.
(15)
Kijken we op figuur 10.9 naar de projectie HˆfΩ , dan zien we links artefacten opduiken, wat de
projectiefout vergroot. Dit fenomeen doet zich voor op de plaatsen waar het object tot tegen de
rand van de detector geprojecteerd wordt. In dit geval raakt de projectie de rand van de detector
in zijn v-richting.
Figuur 10.10 toont dat reconstructie volgens Ω, maar met ∆v = −115 pixels, dit probleem in de
projectieruimte lijkt op te lossen. De projectiefout is er kleiner dan de projectiefout voor Ω. De
reconstructie met ∆v = −115 pixels is echter onscherp. Deze lagere kwaliteit in de beeldruimte
kan dus wel al door de numerieke scherpte voorspeld worden.
(15)
Op figuur 10.11 tonen we de projectie HˆfΩ , maar met ∆z = −5 mm. Ook deze ingreep
lost het vastgestelde probleem in de projectieruimte op (en zelfs beter), en de kwaliteit in de
beeldruimte lijkt behouden: zie figuur 10.12. Merk wel dat een numeriek vergelijk in scherpte
met de reconstructie uit figuur 10.2 moeilijk is, door de extra beelddata die nu zichtbaar geworden
is (beide reconstructies bevatten niet evenveel beelddata).
We moeten besluiten dat ∆v en ∆z hetzelfde probleem in de projectieruimte trachten op te
lossen, maar pas door de scherpte van de bekomen reconstructies te vergelijken, kunnen we
bepalen welke parameter het beste aangewezen is om de projectiefout te verkleinen.
Mocht de projectie van het object zich tegen de boven- of onderrand van de detector bevonden
hebben (zijn u-richting), dan zouden aanpassingen van ∆x en ∆u afgewogen moeten worden.
Dit is hier niet het geval en zoals uit het kalibratieproces zal blijken (sectie 10.5), is het optimaliseren van ∆x hier te verwaarlozen. Ook ∆y laten we buiten beschouwing, omdat de andere
vergrotingsparameters SD en SO al erg moeilijk te bepalen zullen zijn.
Dit alles leidt tot uitermate belangrijke conclusies voor het verder verloop van dit onderzoek. Een
optimum (scherp tomogram) kan bestaan bij een projectiefout die niet noodzakelijk minimaal
is. Het aanpassen van ∆z zal immers meer data in beeld brengen, een ongeveer even scherp
tomogram als reconstructie met Ω opleveren, en tot een kleinere projectiefout leiden. Dat een
optimum bekomen kan worden bij een projectiefout die niet noodzakelijk minimaal is, pleit in
het voordeel van een objectieffunctie die in de beeldruimte werkt. Scherpte bleek echter niet in
staat om SD en SO te optimaliseren. Scherpte is bovendien geen betrouwbare objectieffunctie
wanneer de hoeveelheid beelddata wijzigt, door hier bijvoorbeeld ∆v en ∆z aan te passen. De
projectiefout blijft dus onontkoombaar, al was het maar om SD en SO te bepalen.
10.4.2
Detectorparameters
We trachten nu visueel te verklaren waarom de detectorparameters zich snel laten schatten op
basis van projectiefout en scherpte, zoals we numeriek aantoonden in de voorstudie. In de vorige
sectie lichtten we reeds de invloed van ∆v toe, met betrekking tot ∆z. In deze sectie zullen
we kort de invloed van ∆u en φ op de geprojecteerde data bespreken. Figuur 10.13a toont
HOOFDSTUK 10. KALIBREREN
(a) projectie
73
(b) projectiefout
(c) reconstructie
Figuur 10.9: Projectie volgens Ω (k = 15), reconstructie met k = 100.
(a) projectie
(b) projectiefout
(c) reconstructie
Figuur 10.10: Projectie volgens Ω (k = 15, ∆v = 115 pixels), reconstructie met k = 100.
(a) projectie
(b) projectiefout
(c) reconstructie
Figuur 10.11: Projectie volgens Ω (k = 15, ∆z = −5 mm), reconstructie met k = 100.
HOOFDSTUK 10. KALIBREREN
(a) transversaal
74
(b) coronaal
(c) sagittaal
Figuur 10.12: Reconstructie volgens Ω, met ∆z = −5 mm.
(15)
de geprojecteerde data HˆfΩ . Figuur 10.13b toont eveneens projectie volgens Ω, maar met
∆u = −78 pixels. Figuur 10.13c past in Ω de hoek φ = −4◦ aan.
Het is meteen duidelijk waarom deze parameters zich zo snel laten schatten op basis van de
projectiefout. Het gaat hier immers om verschuivingen van de detector in zijn vlak, waardoor er
meteen een fout ten opzichte van de data in figuur 10.8 optreedt: de contouren van de plastimouse
komen niet meer overeen. Figuur 10.14 toont reconstructies voor k = 15. Scherpte blijft, net
als in de voorstudie, bruikbaar om deze parameters na een beperkt aantal iteraties te schatten,
zolang we de vergrotingsfactor constant houden (cf. figuur 10.4).
10.4.3
Vergrotingsfactor
Op figuur 10.15 merken we dat het aanpassen van SD en SO zichtbaar weinig invloed heeft op
de geprojecteerde data, terwijl in de beeldruimte de vergrotingsfactor meteen zichtbaar is.
Dit maakt duidelijk waarom we de detectorparameters ∆u, ∆v en φ afzonderlijk kunnen optimaliseren. Hun effect op de geprojecteerde data is groter, omdat er een verplaatsing in het vlak
plaats vindt. Deze verplaatsing resulteert in een toenemende projectiefout en een afgenomen
scherpte. De parameters die de vergroting beïnvloeden, SD en SO, tonen op de geprojecteerde
data enkel wazige randen rond de contouren van de plastimouse. Deze wazige randen hebben pas
na veel iteraties een doorslaggevende invloed op de projectiefout. Het aligneren van de contouren
van Hˆf (k) met g heeft een grotere impact dan de randonscherpte op Hˆf (k) die zich voordoet bij
het aanpassen van SD en SO.
We merken tevens op dat SD en SO de artefacten op de geprojecteerde data beïnvloeden. Als
we de projectiefout willen gebruiken om SD en SO te optimaliseren, zal ∆z hoe dan ook geoptimaliseerd moeten worden om SD en SO niet in de verkeerde richting te duwen.
We kunnen nu overgaan tot het formuleren van een stappenplan om de CT-scanner te kalibreren.
HOOFDSTUK 10. KALIBREREN
(a) Ω
75
(b) Ω, ∆u = −78 pixels
(c) Ω, φ = −4◦
Figuur 10.13: Projecties voor k = 15.
(a) Ω
(b) Ω, ∆u = −78 pixels
(c) Ω, φ = −4◦
Figuur 10.14: Reconstructies voor k = 15.
(a) Ω
(b) Ω, SD = 310 mm
Figuur 10.15: Projecties voor k = 15.
(c) Ω, SO = 141 mm
HOOFDSTUK 10. KALIBREREN
10.5
76
Stappenplan
Gebaseerd op de voorstudie, waar we parameters afzonderlijk konden optimaliseren, bouwen we
ook hier een stappenplan om te kalibreren uit, mede gebaseerd op de kennis uit de experimten
van sectie 10.3 en 10.4.
10.5.1
Optimaliseren van ∆u en φ
∆u, ∆v en φ hebben een grote invloed op de reconstructiekwaliteit [27]. Uit sectie 10.4.1 weten
we dat ∆v en ∆z dezelfde fout kunnen compenseren in de projectieruimte, maar dat hun invloed
in de beeldruimte anders is. Vier parameters optimaliseren — ∆u, ∆v, φ en ∆z — maakt de
zoekruimte al aardig complex. Bovendien vereist het optimaliseren van ∆v en ∆z een afweging
tussen minimale projectiefout en maximale scherpte (cf. figuren 10.10 en 10.11), wat in de
richting van een MOEA wijst. We kiezen ervoor om eerst ∆u en φ te optimaliseren, op basis
van scherpte. Tabel 10.3 vat de eerste optimalisatiestap samen. De evaluatie van scherpte
gebeurt na 2 iteraties van het iteratief reconstructiealgoritme. We zetten een bovengrens op de
projectiefout. De projectiefout zelf minimaliseren, zou pas betrouwbare resultaten opleveren na
een viertal iteraties.
Tabel 10.4 geeft de resultaten van het optimalisatiealgoritme, voor drie verschillende tests. Figuur
ˆ 1,1 .
10.16 toont het tussentijdse resultaat met de parameterverzameling Ω
Symbool
Beschrijving
Waarde
Eenheid
Wijzigen
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
0,0
-30
162
0,0
293,0
140,0
mm
mm
mm
pixels
pixels
nee
nee
nee
ja
nee
ja
nee
nee
◦
mm
mm
Tabel 10.3: Startwaarden genetisch algoritme op scherpte, met |Tt | = 12.
ˆ 1,1
Ω
ˆ 1,2
Ω
ˆ 1,3
Ω
∆u
φ
-37,8105 pixels
-37,7601 pixels
-37,9912 pixels
-0,2621◦
-0,2329◦
-0,2765◦
Tabel 10.4: Behaalde resultaten na de eerste kalibratiestap.
HOOFDSTUK 10. KALIBREREN
(a) transversaal
77
(b) coronaal
(c) sagittaal
ˆ 1,1 (tabel 10.4).
Figuur 10.16: Resultaat stap 1: reconstructie volgens Ω
10.5.2
Optimaliseren van ∆v en ∆z
Wanneer het object tijdens het scannen tot tegen de rand van de detector geprojecteerd werd in
de v-richting (de u-richting), dan kunnen ∆v en ∆z (∆u en ∆y) afzonderlijk de fout compenseren
die mogelijks ontstaat in de projectieruimte (cf. figuur 10.9). Het resultaat in de beeldruimte
is evenwel anders, cf. figuren 10.10 en 10.11. Scherpte kan hier uitsluitsel geven over de beste
reconstructie.
∆v en ∆z enkel op scherpte optimaliseren is evenwel niet ideaal. Het aanpassen van deze parameters zal meer of minder data in beeld brengen, wat een vertekend beeld oplevert van de
numerieke scherpte. We kiezen om met een MOEA snel een afweging te vinden tussen oplossingen met kleine projectiefout en grote scherpte, en kleine projectiefout en lagere scherpte. Het
gebruik van een MOEA heeft als voordeel dat de projectiefout de richting van de parameters
∆v en ∆z zal kunnen bepalen, en dat scherpte bepaalt welke parameter best niet te veel de
projectiefout compenseert om geen slechte resultaten in de beeldruimte te bekomen.
Met de oplossingen die het MOEA na een beperkt aantal tableaus oplevert, kunnen we een
genetisch algoritme starten dat enkel de projectiefout minimaliseert, en zich houdt aan de richting
waarin de parameters tijdens de uitvoering van het MOEA evolueerden. Het aantal iteraties
waarna we de objectieffuncties evalueren is eveneens beperkt: k = 4 geeft goede resultaten.
Tabel 10.5 geeft een overzicht van deze kalibratiestap.
Tabel 10.6 toont de bekomen waarden voor twee verschillende tests. Figuur 10.17 toont het
ˆ 2,1 .
tussentijdse resultaat met de parameterverzameling Ω
10.5.3
Optimaliseren van SD en SO
Enkel door het minimaliseren van de projectiefout kunnen we de parameters SD en SO bepalen. Het grote nadeel is dat dit erg veel iteraties van het iteratief reconstructiealgoritme vergt.
HOOFDSTUK 10. KALIBREREN
78
Symbool
Beschrijving
Waarde
Eenheid
Wijzigen
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
0,0
-37,8105
-162
-0,2621
293,0
140,0
mm
mm
mm
pixels
pixels
nee
nee
ja
nee
ja
nee
nee
nee
◦
mm
mm
Tabel 10.5: Startwaarden genetisch algoritme op projectiefout, met |Tt | = 12.
ˆ 2,1
Ω
ˆ 2,2
Ω
∆v
∆z
-151,993 pixels
-150,788 pixels
-4,22209 mm
-4,38001 mm
Tabel 10.6: Behaalde resultaten na de tweede kalibratiestap.
(a) transversaal
(b) coronaal
(c) sagittaal
ˆ 2,1 (tabel 10.6).
Figuur 10.17: Resultaat stap 2: reconstructie volgens Ω
Omdat we op figuur 10.17 al een scherp tomogram waarnemen, maken we een afweging tussen
nauwkeurigheid, bruikbaarheid en tijd. We evalueren de projectiefout na k = 80.
Merk op dat het kalibreren van ∆z strikt noodzakelijk was. Anders zouden SD en SO de fout
in de projectieruimte compenseren die ∆z (in samenwerking met ∆v) kan oplossen, cf. figuur
10.15.
Tabel 10.7 toont de startwaarden voor deze derde kalibratiestap, met in tabel 10.8 de behaalde
resultaten.
HOOFDSTUK 10. KALIBREREN
79
Symbool
Beschrijving
Waarde
Eenheid
Wijzigen
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
-4,22209
-37,8105
-151,993
-0,2621
293,0
140,0
mm
mm
mm
pixels
pixels
nee
nee
nee
nee
nee
nee
ja
ja
◦
mm
mm
Tabel 10.7: Startwaarden genetisch algoritme op projectiefout, met |Tt | = 8.
ˆ 3,1
Ω
SD
SO
305,6862 mm
127,6113 mm
Tabel 10.8: Behaalde resultaten na de derde kalibratiestap.
(a) transversaal
(b) coronaal
(c) sagittaal
ˆ 3,1 (tabel 10.8).
Figuur 10.18: Resultaat stap 3: reconstructie volgens Ω
10.6
Conclusie
Figuur 10.18 toont het finale resultaat na kalibreren. Tabel 10.9 vat de geoptimaliseerde parameters nog eens samen.
10.6.1
Herhaalbaarheid
We hebben elke kalibratiestap een aantal keer herhaald, op het optimaliseren van SO en SD in
stap 3 na, uit tijdbeperkingen. Het genetisch algoritme blijkt inderdaad goed naar de optima te
streven en de resultaten liggen erg dicht in elkaars buurt.
HOOFDSTUK 10. KALIBREREN
80
Symbool
Beschrijving
Waarde
Eenheid
∆x
∆y
∆z
∆u
∆v
φ
SD
SO
Uitwijking object in x-richting
Uitwijking object in y-richting
Uitwijking object in z-richting
Uitwijking detector in u-richting
Uitwijking detector in v-richting
Hoek detector (twist)
Afstand stralingsbron tot detector
Afstand stralingsbron tot object
0,0
0,0
-4,22209
-37,8105
-151,993
-0,2621
305,6862
127,6113
mm
mm
mm
pixels
pixels
◦
mm
mm
Tabel 10.9: Eindresultaat van het kalibreren.
10.6.2
Nauwkeurigheid
Gezien het gegeven optimum Ω slechts een benadering was, is het moeilijk de geoptimaliseerde
parameters ermee te vergelijken. We beperken ons tot een visuele controle, en dan merken we
dat de artefacten ontstaan door foute detectorparameters opgelost werden. De nauwkeurigheid
van de resultaten voor SO en SD is moeilijk in te schatten. Het optimalisatiealgoritme heeft
geen stopvoorwaarde bereikt. Toch konden we met het behaalde resultaat een visueel goede
reconstructie bekomen.
Als het genetisch algoritme voor ∆u de waarde -37,8105 teruggeeft, betekent dit niet noodzakelijk
dat -37,8106 of -37,8104 slechtere resultaten opleverden. De nauwkeurigheid, in termen van
beduidende cijfers, is afhankelijk van de precisie (float) en de standaardafwijking gebruikt om
mutatie uit te voeren tijdens de uitvoering van het genetisch algoritme. Deze standaardafwijking
kan weliswaar kleiner gemaakt worden in de tijd om de resultaten te verfijnen, maar dan bestaat
het risico dat het genetisch algoritme niet meer uit lokale extrema ontsnapt.
10.6.3
Tijd
Het duurt op de GPU ongeveer 3 minuten en 26 seconden om SIRT 1 iteratie van 1 reconstructie
te laten uitvoeren. Dit betekent dat het berekenen van één tableau van het genetisch algoritme
met 12 individuen per iteratie 42 minuten duurt.
10.6.4
Samenvatting
Figuur 10.19 toont het te volgen stappenplan bij kalibreren.
Het optimaliseren van ∆u en φ gebeurt best op scherpte en de evaluatie kan na een beperkt
aantal iteraties: k = 2 geeft al goede resultaten. Een bovengrens op de projectiefout vermijdt
dat de scherpte naar andere maxima streeft.
HOOFDSTUK 10. KALIBREREN
81
Figuur 10.19: Stappenplan voor het kalibreren van CT-scanners.
Vervolgens kalibreren we ∆v en ∆z. Wanneer het gescande object tot tegen de rand van de
detector geprojecteerd werd in de v-richting (de u-richting), dan zullen ∆v en ∆z (∆u en ∆y)
de fout compenseren die mogelijks ontstaat in de projectieruimte (cf. figuur 10.9). Scherpte
moet uitsluitsel geven welke parametercombinatie in de beeldruimte het scherpste beeld oplevert. Om de afweging tussen grote scherpte en kleine projectiefout te maken, hanteren we eerst
zowel scherpte als projectiefout. Wanneer we de richting van de parameters bepaald hebben,
minimaliseren we de projectiefout. Het aantal iteraties is eveneens beperkt, k = 4 geeft goede
resultaten.
Een optionele stap is om ∆u, ∆v en φ nog eens te optimaliseren, op basis van scherpte. De
mogelijkheid bestaat immers dat ∆u en φ moeten bijgestuurd worden door de nieuwe waarde
van ∆v. In de praktijk blijkt deze invloed echter beperkt.
De laatste stap bestaat erin SO en SD te optimaliseren volgens de projectiefout. Het gebruik
van scherpte heeft hier geen nut. Het aantal iteraties moet wel groot gemaakt worden: de
convergentie van de projectiefout is slecht voor deze parameters.
Het blijkt dat optima kunnen bestaan bij een projectiefout die niet minimaal is. Dit doet ons
denken in de richting van een objectieffunctie die enkel in de beeldruimte werkt, maar scherpte is
onbruikbaar bij de parameters die de vergrotingsfactor beïnvloeden. Met de vraagtekens die men
bij de minimalisatie van entropie moet plaatsen [13], bestaat er dus geen algemeen toepasbare
methode in de beeldruimte.
Het gebruik van de projectiefout blijft onontkoombaar. Gelukkig hebben we kunnen aantonen
dat parameters apart gekalibreerd kunnen worden, in aanwezigheid van andere parameters die
mogelijks foute waarden bezitten. Op die manier daalt de complexiteit van de zoekruimte aanzienlijk. Uit tijdbeperkingen is het niet mogelijk om een statistische analyse te maken van de
geoptimaliseerde parameters, zoals bijvoorbeeld een standaardafwijking bepalen op de bekomen
optima. Hoewel strikt genomen niet statistisch relevant, hebben we door een paar keer de experimenten te herhalen toch gezien dat min of meer dezelfde optima worden bekomen. Het
stappenplan uit figuur 10.19 blijkt dus bruikbaar en de optimalisatiealgoritmes presteren goed.
Hoofdstuk 11
Conclusie
11.1
Resultaten
Het minimaliseren van de projectiefout blijkt erg robuust te zijn om alle parameters tegelijkertijd te optimaliseren, zolang het aantal iteraties van het iteratief reconstructiealgoritme maar
voldoende groot gemaakt wordt. In de praktijk betekent dit een al te complexe zoekruimte en
een optimalisatie die te veel tijd in beslag neemt.
Scherpte als objectieffunctie levert geen betere resultaten op zonder ook de projectiefout mee
in rekening te brengen. Deze objectieffunctie is immers ongeschikt om de parameters die de
vergroting beïnvloeden, te optimaliseren. Scherpte heeft wel een richtinggevend karakter: bij
weinig iteraties kan het al de richting aanduiden waarin de parameters moeten evolueren opdat
ze hun optimale waarden zouden bereiken. We toonden dit aan door het gebruik van een MOEA.
We hebben gebruik kunnen maken van de veelbelovende resultaten uit het werk van Kingston
et al. [13] wat betreft scherpte. Toch hebben we situaties kunnen ontdekken waarin scherpte
geen globaal maximum vertoont bij het optimum. Dubbele randen, een typisch artefact bij het
fout inschatten van de detectorparameters, kunnen zelfs door scherpte bevorderd worden. Een
bovengrens op de projectiefout blijkt nodig. Scherpte optimaliseert de detectorparameters wel
iets sneller dan de projectiefout zelf.
Het werk van Wein et al. [15] steunt volledig op het gebruik van de projectiefout. Zij optimaliseren alle parameters tegelijkertijd aan de hand van het Nelder-Mead algoritme. Wij hebben
aangetoond dat dit voor een correcte inschatting veel iteraties vereist, en dat de Nelder-Mead
procedure ongeschikt is wegens de aanwezigheid van te veel lokale extrema. Dit verklaart de
grote afwijkingen tussen de schatting en het optimum in hun werk, ook al konden ze tot een
schatting komen die visueel een goede reconstructie opleverde.
Wij hebben een methode uitgewerkt om de parameters stapsgewijs te optimaliseren. Dit verkleint
de zoekruimte aanzienlijk en bovendien kunnen sommige parameters al na een erg klein aantal
82
HOOFDSTUK 11. CONCLUSIE
83
iteraties geschat worden.
Sawall et al. [3] stellen terecht een genetisch algoritme voor om te kunnen ontsnappen aan
de lokale extrema waarin Nelder-Mead typisch verzeild raakt. Toch optimaliseren ook zij alle
parameters samen, wat de zoekruimte onnodig complex maakt.
Ondanks het gebrek aan de systeemmatrix, waarmee Bronnikov [18] een algemene optimalisatiemethode uitwerkte, hebben we dus toch alle parameters kunnen optimaliseren en hun onderlinge
verbanden in de beeld- en projectieruimte kunnen vaststellen.
11.2
Toekomstig werk
Toekomstig werk moet zich toeleggen op het gebruik van een alternatief voor SIRT dat zich
eveneens goed laat parallelliseren. Kalibreren is immers een tijdrovend proces. Voorts hebben
we het kalibratieprobleem hier in al zijn algemeenheid aangepakt. Het kan nuttig zijn om na te
gaan of het bijvoorbeeld mogelijk is dezelfde principes toe te passen op geschaalde projector- en
reconstructiedata.
We hebben niet kunnen testen of de scherpte ook bruikbare resultaten levert wanneer het object
beweegt tijdens de acquisitie. Een mogelijkheid bestaat erin na te gaan of we het kalibratieproces
ook kunnen toepassen door de beeldruimte te beperken tot het bed waarop het object rust, in
de veronderstelling dat dit stabiel blijft.
Bijlage A
Inhoud digitale drager
Aan deze thesis werd een digitale drager toegevoegd met daarop de (geannoteerde) code van
reconstructie- en optimalisatiealgoritmes. Ook enkele reconstructies van de plastimouse die representatief zijn voor dit onderzoek zijn aanwezig.
A.1
CPU
CPU/code/src/
Reconstructie op de CPU gebruikt C++11 en OpenMP. Compileren vereist volgende syntax:
g++ −fopenmp −s t d=c++11 −O3 −p i p e f o o . cpp −o . . / b i n / f o o
Voorwaartse projectie:
• ForwardProject2D.cpp
• ForwardProject2D_geometry.cpp
Reconstructie:
• SIRT.cpp
• SIRT_geometry.cpp
Kalibreren:
•
•
•
•
•
SIRT_genetics_error.cpp
SIRT_genetics_sharpness.cpp
SIRT_neldermead_error.cpp
SIRT_neldermead_sharpness.cpp
SIRT_moea_error_sharpness.cpp
84
BIJLAGE A. INHOUD DIGITALE DRAGER
CPU/code/include/calibration/
Dataopslag:
• DataTable.hh
• DataLog.hh
Optimalisatie:
•
•
•
•
Simplex.hh
Genetics.hh
Population.hh
MOEA.hh
CPU/data/
Fantomen:
•
•
•
•
phantom.png
phantom.raw
phantom_small.png
phantom_small.raw
Gesimuleerde projectie:
• sinogram_small.raw
Configuratiebestand:
• calibration.conf
A.2
GPU
GPU/code/CTIRGPU/include/calibration/
De voorzieningen voor dataopslag en optimalisatie zijn identiek aan de CPU-versies.
85
BIJLAGE A. INHOUD DIGITALE DRAGER
86
GPU/code/CTIRGPU/include/recon/XOCT/
Reconstructie:
•
•
•
•
•
Detector.hh
Detector.cu
Tube.hh
XOCT_calibration.hh
XOCT_calibration.cu
In de ImagingLibrary werden BinaryData_GPU, BinaryData3D_GPU en Image_GPU
aangepast.
GPU/code/CTIRGPU/include/src/
Kalibreren:
• SIRT_calibration.cu
Compileren kan met de bijgevoegde Makefile.
GPU/data/
Configuratiebestand:
• plastimouse.txt
Reconstructies:
• given_SIRT_100.img
• guess_SIRT_100.img
• result_SIRT_100.img
Bijlage B
Configuratiebestanden
B.1
Principe
Bij het kalibreren van de geometrische parameters heeft het reconstructiealgoritme nood aan
een tekstgebaseerd configuratiebestand. Dit configuratiebestand bevat de eigenschappen van de
CT-scanner en de geometrische parameters die het model gebruikt. De geometrische parameters vormen de initiële schatting waaruit het optimalisatiealgoritme vertrekt. Per individuele
parameter kan met 0 of 1 in het binaire veld calibrate aangeduid worden of deze parameter
aangepast en dus gekalibreerd mag worden.
Het configuratiebestand gaat uit van een genetisch algoritme. Het binaire veld sharpness
laat toe om het reconstructiealgoritme de scherpte te laten maximaliseren. Analoog laat het
veld projection_error toe om de projectiefout te minimaliseren. Beide velden op 1 zetten,
betekent dat het reconstructiealgoritme automatisch een MOEA zal gebruiken. Populatie- en
archiefgrootte zijn in te stellen in de respectievelijke velden population en survivors.
De betekenis van de andere velden zou vanzelfsprekend moeten zijn.
B.2
Formaat
[ CTIR ]
access_order = sequential
i n p u t = /home/ jim / data /001 _ p l a s t i m o u s e /20140326 T144248 /Raw.%04 i
output = /home/ jim / workspace / p l a s t i m o u s e _ o u t p u t / r e c o n
s t o r a g e = /home/ jim / workspace / p l a s t i m o u s e _ o u t p u t /
rays = 1
subsets = 4
airraw_path = /home/ jim / data /001 _ p l a s t i m o u s e /20140326 T144248 / a i r r a w
cudaDevice = 0
spectrum = /home/ i n f i n i t y / workspace / CTspectra / poly_50kVp_triumph . s p c
87
BIJLAGE B. CONFIGURATIEBESTANDEN
subsample = 1
[ Store ]
normalization = 0
backNormalization = 0
estimated = 0
projection_error = 0
image_error = 0
reconstructed = 1
[ Geometry ]
a n g l e F i l e = /home/ jim / data /001 _ p l a s t i m o u s e /20140326 T144248 / a n g l e . b i n
n u m b e r _ p r o j e c t i o n s = 512
skip_projections = 1
t o t a l _ r o t a t i o n = 360
detectorDimY = 2368
detectorDimZ = 2240
p i x e l s i z e = 0.05
imageDimX = 256
imageDimY = 256
imageDimZ = 256
voxelsize = 0.2
f o c a l s p o t s i z e = 0.05
[ Calibration ]
sharpness = 1
projection_error = 0
p o p u l a t i o n = 40
survivors = 2
i t e r a t i o n s = 10
offset = 8
0 = ImageOffsetX
1 = ImageOffsetY
2 = ImageOffsetZ
3 = DistanceDetectorFocalSpot
4 = DistanceObjectFocalSpot
5 = DetectorOffsetU
6 = DetectorOffsetV
7 = Pivoting
[ ImageOffsetX ]
calibrate = 0
value = 0
stdev = 1
[ ImageOffsetY ]
calibrate = 0
value = 0
stdev = 1
88
BIJLAGE B. CONFIGURATIEBESTANDEN
[ ImageOffsetZ ]
calibrate = 0
value = 0
stdev = 1
[ DistanceDetectorFocalSpot ]
calibrate = 1
v a l u e = 293
stdev = 1
[ DistanceObjectFocalSpot ]
calibrate = 1
v a l u e = 140
stdev = 1
[ DetectorOffsetU ]
calibrate = 1
v a l u e = −30
stdev = 1
[ DetectorOffsetV ]
calibrate = 1
v a l u e = −162
stdev = 1
[ Pivoting ]
calibrate = 1
value = 0
stdev = 1
89
Bijlage C
Gebruik optimalisatiealgoritme
C.1
Principe
Elk door ons in C++ geïmplementeerd genetisch optimalisatiealgoritme, of het nu met één
(Population.hh) of meer (MOEA.hh) objectieffuncties werkt, moet alle methodes van de abstracte ouderklasse Genetics.hh implementeren. Deze klasse levert dus een blauwdruk voor
de genetische algoritmes.
De klasse DataLog.hh houdt voor alle tableaus, per parameterverzameling, per iteratie van het
reconstructiealgoritme, de waarde van een objectieffunctie bij. Voor elke beschouwde objectieffunctie moet een object van de klasse DataLog.hh aangemaakt worden.
We tonen nu de volgorde waarin we de methodes uit Genetics.hh en DataLog.hh moeten
oproepen om één of twee objectieffuncties te optimaliseren binnen de context van een iteratief
reconstructiealgoritme.
C.2
Generieke structuur
DataLog<f l o a t > c o s t 1 _ l o g ( N t o t a l p o p u l a t i o n , N i t e r a t i o n s , Nparameters ) ;
DataLog<f l o a t > c o s t 2 _ l o g ( N t o t a l p o p u l a t i o n , N i t e r a t i o n s , Nparameters ) ;
int tableau = 0 ;
Genetics ∗ c a l i b r a t i o n ;
i f ( c o s t 1 && c o s t 2 ) {
c a l i b r a t i o n = new MOEA( N t o t a l p o p u l a t i o n , N s u r v i v o r s , Nparameters , 2 ) ;
} else {
c a l i b r a t i o n = new P o p u l a t i o n ( N t o t a l p o p u l a t i o n , N s u r v i v o r s , NParameters ) ;
}
f o r ( i n t parameter = 0 ; parameter < Nparameters ; parameter++) {
c a l i b r a t i o n −>S e t C a l i b r a t i o n ( parameter , i n i t _ p a r a m e t e r s [ parameter ] [ 0 ] ) ;
c a l i b r a t i o n −>S e t I n i t i a l P a r a m e t e r ( parameter , i n i t _ p a r a m e t e r s [ parameter ] [ 1 ] ) ;
c a l i b r a t i o n −>S e t S t a n d a r d D e v i a t i o n ( parameter , i n i t _ p a r a m e t e r s [ parameter ] [ 2 ] ) ;
90
BIJLAGE C. GEBRUIK OPTIMALISATIEALGORITME
}
c a l i b r a t i o n −> I n i t i a l i z e ( ) ;
while ( true ) {
c a l i b r a t i o n −>CopyParameters(& c o s t 1 _ l o g ) ;
c a l i b r a t i o n −>CopyParameters(& c o s t 2 _ l o g ) ;
f o r ( i n t i t e r a t i o n = 0 ; i t e r a t i o n < N i t e r a t i o n s ; i t e r a t i o n ++) {
// l o g i c a r e c o n s t r u c t i e a l g o r i t m e
// v a l u e 1 en v a l u e 2 v o o r de r e s p e c t i e v e l i j k e o b j e c t i e f f u n c t i e s b e p a l e n
f o r ( i n t i n d i v i d u a l = 0 ; i n d i v i d u a l < N t o t a l p o p u l a t i o n ; i n d i v i d u a l ++) {
c o s t 1 _ l o g . StoreData ( i n d i v i d u a l , i t e r a t i o n , v a l u e 1 [ i n d i v i d u a l ] ) ;
c o s t 2 _ l o g . StoreData ( i n d i v i d u a l , i t e r a t i o n , v a l u e 2 [ i n d i v i d u a l ] ) ;
}
}
i f ( c o s t 1 && c o s t 2 ) {
c o s t 1 _ l o g . Write ( path , " c o s t 1 " , t a b l e a u + 1 , true ) ;
c o s t 2 _ l o g . Write ( path , " c o s t 2 " , t a b l e a u + 1 , true ) ;
f o r ( i n t i = 0 ; i < N t o t a l p o p u l a t i o n ; i ++) {
c a l i b r a t i o n −>S e t C o s t ( i , c o s t 1 _ l o g . GetLastData ( i ) , 0 ) ;
c a l i b r a t i o n −>S e t C o s t ( i , c o s t 2 _ l o g . GetLastData ( i ) , 1 ) ;
}
} else {
i f ( cost1 ) {
f o r ( i n t i = 0 ; i < N t o t a l p o p u l a t i o n ; i ++) {
c a l i b r a t i o n −>S e t C o s t ( i , c o s t 1 _ l o g . GetLastData ( i ) , 0 ) ;
}
} else {
f o r ( i n t i = 0 ; i < N t o t a l p o p u l a t i o n ; i ++) {
c a l i b r a t i o n −>S e t C o s t ( i , c o s t 2 _ l o g . GetLastData ( i ) , 0 ) ;
}
}
}
i f ( c a l i b r a t i o n −>I s T e r m i n a t e d ( ) ) {
break ;
}
c a l i b r a t i o n −>CreateNewPopulation ( t a b l e a u ) ;
t a b l e a u ++;
}
delete c a l i b r a t i o n ;
91
Bijlage D
Reconstructies
In deze bijlage vergelijken we in meer detail de reconstructie volgens Ω (tabel 10.1) en de reconstructie volgens Ωcal (tabel 10.9).
(a) Ω
(b) Ωcal
(c) Ω
(d) Ωcal
Figuur D.1: Transversale doorsnedes.
92
BIJLAGE D. RECONSTRUCTIES
93
(a) Ω
(b) Ωcal
(c) Ω
(d) Ωcal
Figuur D.2: Coronale doorsnedes.
(a) Ω
(b) Ωcal
(c) Ω
(d) Ωcal
Figuur D.3: Sagittale doorsnedes.
Bibliografie
[1] http://www.cancerresearchuk.org/cancer-help/about-cancer/tests/
ct-scan, “Cancer Research UK: CT Scan”, 14/08/2013, geraadpleegd op 20/05/2014.
[2] http://medical.neusoft.com/en/products/1496/, “NeuViz 16”, 25/04/2014,
geraadpleegd op 05/07/2014.
[3] S. Sawall, M. Knaup, en M. Kachelrieß, “An Adaptive Genetic Algorithm for Misalignment
Estimation (AGAME) in Circular, Sequential and Spiral Cone-Beam Micro-CT”, Proceedings of the 11th International Meeting on Fully Three-Dimensional Image Reconstruction in
Radiology and Nuclear Medicine, 2011, pp. 187-190.
[4] http://deredactie.be/cm/vrtnieuws/binnenland/1.1965141,
“Campagne
Waarschuwt Belgen voor Gevaar van Straling”, 12/05/2014, geraadpleegd op 20/05/2014.
[5] http://deredactie.be/cm/vrtnieuws/binnenland/1.1898501, “Meer MRI- en
PET-Scanners Moeten Blootstelling Bestraling Verminderen”, 03/03/2014, geraadpleegd op
20/05/2014.
[6] D. S. Lalush en M. N. Wernick, “Iterative Image Reconstruction”, Emission Tomography,
Elsevier, 2004, pp. 443-472.
[7] C. Mennessier, R. Clackdoyle, en F. Noo, “Direct Determination of Geometric Alignment
Parameters for Cone-Beam Scanners”, Physics in Medicine and Biology, vol. 54, 2009, pp.
1633-1660.
[8] I. S. Kyprianou, S. Paquerault, B. D. Galas, A. Badano, S. Park, en K. J. Myers, “Framework
for Determination of Geometric Parameters of a Cone Beam CT Scanner for Measuring the
System Response Function and Improved Object Reconstruction”, Proceedings of the IEEE
International Symposium on Biomedical Imaging, 2006, pp. 1248-1251.
[9] S. Brandt, J. Heikkonen, en P. Engelhardt, “Multiphase Method for Automatic Alignment
of Transmission Electron Microscope Images Using Markers”, Journal of Structural Biology,
vol. 133, 2001, pp. 10-22.
94
BIBLIOGRAFIE
95
[10] T. Varslot, A. Kingston, G. Myers, en A. Sheppard, “High-Resolution Helical Cone-Beam
Micro-CT with Theoretically-Exact Reconstruction from Experimental Data”, Medical
Physics, vol. 38, no. 10, 2011, pp. 5459-5476.
[11] D. Panetta, N. Belcari, A. Del Guerra, en S. Moehrs, “An Optimization-Based Method
for Geometrical Calibration in Cone-Beam CT without Dedicated Phantoms”, Physics in
Medicine and Biology, vol. 53, 2008, pp. 3841-3861.
[12] Y. Kyriakou, R. M. Lapp, L. Hillebrand, D. Ertel, en W. A. Kalender, “Simultaneous
Misalignment Correction for Approximate Circular Cone-Beam Computed Tomography”,
Physics in Medicine and Biology, vol. 53, 2008, pp. 6267-6289.
[13] A. Kingston, A. Sakellariou, T. Varslot, G. Myers, en A. Sheppard, “Reliable Automatic
Alignment of Tomographic Projection Data by Passive Auto-Focus”, Medical Physics, vol.
38, no. 9, 2011, pp. 4934-4945.
[14] J. Wicklein, H. Kunze, W. A. Kalender, en Y. Kyriakou, “Image Features for Misalignment
Correction in Medical Flat-Detector CT”, Medical Physics, vol. 39, no. 8, 2012, pp. 49184931.
[15] W. Wein, A. Ladikos, en A. Baumgartner, “Self-Calibration of Geometric and Radiometric Parameters for Cone-Beam Computed Tomography”, Proceedings of the 11th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear
Medicine, 2011, pp. 327-330.
[16] C. Debbeler, N. Maass, M. Elter, F. Dennerlein, en T. M. Buzug, “A New CT Rawdata
Redundancy Measure Applied to Automated Misalignment Correction”, Proceedings of the
12th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology
and Nuclear Medicine, 2013, pp. 264-267.
[17] V. Patel, R. N. Chityala, K. R. Hoffmann, C. N. Ionita, D. R. Bednarek, en S. Rudin,
“Self-Calibration of a Cone-Beam Micro-CT System”, Medical Physics, vol. 36, no. 1, 2009,
pp. 48-58.
[18] A. V. Bronnikov, “A New Algorithm for Geometric Self-Calibration in Cone-Beam CT”,
Proceedings of the 11th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 2011, pp. 175-178.
[19] A. V. Bronnikov, “Reconstruction of Attenuation Map Using Discrete Consistency Conditions”, IEEE Transactions on Medical Imaging, vol. 19, no. 5, 2000, pp. 451-462.
[20] F. Stopp, A. J. Wieckowski, M. Käseberg, S. Engel, F. Fehlhaber, en E. Keeve, “A Geometric Calibration Method for an Open Cone-Beam CT System”, Proceedings of the 12th
International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and
Nuclear Medicine, 2013, pp. 106-109.
BIBLIOGRAFIE
96
[21] Peter Gilbert, “Iterative Methods for the Three-Dimensional Reconstruction of an Object
From Projections”, Journal of Theoretical Biology, 36(1), 1972, pp. 105-117.
[22] B. De Man en S. Basu, “Distance-Driven Projection and Backprojection in Three Dimensions”, Physics in Medicine and Biology, vol. 49, no. 11, 2004, pp. 2463-2475.
[23] P. G. Schulz, “Nelder-Mead Simplex”, Creative Design in Optimization: Metaheuristics
Applied to Multi-Modal Continuous Functions, 2006, pp. 10-14.
[24] J. C. Lagarias, J. A. Reeds, M. H. Wright, en P. E. Wright, “Convergence Properties of the
Nelder-Mead Simplex Method in Low Dimensions”, SIAM Journal on Optimization, vol. 9,
no. 1, 1998.
[25] T. Bäck, “Evolutionary Algorithms in Theory and Practice”, Oxford University Press, 1996.
[26] R. Kumar en Jyotishree, “Blending Roulette Wheel Selection & Rank Selection in Genetic
Algorithms”, International Journal of Machine Learning and Computing, vol. 2, no. 4, 2012,
pp. 365-370.
[27] L. V. Smekal, M. Kachelrieß, E. Stepina, en W. A. Kalender, “Geometric Misalignment
and Calibration in Cone-Beam Tomography”, Medical Physics, vol. 31, no. 12, 2004, pp.
3242-3266.
[28] M.A. Gennert en A.L. Yuille, “Determining the Optimal Weights in Multiple Objective
Function Optimization”, Proceedings of the Second International Conference on Computer
Vision, 1988.
[29] E. Zitzler, M. Laumanns, en S. Bleuler, “A Tutorial on Evolutionary Multiobjective Optimization”, Metaheuristics for Multiobjective Optimisation, vol. 535, Springer, 2004, pp.
3-37.
[30] E. Zitzler en L. Thiele, “Multiobjective Evolutionary Algorithms: A Comparative Case
Study and the Strength Pareto Approach”, IEEE Transactions on Evolutionary Computation, 3(4), 1999, pp. 257-271.
[31] E. Zitzler, M. Laumanns, en L. Thiele, “SPEA2: Improving the Strength Pareto Evolutionary Algorithm”, TIK-Report 103, 2001.
[32] R. Sivaraj, “A Review of Selection Methods in Genetic Algorithm”, International Journal
of Engineering Science and Technology, vol. 3, no. 5, 2011, pp. 3792-3797.
[33] https://developer.nvidia.com/about-cuda,
Cuda”, geraadpleegd op 01/03/2014.
[34] Frank Verhaegen, Maastro Clinic, Nederland.
“NVIDIA CUDA Zone:
About

Manwolf pdf free

Should I Stay or Should I Go? - UC Leuven

null

Kadastraal bericht ob ect

complete workshopoverzicht januari - juni 2015

TYXAL DF - Delta Dore

Bloed moet stromen

structurele vergelijkingsmodellen deel 1

Nieuws - meest recent | SOK | Kortrijk

PLAATSINGSNOTA VAN TE SAMENVATTING

PDF

Ebay-3 2008

ebay

NASIONALE SENIOR SERTIFIKAAT GRAAD 12 GRADE 12

Klik hier voor een PDF overzicht

eBayKopen

Ruimtelijke opgave Atelier 1: draagvlak Atelier 2: voorstellen