Automatisch kalibreren van CT-scanners tijdens iteratieve beeldreconstructie

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


. . . . . . . . . . . . . . . . . . . . . . .

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


. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

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 beelden 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


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.


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


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


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.


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 beelden 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


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


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 beelden 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 beelden/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.


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


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


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


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


(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.


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

◦


-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


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-


(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.


(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.


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.


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


(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.


(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)


(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,


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


(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.


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)


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.


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-


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


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


(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


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


(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.


(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


(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.


(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


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.


(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.


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


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


(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


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.


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


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-


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)


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.


(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


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-


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.


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.


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.


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.


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 e s t i m a t e ; P r o j e c t i o n s S e t measured ; ProjectionsSet_GPU 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 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.


Symbool


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

◦


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



0,0 0,0 0,0 -30,0 -162,0 0,0 293,0 140,0


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.


(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.


(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


(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.


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 beelden 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.


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


(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.


(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.


(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


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



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.


(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.


78

Symbool

Beschrijving

Waarde

Eenheid

Wijzigen



0,0 0,0 0,0 -37,8105 -162 -0,2621 293,0 140,0


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.


79

Symbool

Beschrijving

Waarde

Eenheid

Wijzigen



0,0 0,0 -4,22209 -37,8105 -151,993 -0,2621 293,0 140,0


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.


80

Symbool

Beschrijving

Waarde

Eenheid



0,0 0,0 -4,22209 -37,8105 -151,993 -0,2621 305,6862 127,6113


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.


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 beelden 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.

Bĳlage 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

BĲLAGE 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

BĲLAGE 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

Bĳlage 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

BĲLAGE 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

BĲLAGE 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

Bĳlage 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 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 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

BĲLAGE 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

Bĳlage 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

BĲLAGE 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

Automatisch kalibreren van CT-scanners tijdens iteratieve beeldreconstructie

Recommend Documents