Segmentation and Alignment of 3-D Transaxial
Myocardial Perfusion Images and
Automatic Dopamine Transporter Quantification
Reg nr: LiTH-IMT/MI30-A-EX--08/465--SE
Leo Bergnéhr
Segmentation and Alignment of 3-D Transaxial
Myocardial Perfusion Images and
Automatic Dopamine Transporter Quantification
Master’s Thesis in Biomedical Engineering
at Linköping University
by
Reg nr: LiTH-IMT/MI30-A-EX--08/465--SE
Leo Bergnéhr
April, 2008Supervisor:
Abstract
Nuclear medical imaging such as SPECT (Single Photon Emission Tomography) is an imaging modality which is readily used in many applications for measuring physiologi-cal properties of the human body. One very common type of examination using SPECT is when measuring myocardial perfusion (blood flow in the heart tissue), which is often used to examine e.g. a possible myocardial infarction (heart attack). In order for doctors to give a qualitative diagnose based on these images, the images must first be segmented and rotated by a medical technologist. This is performed due to the fact that the heart of different patients, or for patients at different times of examination, is not situated and rotated equally, which is an essential assumption for the doctor when examining the images. Consequently, as different technologists with different amount of experience and expertise will rotate images differently, variability between operators arises and can often become a problem in the process of diagnosing.
Another type of nuclear medical examination is when quantifying dopamine trans-porters in the basal ganglia in the brain. This is commonly done for patients showing symptoms of Parkinson’s disease or similar diseases. In order to specify the severity of the disease, a scheme for calculating different fractions between parts of the dopamine transporter area is often used. This is tedious work for the person performing the quan-tification, and despite the acquired three dimensional images, quantification is too often performed on one or more slices of the image volume. In resemblance with myocardial perfusion examinations, variability between different operators can also here present a possible source of errors.
In this thesis, a novel method for automatically segmenting the left ventricle of the heart in SPECT-images is presented. The segmentation is based on an intensity-invariant local-phase based approach, thus removing the difficulty of the commonly varying in-tensity in myocardial perfusion images. Additionally, the method is used to estimate the angle of the left ventricle of the heart. Furthermore, the method is slightly adjusted, and a new approach on automatically quantifying dopamine transporters in the basal ganglia using the DaTSCANTMradiotracer is proposed.
Sammanfattning
Nukleärmedicinska bilder som exempelvis SPECT (Single Photon Emission Tomogra-phy) är en bildgenererande teknik som ofta används i många applikationer vid mätning av fysiologiska egenskaper i den mänskliga kroppen. En vanlig sorts undersökning som använder sig av SPECT är myokardiell perfusion (blodflöde i hjärtvävnaden), som ofta används för att undersöka t.ex. en möjlig hjärtinfarkt. För att göra det möjligt för läkare att ställa en kvalitativ diagnos baserad på dessa bilder, måste bilderna först seg-menteras och roteras av en biomedicinsk analytiker. Detta utförs på grund av att hjärtat hos olika patienter, eller hos patienter vid olika examinationstillfällen, inte är lokaliserat och roterat på samma sätt, vilket är ett väsentligt antagande av läkaren vid granskning av bilderna. Eftersom olika biomedicinska analytiker med olika mängd erfarenhet och expertis roterar bilderna olika uppkommer variation av de slutgiltiga bilder, vilket ofta kan vara ett problem vid diagnostisering.
En annan sorts nukleärmedicinsk undersökning är vid kvantifiering av dopaminrecep-torer i de basala ganglierna i hjärnan. Detta utförs ofta på patienter som visar symptom av Parkinsons sjukdom, eller liknande sjukdomar. För att kunna bestämma graden av sjukdomen används ofta ett utförande för att räkna ut olika kvoter mellan områden runt dopaminreceptorerna. Detta är ett tröttsamt arbete för personen som utför kvantifierin-gen och trots att de insamlade bilderna är tredimensionella, utförs kvantifierinkvantifierin-gen allt för ofta endast på en eller flera skivor av bildvolymen. I likhet med myokardiell perfu-sionsundersökningar är variation mellan kvantifiering utförd av olika personer en möjlig felkälla.
I den här rapporten presenteras en ny metod för att automatiskt segmentera hjärtats vänstra kammare i SPECT-bilder. Segmenteringen är baserad på en intensitetsinvariant lokal-fasbaserad lösning, vilket eliminerar svårigheterna med den i myokardiella perfu-sionsbilder ofta varierande intensiteten. Dessutom används metoden för att uppskatta vinkeln hos hjärtats vänstra kammare. Efter att metoden sedan smått justerats används den som ett förslag på ett nytt sätt att automatiskt kvantifiera dopaminreceptorer i de basala ganglierna, vid användning av den radioaktiva lösningen DaTSCANTM.
Acknowledgements
This project has been carried out at Exini Diagnostics AB and I would like to thank everyone there for their support and friendliness. Special gratitude goes to my company supervisor, professor Lars Edenbrandt, who initiated the project and invited me to work on it. Thank you; bringer of good confidence!
Additionally, I would like to thank my university supervisor, Joakim Rydell, at the De-partment of Biomedical Engineering. In the beginning of this project I had a lot of questions, to which you always had enlightening answers, despite completing a PhD at the same time.
Finally, I give a special acknowledgement to my family and near friends; not for helping me while formulating mathematical expressions, but for just being there when I was not. Thank you!
Table of Contents
1 Introduction 1 1.1 Background . . . 1 1.2 Objectives . . . 2 1.3 Method . . . 2 1.4 Outline . . . 3 1.5 Notation . . . 4 2 Scintigraphy 5 2.1 Nuclear medicine physics . . . 52.1.1 Decay of radioactive isotopes . . . 5
2.1.2 Detection of radiation . . . 6
2.2 Single Photon Emission Computed Tomography . . . 8
2.2.1 Image quality . . . 8
2.3 Scintigraphic examinations . . . 9
2.3.1 Myocardial perfusion SPECT . . . 9
2.3.2 DaTSCAN Quantification . . . 10 2.4 Summary . . . 11 3 Image Segmentation 13 3.1 Segmentation methods . . . 13 3.1.1 Thresholding . . . 13 3.1.2 Region growing . . . 14 3.1.3 Active shapes . . . 15 3.2 Image registration . . . 15 3.3 Summary . . . 17
4 The Morphon Method 19 4.1 Displacement field estimation . . . 20
4.1.1 Quadrature phase . . . 20
4.1.2 Least square solution . . . 24
4.2 Displacement field accumulation . . . 24
4.3 Displacement field regularisation . . . 25
4.4 The model . . . 26
4.5 Summary . . . 27
5 Rotation of Myocardial Perfusion Images 29 5.1 Pre-processing of images . . . 29
5.2 2-D Morphon segmentation . . . 30
5.2.1 Locating the heart . . . 31
5.2.2 Segmentation . . . 31
5.3 3-D Morphon segmentation . . . 32
5.4 Angling of the heart . . . 33
5.4.1 Angle estimation . . . 34
5.4.2 Rotation . . . 35
5.5 Summary . . . 35
6 Dopamine Transporter Quantification 37 6.1 3-D Morphon segmentation . . . 37
6.1.1 Locating the basal ganglia . . . 38
6.1.2 Segmentation . . . 39
6.2 Quantification . . . 40
6.2.1 Weighted principal component analysis . . . 41
6.3 Summary . . . 43
7 Results 45 7.1 Segmentation and alignment of 3-D transaxial myocardial perfusion im-ages . . . 45
7.1.1 Segmentation . . . 45
7.1.2 Alignment . . . 46
7.2 Automatic dopamine transporter quantification . . . 46
7.2.1 Segmentation . . . 52
7.2.2 Quantification . . . 52
8 Discussion 59 8.1 Segmentation and alignment of 3-D transaxial myocardial perfusion im-ages . . . 59
8.2 Automatic dopamine transporter quantification . . . 61
1
Introduction
1.1 Background
Scintigraphy is a widely used medical imaging modality for observing physiological properties in the human body. Two examples that are considered in this thesis are scinti-graphic heart images and images of the basal ganglia in the brain.
At the time of writing, a scintigraphic examination of the heart is at hospitals carried out by a technologist. The main tasks for the technologist is to prepare the patient for the diagnostic procedure, and later on to perform a sequence of manual operations on the acquired images. These operations includes positioning the location of the heart to a predefined area in the images, determining the orientation of the heart and concluding whereas the image quality is sufficient for the doctor to base a diagnosis on. The reori-entation of the heart is necessary due to the simple fact that different individuals do not look the same; not on the outside, nor on the inside of the body. In heart scintigraphy, this presents a problem.
For a doctor to give a qualitative diagnosis, the result of an examination is compared to other images. This may occur directly, by looking at different images from other patients, or indirectly, by taking into account the vast amount of images in the past that
is the experience of the doctor. However, both these methods depend on the presumption
that the images have been observed from the same point of view, which would not be the case if the technologist did not manually compensate for angular and rotational differences of the orientation of the heart, between different scintigraphic examinations. However, as different technologists have different experience and slightly different ways of making these adjustments, the result will not be consistent for all examinations. An automatic process for doing this would not only hopefully make the result more con-sistent, but also save time that the technologist has to spend on performing these tasks manually. One method for doing this has been proposed by Germano et al. [1], hoping
to move more and more into a fully automated process [2]. However, as this method is based on intensity of the images, a more insensitive method may be desired.
Moreover, scintigraphic images of the basal ganglia in the brain are beginning to be more and more popular to diagnose Parkinson’s disease. At the time of writing, quan-tification is used to measure a number of properties from the images of the basal gan-glia, thus being able to determine the state of the patient. This process is tedious and can present difficulties if the state of the basal ganglia is poor. An automatic method for doing this quantification would benefit in much the same ways as the case with scinti-graphic heart examinations. The process would be much faster compared to manual quantification, and the variability between medical staff would be eliminated. Addi-tionally, manual quantification is performed in two dimensions and hence it could be beneficial to quantify in 3-D. Morton et al. [3] has evaluated two computer driven meth-ods in the same field and found that a 3-D based method gives better concordance with visual assessments.
1.2 Objectives
The primary aim of this project is to investigate if an automatic reorientation of three-dimensional, transaxial myocardial perfusion SPECT images is possible, and if so, to design and implement an automatic segmentation method that can distinguish, and com-pensate for, differences in anatomical structure of the heart; particularly the angle and position of the left ventricle. The method should be able to segment and reorient the left ventricle from transaxial images into short-axis (oblique) images.
Secondarily, the same method should, after little modification, be applicable to scinti-graphic images of the basal ganglia. These should be segmented, and a variety of useful parameters should be extracted. These parameters should correspond to parameters that are extracted in manual quantification.
1.3 Method
The work in this project has been divided into two different parts. First, the segmen-tation algorithm was implemented, and secondly the algorithm was applied on the two problems at hand, which first required some pre-processing in each case. Computa-tions and simulaComputa-tions have been carried out with the help of MATLAB!R, version 6 and
7, from The MathWorks, Inc. The main part of the functions and methods used has been implemented with the help of MATLAB-included primitive functions. However,
to decrease computation time, a C-compiled version of convolution has been used, with courtesy of Gunnar Farnebäck at the Department of Biomedical Engineering, Linköping
1.4 Outline 3 university. To decrease the time of implementation, some functions from the MAT -LAB Image Processing Toolbox have also been used. These are the bwlabel() and
regionprops() for some simple morphological operations, and tformarray() for rotating 3-D volumetric objects.
All the images used for development of the algorithm have been provided by Swedish hospitals. The scintigraphic heart images originate from Sahlgrenska university hospital in Gothenburg, Sweden, while the DaTSCAN-images of the basal ganglia have been provided by The university hospital of Norrland in Umeå, Sweden.
1.4 Outline
Here follow a short summary of the contents of each chapter of this thesis.
Chapter 2 The basics of nuclear physics is presented along with the essential parts of a nu-clear medical examination. In this chapter, relevant facts about myocardial perfu-sion examinations and DaTSCAN examinations are also introduced.
Chapter 3 In this chapter, image segmentation and registration is introduced. Some com-monly used techniques are presented along with their good and poor properties. Chapter 4 The main method in this thesis, the Morphon method, is gone through in some
detail in this chapter.
Chapter 5 Here the problem with the estimation of angles of the heart in myocardial perfu-sion images is addressed with the help of the earlier presented methods and some complementary techniques.
Chapter 6 This chapter involves applying the Morphon method on the problem of automatic quantification of dopamine transporters and introduces some additional tools for solving this.
Chapter 7 The results of both parts of the thesis are presented in this chapter.
Chapter 8 The results of the previous chapter are discussed and some future development of the methods is proposed.
Chapter 9 This chapter concludes the whole thesis, its results and possible effects on the problems at hand.
1.5 Notation
This thesis inevitably involves some mathematics. The notations used are presented below. s, S Scalar v Vector ˆ v Unit vector |v| Norm of vector M Matrix.
vT, MT Transpose of vector and matrix
M−1 Inverse of matrix
s∗ Complex conjugate of the complex signal s
(s1∗ s2)(x) Convolution of s1 and s2
2
Scintigraphy
To be able to understand the process of generating scintigraphic images using nuclear medical methods such as SPECT (Single Photon Emission Computed Tomography), the basic physics behind it is laid forward in this chapter. Details on the physical aspects and more information regarding different isotopes can be found in [4] and [5, 6] respectively.
2.1 Nuclear medicine physics
Nuclear physical applications in medicine have been in clinical use for many decades. The first experiments started somewhere around 1935, when Chiewitz and de Heversy [7] tested the technique on rats. At first, nuclear medicine was mainly used in thera-peutical treatments, such as that of cancer, introduced by Lindgren [8]. However, as acquisition instruments and radiotracers improved, it also became a very popular way of examining physiological and anatomical properties of many parts of the human body, for example the heart [9] and the brain [10].
2.1.1 Decay of radioactive isotopes
In all imaging modalities, some kind of conceptual probe is used to gather information about the observed object at hand. This probe can be very different in different situa-tions. For example, at an ultra sound examination, mechanical waves are used to gather information about the object, by propagation and reflection between different acous-tic impedances, and in magneacous-tic resonance imaging, radio frequency fields, interacting with fundamental properties of matter, are the probes at hand. In nuclear imaging,
radi-ation is used to acquire informradi-ation about an object. Radiradi-ation, which originates from
radioactive isotopes1.
1An isotope of an element is a specific variation of that element, where the number of neutrons in the
There are a variety of radioactive isotopes used in nuclear medicine. Depending on the type of examination and the region of interest in the body, the energy of the decaying particles from the isotopes has to be taken into consideration. For example, in SPECT, which is the modality of interest in this thesis, the most common isotope used is99mTc
(technetium), which is suitable due to, among other things, its half-life and the energy of its emitted photons [11]; that is about the same as for conventional X-ray γ-energies2.
A radioactive isotope (radioisotope) is, by definition, unstable. That means that the nucleus of that element can, and will, spontaneously decay with a certain probability, given a period of time. The faster the nucleus decays, the faster will it reach its half-time value; that is when 50 percent of the amount of isotope has vanished in favour of decay radiation of some sort. There is a tremendous amount of different isotopes, both occur-ring naturally on Earth and synthetically produced with the help of nuclear reactions. There are also many different types of decaying processes for different types of nucleus. These different types of decay have been categorised after the type of particle emitted in the process. The most known particle emission, or radiation, is γ-radiation, where the emitted particles are photons, or light. Other common types of radiation are α-, β−- and
β+-radiation [4]. β- and γ-radiation is of much interest from a nuclear medical point of
view since they are the key particles in SPECT, while α-radiation is attenuated quickly in most media and is more harmful to human tissue, and therefore is used sparsely in medical diagnostics [5]. β−/+-particles are more commonly known as electrons (e−)
and positrons3(e+).
2.1.2 Detection of radiation
As was stated earlier, γ- and β-radiation is preferable in nuclear medical diagnostics, however the methods for collecting images with the two methods differ. In SPECT imaging, isotopes like 99mTc, mentioned in the first paragraph of 2.1.1, are used for
different purposes. While some have been stated, the one more important for the basic principles of this imaging modality is that these kind of isotopes decay with a single emission of a γ-particle.
Since the nucleus decays spontaneously, it will have an equal probability of emitting a photon in any direction [4], which is illustrated in Figure 2.1. The emitted photon will have a total energy of about 140 keV, which is suitable for many photon detectors at the time of writing [6]. To detect the emitted γ-rays, different types of detectors can be used. However, particularly in nuclear medicine scintigraphy, a common method is to use a so called scintillator for this purpose. A scintillator is basically a type of
2The Greek letter γ is in physics often used as a symbol for photon or light.
3The positron is the antiparticle of the electron, which means it has exactly the same properties, except
2.1 Nuclear medicine physics 7
!
Figure 2.1: The one-step process when a radioactive isotope decays with γ-radiation.
crystal4with specific characteristics. These crystals have in common that they undergo
a process called fluorescence when interacting with γ-rays of an energy matching the characteristics of the crystal. Fluorescence is a process where an element will absorb a γ-ray and immediately emit light within the visible range. This is due to the photo electric effect, where electrons in an atom can be pushed to a higher energy state and then, when falling back to the original state, send out a flash of light [4]. This light can later be amplified, for example with the help of a photomultiplier, and then detected. Guerra [6] points out that while the scintillator is more sensitive and therefore preferable in many applications, other detectors are used, such as gas-filled chambers and solid state semiconductors.
!
1
2
e+ e+e-Figure 2.2: The two-step process when a radioactive isotope decays with a positron (β-decay), which shortly thereafter reacts with an electron, creating γ-radiation.
While single photon emitters (SPE) such as99mTc are used to detect one γ-ray at a time,
4Common crystals that are used in nuclear medicine are: Cerium doped lutetium orthosilicate (LSO),
there is another method used in nuclear medicine that utilise positron emitters (PE), such as Oxygen-15. In accordance to the name, positron emitters decay by emitting positrons; that is with β+-decay (see figure 2.2) [5]. However, as an antiparticle to the
electron, the positron will readily react with electrons it encounters, annihilating both particles and sending out two photons near exactly in the opposite direction to each other. In vivo, the positron will react almost instantly with an electron, producing two photons that will radiate from the same source [4]. Consequently, even though another type of nuclear decay is the source of radiation, photons will, in the end, be detected by both type of radioisotopes. In nuclear medicine, the imaging modality using PE is referred to as PET (Positron Emission Tomography). PET-images are very similar to SPECT-images, but they will not be studied any further in this thesis due to the lack of PET-image material.
2.2 Single Photon Emission Computed Tomography
The imaging modality that is the source of all images in this thesis is SPECT. In SPECT, photons from single photon emitters are detected in a manner as described in 2.1.2. SPECT is a widely used modality and has been used for some time. One of the great ad-vantages of SPECT compared to projectional modalities like X-ray imaging is its ability to acquire images in 3-D, thereof the word tomography in the abbreviation. Larsson [12] presented one of the first working tomographic systems and the image quality has been much improved ever since.
2.2.1 Image quality
What should be mentioned about the images generated in a SPECT examination is that they are low in resolution and signal to noise ratio compared to images from other imaging modalities, such as CT (computed tomography) and MRI (magnetic resonance imaging). A typical SPECT-image has a resolution of 643 voxels5, while the other
modalities can have resolutions of 5123 voxels and above [5]. The low resolution in
SPECT-images originates from the fact that the photon count, i.e. the number of photons that are detected, is quite low and that SPECT registers physiological processes in the body, not anatomical structures as CT and MRI are mainly used for. Further more, even though the photon count could probably be increased with a more radiating radioisotope, it would consequently increase the dose of radiation applied to the patient.
2.3 Scintigraphic examinations 9
2.3 Scintigraphic examinations
There are many steps involved in a clinical scintigraphic examination. In this section, the process of examining myocardial perfusion as well as dopamine activity in the basal ganglia, and some problems that arrises with manual diagnoses in such examinations, will be presented.
In most nuclear medical applications the radioisotope that generates the image signal is tagged to a compound of some sort, called a pharmaceutical agent, which is injected into the patient. Together, the isotope and the agent are named a radiotracer, which is used to attach itself to the organ that is to be investigated, and to make the radioisotope less hostile to the body. For example, as exemplified by Wilson et al. [13], in myocardial perfusion SPECT examinations, tert-butyl is commonly used, and for the brain images studied in this thesis, the popular DaTSCANTM radiotracer from GE Healthcare [14]
has been used.
2.3.1 Myocardial perfusion SPECT
Myocardial perfusion SPECT examinations are widely used as a routine examination for patients with chest pain. Some time after the technologist has injected the radiotracer into the patient, he/she will be placed in a camera station, often referred to as γ-camera, with a detector suitable for SPECT γ-energies [6]. After the camera has finished with the acquisition the image is reconstructed using some reconstruction technique, e.g. as proposed by Perkins [5]; filtered back-projection or iterative reconstruction. This will
Figure 2.3: A 3-D volume image displayed as iso-surfaces (left) and two slices (center and right), with the left ventricle of the heart marked with an arrow.
yield the images, which are tomographic 3-D volume data, i.e. each voxel in space has three space-coordinates (x, y, z) as well as an intensity level. In practice, the 3-D volume consists of a number of slices. This makes it convenient to display the volume in the form of a stack of 2-D images, which is often used clinically. Strauss and Miller
[15] describes the procedure in more detail. An example of a 3-D volume acquired using SPECT can be found in Figure 2.3.
The next task for the technologist is to rotate the acquired images of the heart into predefined angles. This is done by choosing two slices from the volume, which the technologist finds good for estimating the angle of the heart [15]. An example is given in Figure 2.4.
Figure 2.4: The manual rotation of the heart performed by the technologist. The heart is here illustrated before (first row) and after (second row) rotation, for slices through each axis respectively.
It is obvious that the centre line of the left ventricle, which is commonly seen as the line that should align with the axes after rotation [15], cannot be exactly defined for such an asymmetrical object as the left ventricle. Even more so when the images acquired are of poor quality, or if the uptake of radiotracer in the left ventricle is poor, due to, e.g., myocardial infarction. This poses a problem for myocardial perfusion SPECT images, and is discussed in this thesis.
2.3.2 DaTSCAN Quantification
DaTSCAN quantification is another SPECT-application that is very commonly used to diagnose Parkinson’s disease. If a patient suffers from Perkinson’s disease, the basal ganglia are often degenerated, i.e. the number of dopamine receptors are reduced. More-over, the degeneration is generally first observed from the same place of the ganglia. However, the quantification methods that are used introduces great inter- and intraop-erator variability, which is discussed by Tatsch et al. [16]. Often, a two-dimensional quantification method with predefined or adjustable regions of interest (ROI) is used,
2.4 Summary 11 which requires operator interaction throughout the whole quantification process. Visual interpretation of the images are also common, which introduces even more sources of error, e.g. the use of different colour maps; noted in [16] as:
Use of non-continuous color tables may overestimate findings due to abrupt color changes.
Tatsch et al. [16] As with myocardial perfusion SPECT images (section 2.3.1), the variability is great between different images. An example of how the basal ganglia are manually quantified is given in Figure 2.5, where a typical ROI has been drawn around each of the basal ganglia, as well as for the background signal of the image. The quantification of the
Figure 2.5: Illustrating a typical ROI of the basal ganglia and background area (4) of a manually quantified DaTSCAN-image.
images is commonly based on measuring the signal inside the ROI (or different parts of it) of the basal ganglia, and comparing it to some kind of background signal. The background signal has traditionally, with 2-D quantification, been chosen as the visual cortex (see nr. 4 in Figure 2.5), but there is no absolute standard [16]. The quantification procedure can take a significant amount of time for an inexperienced operator and there is little standardisation in the way that the different ROI:s are selected and adjusted.
2.4 Summary
After this chapter, the reader should be sufficiently informed on how SPECT-imaging works, both theoretically, physically and on a more practical level, which has been ex-emplified with the two applications that are involved in this thesis; myocardial perfusion
SPECT and DaTSCAN quantification. Problems that both applications have in common regarding inter- and intraoperator variability as well as dimensional reduction approxi-mations from 3-D to 2-D, have also been introduced.
3
Image Segmentation
Image segmentation is in principle a method for bringing out certain, more interesting, parts of a signal, and has been used in medical applications for many years. The signal may for example consist of a two dimensional image, or as in the case of this thesis, a three dimensional volume. Both cases are common situations in medical informatics. In this chapter, a brief overview of common image segmentation techniques will be cov-ered, and since the principles are the same for signals of different dimension, examples will mainly concern two dimensional images, since they are easier to visualise. Image segmentation in more detail can be found in [17].
3.1 Segmentation methods
There are many different image segmentation methods. From the very simple and in-tuitive, to techniques based on advanced mathematical concepts. An important issue to consider is therefore what type of segmentation method to apply, given a certain situa-tion. This often depends on the image modality used to acquire the image, as well as the complexity of the part of the image that is to be segmented.
3.1.1 Thresholding
In a very simple and clear image, simple and intuitive segmentation methods can often be sufficient. In such images, with a good signal to noise ratio, so called thresholding is often a good choice. Segmentation by thresholding works simply by keeping the part (pixels/voxels) of the image that lies in a certain intensity interval. This is demonstrated in Figure 3.1 b), where the pixels above a certain intensity level in the original image has been set to 1, and the rest of the image to 0, and can also be described mathematically
with
Isegmented =
!
1, a < I < b
0, o.w. , (3.1)
where I is the original image and a and b are the thresholds. However, if the image signal is bad and mixed with noise, thresholding will often yield a bad result. This is shown in Figure 3.1 c, where noise has been added to the original image.
Figure 3.1: a) (left): An X-ray image of the hand of Wilhelm Röntgen’s wife. b) (centre): The image after thresholding. c) (right): The image after thresholding where noise has been added to the original image.
Furthermore, another problem with thresholding is choosing a good threshold level that is not specific to one image. This can be done by for example examining the histogram of the image, or by other more or less sophisticated guesses. At last, segmentation by thresholding is often not a good choice where the interesting parts of an image has the same intensity levels as the rest of the image, and is only separated by structure and shape. In such situations, some kind of image registration is often required for good results.
3.1.2 Region growing
A segmentation method that resembles thresholding is region growing [17]. However, in contrast to thresholding, which operates on each point in an image independently of neighbouring points, region growing takes into account the surrounding connected pixels. Region growing requires a starting point in the image. This point should be a known point within the area of interest, and consequently the method requires manual
3.2 Image registration 15 input or some kind of pre-processing. When the point is chosen, the method works iteratively. A neighbouring point is chosen. If the intensity of that point is within a certain intensity level from the starting point, the point is considered part of the final segmentation, and the process starts over, now starting at this point. There are several more or less advanced techniques for implementing region growing, however, it often gets very slow with higher image dimensions. More over, region growing can in its basic form handle images with much noise very poorly. There are ways of reducing noise, but nevertheless, region growing will not take into account the structure of the interesting parts of an image and is highly dependent on image intensity.
3.1.3 Active shapes
The thought behind active shapes takes region growing a step further by introducing prior knowledge of the the object that is to be segmented. The goal in active shapes is to find a line in 2-D or surface in 3-D, that encloses the object of interest after segmen-tation. Active shapes works by first defining a model. The model is built from a data set of other images of the same kind of object that is going to be segmented. In each one of the images in the data set, the same landmarks (often anatomical) have been marked. By observing how the position of the landmarks differ throughout the data set, the vari-ance of each landmark can be calculated [17]. In turn, this is used to define how the model is allowed to move while it is segmented. The actual deformation of the model can be generated in various ways. One common way is to use the energy in the image to calculate edges and borders, to where the model should deform. Compared to region growing techniques, active shapes are much more robust and insensitive to noise. This is due to the prior knowledge introduced about the object, which prevents the model from deforming in an unnatural way. The constraint of active shapes is off course the data set. If the data set is small, the model will not be able to deform sufficiently when introduced to cases that are not represented in the data set. However, with a large data set, active shapes is often a good method in more complicated segmentation situations.
3.2 Image registration
In many situations, particularly in medical imaging, simple segmentation methods such as thresholding is not sufficient to yield good results. The segmentation process can then often be helped by adding image registration to the segmentation process. Image registration is in principle a method for deforming one image to align to another image. In medical informatics this is often used to map one image, with one type of informa-tion, to another image, containing other types of information. One example is to map an image of a brain acquired with computed tomography (CT) or magnetic resonance
imaging (MRI) to another image of the same brain, acquired with SPECT. Since the brain may not be in the same position at both times of examination, image registration is used to make the same areas of the two images correspond to the same areas of the brain. In this example, information about the anatomy of the brain, in the CT or MRI images, will then be complemented with physiological information from the SPECT image.
Mathematically, image registration can be described by finding a deformation field, or a displacement field, v(x), that after acting on an image, I2(x)will deform the image in
such a way that it will map to the target image1, I
1(x). The mapping in n dimensions
can be described by
v(x) : Rn"→ Rn. (3.2)
Moreover, to be able to describe the quality of the registration process, the error
$2 =$I
2(x + v(x))− I1(x)$2 (3.3)
is introduced. If this error is zero, the registration process has been completely success-ful. This is not a realistic goal in practice, since the two images do not contain exactly the same information, in which case image registration would not be necessary. Conse-quently, v(x) describes how each point in I2 should be moved to minimise $. There are
many attempts of solving this problem, and many involve minimising some variation of a least square problem. For example, if the following assumptions are made:
1. v(x) can locally be described solely by a spatial displacement.
I(x, t) = I(x + ∆x, t + 1), (3.4)
where t denotes the time in the discrete registration process.
2. The image can in each point be approximated with a first order Taylor expansion.
I(x + v(x), t + 1) = I(x, t + 1) +∇ITv, ∇I = [∂I/∂x, ∂I/∂y]T. (3.5)
If equations (3.4) and (3.5) are combined, the equation for the so called optical flow, denoted by v, can be found.
∇ITv
− ∆t= 0, ∆t = I(x, t)− I(x, t + 1). (3.6)
Optical flow is one of the oldest concepts in image registration and computer vision and a lot of research has been done around it. Consequently, there are many methods for
3.3 Summary 17 solving the differential equation (3.6). As a matter of fact, there are many solutions for
vin each point in the image; solutions that are not practically usable, since image points
can then have completely independent solutions of each other. This is described in more detail in [17], chapter 6:6. However, one robust solution that consider the optical flow as a continuous function is a least square problem formulation:
$2 ="
i
([∇I(xi)]Tv(xi)− ∆t(xi)) (3.7)
By assuming different properties for v in equation (3.7), different types of solutions for the displacement field can be found, with different amount of rigidity [18]. For example, in some situations a mere translation or rotation of v is a good enough reg-istration model, while in other more complex applications, a non-rigid model must be used, which allows for any kind of deformation of I2. The parameterized solutions for
different levels of rigidity can be found in [17].
In this section, the optical flow was found through differences in the intensity of im-ages. While straightforward, this approach is obviously dependent on the assumption in (3.4). There are, however, other more robust procedures for estimating the optical flow in a registration process between two images. Differences in local phase is one such approach, and will be described in more detail in chapter 4.
3.3 Summary
In this chapter, an overview of more and less advanced segmentation methods have been presented. From thresholding, which is easy to understand and implement, to registration based methods; derived within a rich mathematical framework and still issue to much research and algorithm optimisation.
4
The Morphon Method
The Morphon method is a newly developed method for image segmentation and registra-tion, first introduced by Knutsson and Andersson [19]. Since then, the method has been under refinement and have been applied in several medical applications [20, 21, 22, 23]. This chapter is mainly an overview of the Morphon method, based on [19] and [20]. The method can be brought to a much more advanced level, nevertheless, this simplified (and hopefully faster) version is sufficient for the subject of this thesis.
The Morphon method is essentially a registration driven segmentation technique. It works by first estimating a displacement field through differences in local phase be-tween the segmentation target image and a model image, which resembles the target image in e.g. intensity variations. The displacement field is further processed to ensure robustness of the method. The model of the object that is to be segmented is deformed according to the displacement field, and the whole process is iterated until the process is satisfactory and the model has mapped on to the target image.
In the beginning of the work of this thesis, a method for solving the stated problem in a satisfactory way had not yet been proposed. After investigating several image segmen-tation and registration methods, the conclusion was finally that a phase based approach such as the Morphon method had several advantages in comparison, given the prob-lematics of this project. The primary reasons for choosing the Morphon method can be summarised with the following arguments.
1. Phase: As mentioned in section 3.2, a phase based approach on image registration can be advantageous compared to classical techniques with intensity based optical flow estimation. Even more so in the images that is considered in this thesis. Section 2.3.1 shows that SPECT images of the heart can vary considerably in intensity. Not only among different images, but also between the heart and other organs. In such situations, intensity based techniques are inferior to a method such as the Morphon, which, as presented in 4.1.1, is insensitive to intensity variations.
2. Noise: If the images are of low quality1, which is not an unusual case with SPECT
images (see section 2.2.1), assumption 2 (equation (3.5) might very well not hold. However, the Morphon method uses robust methods such as displacement field accumulation and regularisation to prevent the displacement field from diverging. 3. Model: The Morphon framework includes an easy and dynamic way for includ-ing prior knowledge about the object that is to be segmented. This makes it easy to develop further for other applications, as will be presented in chapter 7,
Seg-mentation and quantification of DaTSCAN images.
4.1 Displacement field estimation
A displacement field is essentially a vector field. In this approach, the field is dense, which means that in each point in space, there exists a vector with a certain direction and magnitude. Since the image space is discretised with pixels or voxels, the vector field will have the same size as the image that is to be deformed by it. There are, however, other applications, where the deformation field is not dense, and a vector describes the movement of a block of discrete points. To estimate a displacement field in the Morphon method, analysis of the local phase of the images is made.
4.1.1 Quadrature phase
Phase can principally be described as the structure of the image. Mathematically, the phase, as it is used in this thesis, is the response from a certain set of filters that respond to edges and borders in the image. These filters are known as quadrature filters and are defined as follows.
Quadrature filter 4.1.1. A filter, f(x), is a quadrature filter if its Fourier transform,
F (ω), is zero on one side of a hyper plane through the origin in the Fourier domain. That is, there exists a vector, ˆn, such that
F (ω) =
!
0, nˆTω≤ 0
F (ω), o.w. . (4.1)
In this thesis, one special type of quadrature filters have been used, namely lognormal2
filters, which are described, for example in [24], in the Fourier domain as
F (ω) =e−B24ln2ln2(
ω
ω c) (4.2)
1Quality can be hard to define. Here, quality mainly refers to signal-to-noise ratio and the amount of
artifacts.
4.1 Displacement field estimation 21 for positive values of ω. Here, B is the relative bandwidth of the filter in octaves and ωc
is the centre frequency of F (ω). In Figure 4.1, a lognorm filter is shown in the spatial and frequency domain. For displaying reasons it is shown for only 1 dimension, but the principle is the same for arbitrary dimensions. Due to the quadrature filter definition in the frequency domain, quadrature filters are always complex in the spatial domain. This can be understood since
F (ω)'= F (−ω), (4.3)
while a necessary requirement for f(x) to be real is that
F (ω) = F (−ω). (4.4) ! " #! #" $! $" %! !!&# !!&!" ! !&!" !&# !&#" '()* +,)-./)01 !! !" !# $ # " ! $ $%# $%" $%! $%& $%' $%( $%) $%* $%+ #
Figure 4.1: Lognormal filter with real and imaginary parts in the spatial domain (left) and the Fourier domain (right)
In computational practice, filters are discrete, and it is often desirable to adjust the fi-nite length of the quadrature filter. A smaller filter will imply less computations while making calculations with the filter. Moreover, the size of the filter will have impact on the size of the items that the filter is able to detect. The size of the filter can, and will, therefore depend on the application at hand. A good way to design a filter kernel, that introduces much control over the filter properties, is suggested by Knutsson et al. [25] and is based on a least square formulation of the filter design problem. The main idea is to introduce a number of ideal filters in each domain where they can be easily described. Here the quadrature filter has been defined in the Fourier domain (4.1.1), which is suit-able for this technique. However, it is also desirsuit-able to achieve a filter which is small in the spatial domain, i.e. ideally resembling δ(x), an impulse3. In short, the problem can
be formulated as
$2 ="
j
$Wj(fj$ − Bjf )$2, (4.5)
where f is the desired filter in the spatial domain, Bj is the base matrix that when
multiplied by f will give f described in domain i, and f$
j is the ideal filter for domain
i. Wj is a weight matrix which can be used to make certain parts of the filter more
important, i.e. more sensitive to errors. For this case, the only domain other than the spatial is the Fourier domain. Thus, Bj will here either be the Fourier base matrix or
the unity matrix. In each case respectively, Bjf will give F (ω) and f.
As mentioned above, the phase, φ, of a signal, s(x), is related to the response from these filters. More specifically, phase can be defined as the angle of the filter response between the signal and the filter, derived from
q(x) = (s∗ f)(x), (4.6)
where q is the filter response and ∗ denotes the convolution operator. Assuming s(x) is real valued4, q is complex and can be rewritten on polar form as
q(x) = A(x)eiφ(x). (4.7)
Here A(x) is the magnitude of the filter response, and the phase can now be found through
φ(x) =arg(q(x)). (4.8)
Finally, what is left is to generate a displacement field from the phase of two images. Assuming the two images I1 and I2, we get two different responses from the same set
of filters. Since lognorm filters, from definition 4.1.1, are dependant on the direction of ˆ
n, we actually get one response from each filter in the filter set. A filter set can consist of as many filters wanted, however, there is a minimum of filters needed given a certain dimension of the images to be able to cover the whole Fourier domain. This is due to the fact that the filters have a direction and a limited spread. For example, in two dimensions, at least three filters are needed and in three dimensions, six filters is the minimum to span over the whole domain. This is discussed in more detail for example in [26]. The six directions for 3-D given in [26], for a second order filter, and which are used in this thesis, are
ˆ n1 = c[b, −a, 0]T ˆ n2 = c[b, a, 0]T ˆ n3 = c[0, b, −a]T ˆ n4 = c[0, b, a]T ˆ n5 = c[−a, 0, b]T ˆ n6 = c[a, 0, b]T ,
4The signal in this context is a two or three dimensional image, which in most circumstances is real
4.1 Displacement field estimation 23
a = 2, b = 1 +√5, c =
#
10 + 2√5,
where a, b and c is only for normalisation. From above, it is therefore necessary to apply convolution between the signal and all the n number of filters in the filter set. This leads to n filter responses for each image, each one representing the phase and magnitude in the direction of filter j (see Figure 4.2 for an example), i.e.
qj = (s∗ fj)(x). (4.9)
From this we have one response from I1 and one from I2 for each filter in the filter set:
Figure 4.2: The original X-ray hand (left) along with the phase (centre) and magnitude (right) of it after being filtered with a quadrature filter. Here it is seen that the filter captures edges corresponding to the filter direction, giving more certainty (magnitude) to structures with a frequency within the passband of the filter.
qj(I1) = (I1∗ fj)(x), qj(I2) = (I2∗ fj)(x). (4.10)
Further more, the phase difference at position x between the images can be found through the product between qi, I1 and the conjugate of qi, I2, or
qj(I1)q∗j(I2) = Aj(I1(x))Aj(I2(x))ei(φj(I1)−φj(I2)). (4.11)
From this, the phase difference is, analogue to equation (4.8), equal to
∆φi =arg(qj(I1)q∗j(I2)) = φj(I1)− φj(I2). (4.12)
Since this is the difference between the images, it is also locally proportional to the displacement field estimate at that point:
dj(x)∝ ∆φj. (4.13)
Moreover, what we now have is a set of n displacement field estimates, each one repre-senting the estimate in direction i.
4.1.2 Least square solution
To bring the information from each one of the estimates together into one displacement field estimate, a least square problem is a suitable way of postulating this problem. Since least square problems have the great advantage of being able to weight different points with different importance, as is common in filter design, for example described in [25], this can also be used for the derivation of the total displacement field. In this case, the weight function is called the certainty of the filter responses and is defined as
cj = Aj(I1(x))Aj(I2(x)). (4.14)
The certainty will then be high where the two filter responses match each other in struc-ture, which is sensible, and this makes for all the parts in the least square problem to calculate the whole displacement field, d = (d1, . . . , dN), for N dimensions.
min
d
"
j
[cj(ˆnTjd− dj)]2, (4.15)
for each direction, ˆnj.
Further more, since the Morphon method is an iterative registration scheme, the dis-placement field estimate is only valid for the current iteration, and must be updated with each iteration, thus generating a new field, dk, for iteration k, which brings us to a
central part of the Morphon method.
4.2 Displacement field accumulation
As presented above, each iteration generates a new displacement field, solved from equation (4.15). However, the field for iteration k might not be a good estimate com-pared to estimates prior to that iteration. To prevent bad estimates to disrupt the registra-tion process, a field accumularegistra-tion is therefore used to make good estimates gain greater influence on the final displacement field.
The accumulation also have another effect on the registration process as it prevents interpolation smoothing to be spread throughout the iterations. Since the model (I2here)
has to be deformed in each iteration and a new displacement found from the deformed model, smoothing would occur if dk would be applied directly on the model for each
iteration. Instead, the newly found field is added to the field found and accumulated in iterations before, removing interpolation artefacts since the original model is now used for each iteration. The accumulated displacement field is updated through
d$a= cada+ ck(da+ dk)
ca+ ck
4.3 Displacement field regularisation 25 where da denotes the accumulated field for the current iteration and d$a is the next,
updated accumulated field. caand ckboth describes certainty measures, associated with
the ones mentioned in equation (4.14). ck is directly coupled with the displacement
estimate for each iteration, dk, and is defined as
ck= n
"
j
cj, (4.17)
i.e. the sum of all certainty measures for each quadrature filter. ca, on the other hand, is
another accumulated entity and is a certainty measure weighted with its own certainty, expressed as c$a= c 2 a+ c2k ca+ ck , (4.18) where c$
a, analogue with above, is the updated certainty measure. To summarise; the
accumulation of the displacement field implies that good displacement field estimates will have a higher effect on the accumulated (final) displacement field than poor ones. This creates a robustness of the algorithm as well as increases its reliability.
4.3 Displacement field regularisation
As mentioned in section 4.1.1 about quadrature phase of images, the phase is only a local description of the image. If the accumulated field estimates from equation (4.16) were to be directly applied to the model image, it would most likely tear the image apart, diverging into something far from a successful registration operation. However, since the field represents local changes in the image, it is sensible to think that an averaging of sorts would represent a change on a more global scale. This averaging can be made more or less complicated. Often, it is done by means of a basic Gaussian averaging:
dr = (d∗ g)(x), (4.19)
where g is a Gaussian kernel and dr is the regularised field, but can be made more
ad-vanced in many ways. In this thesis, a special form of normalized convolution, as first introduced by Knutsson and Westin [27], called normalized averaging, is used. Nor-malized averaging introduces the great advantage of allowing to include the certainty measure in the regularisation process as well, adding even more robustness and reliabil-ity to the registration. Normalized averaging is defined as
dr =
((cad)∗ g)(x)
(ca∗ g)(x)
, (4.20)
where the division and multiplication are taken element wise. The effect of a small and a large standard deviation of a Gaussian kernel is presented in Figure 4.3.
Figure 4.3: A displacement field during a Morphon registration process with a small (left) and large (right) regularisation of the field.
4.4 The model
The model in the Morphon method is the image that is deformed by the displacement field, and hence contains the information about the registered object after the registration has finished. However, to perform a successful registration, an appropriate model is needed. The type of model often depends on the context or problem. Nevertheless, a good model should contain approximately the same information regarding shape and intensity variations as the target object. For example, in this thesis segmentation is performed on SPECT images of the heart. A suitable model could therefore reside in a manually segmented heart from another SPECT image.
4.4.1 Deformation
The model is deformed, using a suitable interpolation method, by the accumulated dis-placement field at the end of each iteration and the deformed model is used to generate the new accumulated deformation field, which in turn deforms the original model. This is iterated until some kind of accuracy, δ, or a maximum number of iterations is reached. The accuracy could for example be the change of the accumulated displacement field between two, in the iteration process, adjacent fields, i.e.
δ = $ $ $ $dr, k− ddr, kr, k−1 $ $ $ $ , (4.21)
4.5 Summary 27
4.5 Summary
Here, a simpler version of the Morphon method has been described in some depth. The method takes advantage in its local phase based registration scheme, which is iterated through different resolutions. A displacement field is found for each iteration and is accumulated and regularised to make the registration more robust, as well as adding the option to make the model more or less rigid and deformable.
5
Rotation of Myocardial Perfusion Images
As presented in section 2.3.1, there are some problems with the image analysis in scinti-graphic heart examinations as they are carried out at the time of writing this thesis. The main problems are:
• Interoperator and intraoperator inconsistency, that ends in difficulties in
compar-ing images acquired from different technologists, or even from from the same technologist.
• Time consuming activities for the technologist after images have been acquired. • Difficulties in manual alignment for images with poor uptake.
The approach in this thesis for designing an automatic method for doing all this is based on the Morphon from chapter 4. The automatic process involves some pre-processing, to help the later registration scheme, followed by a 2-D registration and segmentation to decrease the amount of data. After this, the 3-D image should have been decreased to a small volume, containing the left ventricle of the heart. A 3-D Morphon is then applied to the image, from which a displacement field is generated. This displacement field, with enough regularisation, is a rough description of how the model have been rigidly transformed to match the target image. By estimating the angle with the help of two lines, the rotation of the left ventricle can be known, and the target image rotated correctly. The sequence is in this chapter described in detail.
5.1 Pre-processing of images
The first problem in this application that must be dealt with is that the amount of data has to be decreased, in order to increase performance of the analysis. This is primary done in a quite primitive way by cropping the 3-D image along all of its dimensions. The cropping assumes that the technologist has placed the patient correctly in the
γ-camera and that the heart is not positioned in the right side of the body. The cropping is illustrated in Figure 5.1.
Figure 5.1: The pre-cropping that will help the registration process. To the left, the original image, and to the right the cropped image. The left ventricle is marked with a black rectangle in both images.
5.2 2-D Morphon segmentation
After the image has been cropped, the first Morphon based segmentation is performed. As this is a 2-D segmentation, the 3-D volume must first be reduced. This is done by transforming the volume into a max intensity projection (MIP) image, i.e. by generating a 2-D image where each pixel is the maximum of all the voxels that after a projection would result in that pixel. This is performed along the y-axis, which is the axis in the same direction as the nose of the patient, or that is pointing out from the patient’s chest. This is illustrated in Figure 5.2 a).
Figure 5.2: a) Maximum intensity projection along the y-axis (green). b) Maximum intensity projection along the x-axis (red). The resulting projections are symbolised by the purple plane.
5.2 2-D Morphon segmentation 31
5.2.1 Locating the heart
Since other organs can have very similar shape and intensity variation compared to the left ventricle of the heart, it is important to help the registration as much as possible. This is done by finding a point of the heart where the Morphon can be initiated. By assuming that there are only uptake in form of noise in the parts above (in the negative z-direction, which is symbolised by the blue vector in Figure 5.2) the heart, a thresholding can be performed to remove any noise, after which the heart should contain the lowest point in the z-direction. This assumption has been made after investigating a great number of images (over 600). As presented in section 3.1.1, the threshold level can be found in different ways. Here, a dynamic way is needed to keep the thresholding from removing the heart completely. This is done by calculating the histogram of the image. This will result in a curve with some maxima and minima, where the maxima are a representation of intensity levels that are more common in the image. As this do not apply to noise to the same extent, a suitable threshold level that contains the heart can be found. However, this may not be enough if the noise level is very high. In order to separate noise from the heart, the size of the binary clusters that are found is calculated. If the size is too small, it is likely that it does not represent the heart. Additionally, a distance test is performed. This test calculates the distance from the centre of mass (COM) of the image and the point that has been found, where the COM is defined as
%
iI(xi)xi
%
iI(xi)
, ∀i, (5.1)
for a pixel at position xi with intensity level I(xi). As the intensity levels in the MIP
image are much higher in the lower parts of the body, it is there that the COM will be found. The distance from the COM to the heart is therefore relatively short, and points that are found too far away from the COM can be regarded as noise.
5.2.2 Segmentation
When a point of the heart has been found, the 2-D Morphon can be initiated and per-formed. The model used is a MIP image from another heart, which represents the intensity variations in a suitable way. Moreover, as the model may still, despite the tests from above, be initiated quite far from the heart, the registration is first performed on a lower scale of the image (naturally, with a down sampled model as well). Still, the same size of the quadrature filters is used, which implies that the model will "see" longer than on full scale. To perform the registration on different scales is a good way of preventing the registration to converge locally on details, before it has reached the heart. When the model has converged, it should be known where the heart is situated in the 2-D MIP-image. It can then be cut out along the third dimensions, resulting in a
Figure 5.3: The result after two 2-D Morphon segmentations.
cylinder like volume which is oriented along the axis of the maximum intensity projec-tion (y in this case). The resulting volume is then made a MIP along the x-axis, which is illustrated in Figure 5.2 b), and another 2-D Morphon is initiated. However, as much of the interfering organs around the heart has been cut away already, it is sufficient to initiate the model at the COM of the MIP image, which will result in a good registration. Consequently, when the registration has converged, even more can be cut away in the same manner as above. The result is a relatively small box, containing the heart and possibly other organs that lies very close to the heart. An example of the result can be seen in Figure 5.3.
5.3 3-D Morphon segmentation
From the segmentation covered in the last section, the heart is now separated from the majority of interfering organs. The last segmentation process is now made in three di-mensions to be able to gather information about the orientation of the heart. As the co-ordinates of a point that is known to have its position somewhere in the heart, this point is a sensible starting point for the model in the 3-D Morphon registration. How-ever, as the starting point will probably not correspond exactly to the same point in the model, the 3-D registration is here as well performed in multiple resolutions. The
reg-5.4 Angling of the heart 33 ularisation of the displacement field is relatively large on the lower scales to keep the model from deforming incorrectly in the absence of the details of higher resolutions. At full resolution, on the other hand, the standard deviation of the regularisation filter is relatively smaller. This means that the model can deform into more detail, which is nec-essary due to the great variation between images from scintigraphic heart examinations. The model in the D segmentation is as in the 2-D case a manually segmented 3-D heart. The model will deform significantly between different types of hearts. An example of this is illustrated in Figure 5.4.
Figure 5.4: The original model (top-left) and three examples of it after being deformed in the registration process.
5.4 Angling of the heart
After the 3-D Morphon registration has converged, the model should have aligned itself to the left ventricle of the heart. The model can then be used to segment the left ventricle, which is one of the aims of this thesis. Moreover, the rotation of the heart requires some thought on how to represent the rotation of the model.
5.4.1 Angle estimation
To represent an angle between two objects, a minimum of two lines is needed. The first line has been manually defined from the original model. The line does not need to be represented by voxel coordinates in the image, but can be mathematically defined as desired. Here, the line is situated in the middle of the left ventricle, in the cavity called the lumen, and has a length that does not exceed the length of the heart in the same direction. Mathematically, the line can be expressed as
x0+ tv = xy00 z0 + t vvxy vz , t = 0 : l, (5.2)
where x is the starting coordinate of the line, v is the unit direction vector and t is a parameter which determines the length, l, of the line. The deformation field that has been found through the 3-D registration can now be used to deform the line accordingly. However, as a small regularisation has been used on the full scale registration, the line is likely to be deformed into a discrete curve. This is not desirable and hence the field is further regularised. With a large regularisation filter, the global displacement will more and more correspond to a rigid transformation field. However, the deformed line will nonetheless become slightly bent. This is solved by fitting a line to the deformed points of the line, resulting in a linear regression problem:
min
α, β
"
i
[zi− (αxi+ β)]2, x = (x, y), (5.3)
with parameters α and β, or min
˜
x $A˜x − b$, (5.4)
where A is a known basis matrix, ˜x is an unknown parameter vector (here, containing α and β) and b is a known vector containing the function values (here, the z-axis). From linear algebra we know that the solution for a least square problem is:
ˆ
x = (ATA)−1ATb. (5.5)
This will yield a straight line that has been approximately rigidly transformed from the original line, and the angle, θ, for each axis respectively can be calculated from the definition of the scalar product between two vectors, u and v, given as
5.5 Summary 35
5.4.2 Rotation
As the angles of rotation around each axis have now been estimated, the remaining task is to rotate the image heart with the corresponding angles. This can be done in a multitude of ways, e.g. applying a rigid rotation matrix and interpolating the image to the new position. The segmented and rotated heart can then be displayed in any way desired, however, the common clinical way is in slices through each axis. An example of a segmented heart, before and after rotation, can be found in Figure 5.5, where a surface of the volume object is displayed, and in Figure 5.6, where one slice through each axis is shown.
Figure 5.5: Iso-surface plot of an automatically segmented heart (left) and the same heart after rotation (right).
Figure 5.6: Three slices through the rotated heart as they are commonly displayed in clinical use.
5.5 Summary
The segmentation and alignment of myocardial SPECT-images is performed with the aid of two 2-D and one 3-D registrations, based on the Morphon method from chapter 4. The angles are consequently found through deforming a line along with the deformation of the model in the method.
6
Dopamine Transporter Quantification
From section 2.3.2, it is clear that there are currently some problems involved in quan-tification of 3-D scintigraphic brain images which aims at diagnosing e.g. Parkinson’s disease. This relies both in the variability between operators and the fact that much quantification is made in two dimensions, while the acquired images are in fact in 3-D. Thus, The 3-D segmentation method used in chapter 5 has been somewhat modified and complemented, and then applied on the quantification problem at hand.
The algorithm for automatic 3-D quantification works by first dividing the SPECT-brain image into two parts and quantifying each part individually. This is followed by an initi-ation of the 3-D Morphon segmentiniti-ation process. The segmented basal ganglia are then further analysed and quantified with methods comparable to the manual 2-D quantifica-tion commonly used clinically. To be able to do this, a method called weighted principal component analysis (wPCA) is introduced.
6.1 3-D Morphon segmentation
Compared to the segmentation process of the heart in chapter 5, specifically section 5.3, very little pre-processing is necessary in the case of DaTSCAN SPECT images. In myocardial perfusion SPECT images, the image signal from the left ventricle of the heart constitutes a small part of the total image signal. This is in contrast to DaTSCAN SPECT images of the brain, where the basal ganglia (the region of interest) stand out in a much clearer way. This is illustrated in Figure 6.1. Consequently, the initiation of the 3-D segmentation can be greatly simplified compared to the pre-processing in chapter 5; e.g. the two 2-D segmentations from section 5.2.2.
