1. health and fitness
2. disease
3. cancer

Worksheets for the paper:
Quantification of Colon Content from MRI
Using a Semi-automatic Approach
Group 774
December 21st, 2012
School of Medicine and Health
Aalborg University
Worksheets for the article:
Quantification of Colon Content from MRI Using a
Semi-automatic Approach
1. semester master, Biomedical Engineering and Informatics
Group 774
Authors
Astrid Munk
Thomas Holm Sandberg
Søren Sohrt-Petersen
Daniel Vyskovsky
Lin Wang
Supervisors
Lasse Riis Østergaard,
Associate Professor
Jens Brøndum Frøkjær,
Staff Specialist and Head of Research
50 pages
8 copies
Number of Pages:
Number of Printings:
21st of December 2012
The content of this report is freely available, but publication is allowed only with complete reference.
Abstract
Introduction: Constipation is a common problem in the Western world [1]. No golden standard to
assess constipation exists [2]. The aim of this paper was to develop an algorithm for semi-automatic
segmentation of the colon content from magnetic resonance imaging (MRI) in order to estimate
the volume of the colon content.
Methods and Materials: T1 weighted MRI data from liver acquisition with volume acquisition
(LAVA) sequence was used. The initial seed points for region growing were selected using a semiautomatic method based on fuzzy c-means clustering. Growing criteria for a seeded region growing
featured edge detection and intensity thresholds. Total and partial volumes of the segmentation
were calculated. The accuracy was validated by comparison to manual expert segmentations using
the Dice Similarity Coefficient (s) (7.1). The reliability was validated by segmenting continuous
scans.
Results: Five series were segmented with total volumes of 0.8667, 0.0772, 0.6542, 0.8799, and
0.6772 litres, respectively. The comparison with ground truth produced Dice coefficients of s = 0.40
and s = 0.63, respectively. In single coronal slices Dice coefficients of s = 0.60 − 80 were obtained.
The reliability of the algorithm was expressed in a standard deviation of σ = 0.0234 litres from a
mean of µ = 0.5394 litres.
Discussion: The proposed algorithm shows successful and robust results in segmenting bright
colon content. The algorithm proposes a way to semi-automatically quantify volume of colon
content from MRI. Segmentation results show differences between automatic segmentation and
ground truth. Differences are mainly due to the nature of the algorithm being able to segment
content with differences in contrasts between bright colon content and surrounding tissues. Manual
outlining of colon content integrates prior anatomical knowledge. This is not incorporated in
the algorithm that relies on universal segmentation criteria. This results in some parts of the
segmentation leaking to unwanted areas. The empirically selected criteria could be learned by
training the parameters of the criteria. The proposed algorithm may have a clinical application by
helping more objectively assess constipation in patients.
Astrid Munk
Thomas Holm Sandberg
Søren Sohrt-Petersen
Daniel Vyskovsky
Lin Wang
Preface
This is the worksheets describing the development of an algorithm for quantification of colon content
from magnetic resonance imaging (MRI) using seeded region growing. The algorithm utilises acknowledged methods from the fields of image analysis and pattern recognition, and is able to segment colon
content and calculate one total volume and two partial volumes.
The project uses MRI data of human patients acquired from Department of Radiology, Aalborg Hospital. The project group would like to thank Staff Specialist and Head of Research Jens Frøkjær for
providing the MRI data and manually segmenting two series of MRI.
Present project suggests an algorithm for quantification of colon content from magnetic resonance
imaging. In order to provide the reader with requisite knowledge to understand the underlying construct in this project, terms as constipation, diagnostic criteria, and the anatomy of the colon are
described in Chapter 1. Afterwards a review on medical imaging techniques of the colon follows in
Chapter 2. Finally, the aim of this project is set in terms of a problem statement in section 2.3.
The considerations and complications that came up during the development of the proposed algorithm is described in Chapter 3-6. Chapter 3 contains an overview of the image analysis system and
a description of the data set. Chapter 4 deals with an initial segmentation, and Chapter 5 is about
the actual segmentation. In Chapter 6 parameters for clinical use are extracted, while Chapter 7
contains the results of the segmentation and a validation of the algorithm with regard to accuracy and
reliability. Finally, a discussion on the segmentation and other issues is conducted in Chapter 8.
The term colon content is used throughout this project to describe solids in the colon. Colon content
comprises chyme in cecum and colon ascendens, while in the colon transversum and colon descendens
the content is faecal.
v
Contents
Contents
vi
List of Figures
ix
List of Tables
xi
1 Constipation
1.1 Anatomy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 Medical Imaging of the Colon
2.1 Radiography . . . . . . . . .
2.2 Magnetic Resonance Imaging
2.2.1 MRI Parameters . . .
2.2.2 MRI Artifacts . . . . .
2.3 Problem Statement . . . . . .
1
2
.
.
.
.
.
7
7
9
10
11
11
3 Image Analysis System
3.1 Data Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
13
4 Initial Segmentation
4.1 Purpose . . . . . . . . . . . . . . . . . . . . . . .
4.2 Different Methods for Initial Segmentation . . . .
4.2.1 Filtering . . . . . . . . . . . . . . . . . . .
4.2.2 Region Filling . . . . . . . . . . . . . . . .
4.2.3 Edge Detection . . . . . . . . . . . . . . .
4.2.4 Global Thresholding Based on Intensities
4.2.5 Manual Intervention . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
15
15
15
15
15
15
17
21
5 Segmentation
5.1 Purpose . . . . . . . . . . . . . . . . . . . . . .
5.2 Texture analysis . . . . . . . . . . . . . . . . .
5.2.1 Micro-features . . . . . . . . . . . . . .
5.2.2 Co-occurrence matrix . . . . . . . . . .
5.2.3 Difficulties Regarding Texture Analysis
5.3 Region Growing . . . . . . . . . . . . . . . . . .
5.3.1 Common Approaches . . . . . . . . . .
5.4 Region Growing Algorithm . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
23
23
23
23
24
24
24
25
25
6 Post Processing
6.1 Purpose . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Estimation of Cross-sectional Areas . . . . . . . . .
6.2.1 Approach . . . . . . . . . . . . . . . . . . .
6.2.2 Difficulties Regarding Cross-sectional Areas
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
29
29
29
29
32
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
vi
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
CONTENTS
6.3
.
.
.
.
.
.
.
.
.
32
32
33
34
34
35
36
37
37
7 Results and Validation
7.1 Results from Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 Test of Accuracy: Comparison to ground truth using Dice Similarity Coefficients
and a Qualitatively Validation . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Results from Volume Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 Test of Reliability: Comparing Consecutive Scans of one Subject . . . . . . . .
7.2.2 Computational Times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
39
8 Discussion
8.1 Segmentation . . . . . . . . . . . . . . . .
8.2 Automatic Versus Manual Segmentation .
8.3 Sensitivity to Selection of Seed Point . . .
8.4 Optimising the Algorithm for Clinical Use
45
45
46
46
46
6.4
6.5
6.6
Diameter of Inner Joined Spheres . . . . . . . . . . . . . . . . .
6.3.1 Approach . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.2 Difficulties Regarding Diameter of Inner Joined Spheres
Model-based Approach . . . . . . . . . . . . . . . . . . . . . . .
6.4.1 Approach . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.2 Difficulties Regarding the Model Based Approach . . . .
Downgrading the Complexity of the Measurements . . . . . . .
Estimation of Volumes . . . . . . . . . . . . . . . . . . . . . . .
6.6.1 Approach . . . . . . . . . . . . . . . . . . . . . . . . . .
Bibliography
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
39
42
43
44
49
vii
List of Figures
1.1
1.2
1.3
1.4
2.1
2.2
2.3
2.4
2.5
3.1
3.2
3.3
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
5.1
5.2
5.3
The human colon comprises four parts; colon ascendens, colon transversum, colon descendens, and colon sigmoideum. The transitions between the ascending, transverse, and
descending segments are flexura coli dextra and flexura coli sinistra, respectively [3]. . . .
Interpatient variations in shape and location of colon. The variation in distribution, shape,
texture, location and intensities of colon content in between patients are high, and even
within patients, due to time of scan, degree of constipation, and type of food intake [4, 5].
The colon wall [6]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Peristalsis something something. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
3
4
5
A speeding electron hits a tungsten atom. Energy is released in form of a X-ray photon [7]. 7
Abdominal Radiography. The air bubbles in the colon appears as black spots [8]. . . . . .
8
Precession: A rotation of an axis on rotation [9]. . . . . . . . . . . . . . . . . . . . . . . .
9
Left: Gradient coils creating the static fields [10]. Right: RF coil used for abdominal MRI
[11]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Differences in signal intensity between water and fat with varying repetition time, TR [12]. 11
Traditional division of an image analysis system [13]. The components are processed in the
following chapters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Two coronal MR images from a LAVA series. The colon is labelled with green arrows while
artifacts are labelled with red arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Axial representation of the data set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Left: Edge detection using a Sobel kernel. Right: Edge detection using a Gaussian smoothing filter and a Laplacian filter kernel. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Edge detection using a more sophisticated method, the Canny edge detector. . . . . . . .
Left: Coronal slice from a LAVA series. Right: The associated histogram. . . . . . . . . .
Left: Cropping colon content in a coronal slice. The cropped image is labelled with white.
Right: The histogram from the cropped image. . . . . . . . . . . . . . . . . . . . . . . . .
Histogram of one complete series. The lines indicate the cluster centers, and the calculated
threshold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Above: Histogram of original image (coronal slice of data set). Below: Histogram of liver
tissue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Left: True positives from FCM using 7 clusters. Right: True positives from FCM using 5
clusters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Left: Candidate areas for initial seed point selection based on FCM clustering and morphological filtering. Right: Areas of colon content inside the manually drawn region of
interest. These areas are used for final seed points for region growing. . . . . . . . . . . .
Co-occurrence matrix, [14] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4- and 8-connected neighborhoods for region growing [15]. . . . . . . . . . . . . . . . . . .
Crop of the right flexure at the transition from the ascending to the transverse colon. The
contrast between bright colon content and the colon lumen and surrounding tissue is quite
clear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ix
13
14
14
16
17
18
18
20
21
21
22
24
25
26
List of Figures
5.4
5.5
Crop of the resulting gradient filled image after applying a sober filter and morphological
flood-fill operation. The colon content of areas of high gradient peaks results in high
contrast between content and surrounding tissue. . . . . . . . . . . . . . . . . . . . . . . .
Segmentation after region growing has been applied. The pink color shows the initial
seed points segmented based on FCM clustering, and the blue color shows the expanded
segmentation from region growing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1
27
28
Cross-sectional area of a component (yellow). The centroid is marked with a dot, and the
direction of the component is marked with an arrow. . . . . . . . . . . . . . . . . . . . . .
6.2 The direction of each component depends on the components of the eigenvector having the
largest absolute eigenvalue: When the absolute value of the component in the x-direction
was largest, the component was classified into the transverse segment (red). Conversely,
when the absolute value of the component in the y-direction was largest, the component
was classified into either the ascending or descending segments (purple). . . . . . . . . . .
6.3 Segmented colon content. The red lines indicate the plane through each component. . . .
6.4 Axial view of a disk (blue) in a component (not shown). The interslice area is calculated
assuming that the colon does not change in shape in between the slices. . . . . . . . . . .
6.5 Left: Distance map of the segmented colon content. Right: Enlarged segment of the
distance map. Peak points are brightest resulting in a rough skeleton of the shapes. The
values of the distance map (not shown) are radii of inner joined spheres. . . . . . . . . . .
6.6 From left to right: Different estimates of the diameter of the colon content. 1) A reliable
estimate. The diameter of the inner sphere is close to that of the colon content. 2) Less
reliable estimates. The irregular shape limits the size of the diameter of the inner sphere.
3) Not reliable estimates. 4) Less reliable estimates. If this is the situation in one of the
flexures, the reliability might increase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.7 Since the components of the flexures can be composed of both ascending/descending and
transverse part, the majority of the δ points will be located in the flexures. In these cases
the distance map will lack reliability as an estimation of the diameter of the colon content.
6.8 Estimating a rough distribution of the colon using lines created by choosing flexure points.
6.9 Two dividings of colon according to different placements of corner points. . . . . . . . . .
6.10 Problematic position of transverse segment for use of this method. . . . . . . . . . . . . .
34
34
35
36
7.1
7.2
7.3
7.4
7.5
7.6
7.7
7.8
39
40
40
41
41
42
43
44
Manual segmentation in ImageJ indicated by yellow. . . . . . . . . . . . . . . . . . . .
s values in each slice, series 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
s values in each slice, series 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bright intensity terminal ileum segmented as part of colon content . . . . . . . . . . . .
Manual (green) and automatic (red) segmentation resulting in s = 0, 74 . . . . . . . . .
Manual (green) and automatic (red) segmentation resulting in s = 0, 39 . . . . . . . . .
Total volume segmented from data set # 5 . . . . . . . . . . . . . . . . . . . . . . . . .
A leak in the right side of the patient was seen in all 3 series used to test the reliability.
x
.
.
.
.
.
.
.
.
29
30
31
32
33
33
List of Tables
4.1
*) The FCM clustering function in Matlab has a (default) maximum number of iterations.
**) Proportion of false positives to the total number of connected components. . . . . . .
6.1
Percent of total volume in each part of colon. µ is the mean, and σ is the standard deviation. 35
7.1
7.2
7.3
Estimation of volumes. *) Manual segmentation, ground truth. . . . . . . . . . . . . . . .
Test of reliability. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Total computational times for each series. . . . . . . . . . . . . . . . . . . . . . . . . . . .
xi
20
42
43
44
1
Constipation
Constipation is a complex of symptoms of many different pathophysiological processes. Due to a variety of symptoms and risk factors associated with constipation, the etiology might be multifactoral, and
this calls for both a determination of the epidemiology of constipation and a measurement to assess
constipation. Many definitions of constipation exist and this is directly reflected in the prevalence of
constipation. According to a review made by Connel et al the prevalence ranges from 3.4-27.2%. The
wide range is also due to different methods of ascertaining constipation. [16]
The first objective definition of constipation dates back to the 1960s where the bowel frequency was
used as an indicator for constipation. This definition was used for decades until Sander and Drossman
in 1987 suggested that also other symptoms should be taken in consideration. They also found that
physicians and patients had different understandings of constipation. This lead to the definition of the
Rome Criteria with the intention to be able to compare patients. Today the Rome II criteria which
comprise two or more of the following abnormalities for at least three months is used. [2, 17]
The Rome II Criteria:
• Fewer than three bowel movements (BMs) per week
• Hard stool in more than 25% of BMs
• A sense of incomplete evacuation in more than 25% of BMs
• Excessive straining in more than 25% of BMs
• A need for digital manipulation to facilitate evacuation
In the literature three subgroups of pathophysiology regarding constipation exits: Functional constipation (slow-transit constipation), Irritable Bowel Syndrom (IBS), and Outlet Obstruction. Being
able to distinguish between these subgroups often improves the outcome of treatment.
The known etiologies of constipation include mechanical obstruction, metabolic disturbances, neurological disorders, and medication side effects. Furthermore, some patients with constipation lack a
known cause, that is they suffer from idiopathic constipation. A differential diagnosis is Hirschsprungs
disease (Megacolon). Due to a lack of nerve cells in the myentheric plexus in the sigmoid colon,
this part becomes spastic. This results in an accumulation of feces proximally to this point. The diameter could expand to 8-10 cm. [18] A description of the anatomy of the colon is found in Section 1.1.
How is constipation being diagnosed in the clinic? According to World Gastroenterology Organization, no golden standard for assessing constipation exists. However, a recommendation suggests that
history taking including checking the Rome Criteria, physical examination including checking for palpable faeces and gas, and finally diagnostic techniques should be utilised. The latter is often done
by stool analysis or medical imaging of the colon to assess seriousness, or an analysis of colon transit
time depending on the outcome of the history taking and physical analysis.
1
CHAPTER 1. CONSTIPATION
Treatment of constipation is symptomatic. Change of diet and life style and stopping of medication
causing constipation are often the first steps. Further steps are use of fibers and different laxatives
and pelvic floor physiotherapy. [2]
1.1
Anatomy
This section serves to provide the reader with anatomy of the colon (Intestinum Crassum) and its
position and function. The digestive tract is divided into the upper and lower gastrointestinal tracts.
The upper part is composed of the mouth, esophagus, stomach and duodenum, the latter being the
first segment of the small intestine. The lower part is composed of the jejunum and ileum segments of
the small intestine and the entire large intestine. In the small intestine all nutrients from chyme are
digested, and in the large intestine water, sodium and some fat soluble vitamins are absorbed.
The large intestine is composed of the cecum, colon, rectum and anal canal. Hence, the colon is
a segment of the large intestine, and its function is storage of fecal matter. Colon is composed of
four further segments; colon ascendens, colon transversum, colon descendens, and colon sigmoideum.
These segments encircle the small intestine as shown in Figure 1.1. [19]
Fig. 1.1: The human colon comprises four parts; colon ascendens, colon transversum, colon descendens, and
colon sigmoideum. The transitions between the ascending, transverse, and descending segments are flexura coli
dextra and flexura coli sinistra, respectively [3].
Cecum is 6 cm in length and has a diameter of 7 cm. Cecum receives chyme from ostium ileocaele,
from here it continues to colon ascendens. This segment is 0.15 m and has a relatively large diameter,
though smaller than that of cecum. The uppermost part lies close to the lower part of the right
kidney. Flexura coli dextra is the transition from colon ascendens to colon tranversum. This flexure
is often obtuse-angled and is also termed flexura hepatica, which indicates that it lies closely to the
liver. Colon transversum is approximately 0.5 m and is a U-shaped sling with the point of folding
just reaching umbilicus (navel). It can also be V-shaped with the tip reaching the pubic symphysis.
Different manifestations of colon in terms of shape and location is shown in Figure 1.2. Flexura coli
sinistra is the transition from colon transversum to colon descendens colon. This flexure has an acute
angle and is also termed Flexura splenica, which indicates that it is close to the spleen. Colon descendens is 0.25-0.3 m and continues to crista iliaca (wing on the hip bone). This segment is often
contracted with a small diameter. The last segment of the colon is the sigmoid colon. This segment
is 0.4 m in average but can vary from 0.15-0.8 m. Colon sigmoideum has many appearances. It can
wind in various ways in the abdominal cavity.
2 of 50
1.1. ANATOMY
Fig. 1.2: Interpatient variations in shape and location of colon. The variation in distribution, shape, texture,
location and intensities of colon content in between patients are high, and even within patients, due to time of
scan, degree of constipation, and type of food intake [4, 5].
3 of 50
CHAPTER 1. CONSTIPATION
Fig. 1.3: The colon wall [6].
Figure 1.3 shows the physiologic anatomy of the intestinal wall. It comprises four layers. The innermost layer is Tunica mucosa, which epithelial lining amongst others contains Goblet cells. These cells
produce mucus, which provides lubrication for the transportation of chyme and protects the intestinal
mucosa. Tela submucosa is loose connective tissue, storing the submucosal plexus that controls
secretion. Tunica muscularis is smooth muscle characterised by a inner circular layer and three
longitudinal bands, teaniae coli. Because teaniae coli are shorter than the intestinal wall, Haustrae
coli emerges. On the surface of the circular muscle layer is the myentheric plexus, which contributes
to peristalsis. This phenomenon is explained in the next section. The outermost layer of the intestinal
wall is Tunica serosa. [18]
4 of 50
1.1. ANATOMY
Fig. 1.4: Peristalsis something something.
For transportation of the chyme different movements are performed. The peristalsis has already been
mentioned. This is a series of wave-like muscle contractions as shown in Figure 1.4. Some medicaments are known to inhibit the peristalsis resulting in reduced bowel movement which can lead to
constipation [20]. Another important type of movement is mass movement. This movement usually
occurs 1-3 times a day and very often shortly after breakfast has been ingested. In response to a irritated point in colon transversum, a constricted ring appears. 0.2 m distal to this point the haustrae
contract as one unit. This contraction lasts 30 seconds with increasing intensity. After a relaxation
lasting 2-3 minutes, another mass movement will occur. Mixing movements is a category of different
movements. Local intermittent constrictive contractions, which last for 5-30 seconds, are examples
of mixing movements. Due to cocontraction of coli teaniae and the circular muscle in the Tunica
muscularis haustrations proximal and distal to the contraction occur. These haustrations expose the
chyme to the epithelial lining, which enhances the resorption of fluids and dissolved substances. For
a healthy person it takes 8-15 hours for the chyme to pass from ostium ileocaele through the colon. [21]
The next sections are concerned with two methods of how to perform imaging of the colon as a
technique for diagnosing constipation.
5 of 50
2
Medical Imaging of the Colon
2.1
Radiography
A diagnostic technique to assess the degree of constipation is abdominal radiography (AR). Radiography is obtained by X-rays, which are ionising radiation. X-rays are produced by movements of
electrons in atoms. Orbitals are energy levels around an atom’s nucleus, and when an electron goes
to a lower orbital, energy is released as a photon. When this photon collides with another atom,
the energy of the photon is absorbed by the new atom. This forces an electron to a higher energy
level. When the electron falls back to its original energy level, energy is released in form of a photon. Hence, X-rays are electromagnetic energy carried by photons. The energy levels of X-rays bring
shorter wavelengths than visible to the human eye. The equation below describes the relationship
between wavelength λ and photon energy hω:
hω · λ = hc = 1239.842eV · nm
where h is Planck’s constant, h is the reduced Planck’s constant, and c is the speed of light. [22]
Fig. 2.1: A speeding electron hits a tungsten atom. Energy is released in form of a X-ray photon [7].
In practice X-rays are made from a X-ray machine by means of two factors. The first factor is the
voltage difference between a cathode and an anode. By heating the cathode, electrons are discharged.
The anode is made of tungsten and attracts the electrons. When a speeding electron collides with a
tungsten atom, an electron is forced to a lower orbital. As shown in Figure 2.1 As aforementioned,
energy is then released as a photon. When the energy drop is large enough, the photon has a high
energy level making it a X-ray. [23, 24]
When X-rays hit a material a proportion of the X-rays is absorbed, depending on the density of the
material. If placed on a film sensitive to X-rays, the proportion of X-rays passing through the subject
is detected and can be used for a 2D representation of the subject. Since X-rays travel in straight
lines, only small deformations occur in X-ray images. [23, 24]
7
CHAPTER 2. MEDICAL IMAGING OF THE COLON
In the case of assessing constipation from radiography, the amount of colon content is interpreted from
the air bubbles on the images. The fact that colon content is not visible implies some difficulties in
interpreting these images. A study made by Reuchlin-Vroklage et al proved that using plain AR to
assess the amount of colon content is highly subjective and has high interobserver and intraobserver
variability. This impairs the diagnostic accuracy. [25]
To improve the readings from plain AR, different rating tools have been developed, e.g. the Barr,
Leech and Blethyn methods. The rating tool described by Barr et al. evaluates the amount and nature
of content in the ascending, transverse and descending colon, and the rectum. Leech et al. divide the
abdomen into three zones and then evaluate the amount of colon content and bowel dilation in each
zone. Finally Blethyn et al. base their rating tool on the amount of colon content and whether it
is continuous or discontinuous. Moylan et al. evaluated these tools in terms of reliability and found
a poor reliability among multiple examiners [26]. A study made by Bongers et al. agrees that plain
AR has limited value in diagnosing children with functional gastrointestinal disorders [27]. When a
patient is referred, an objective method is required, when the knowledge of ground truth is absent [13].
The biological effects of exposure to ionising radiation are caused by deterministic and stochastic
effects. Deterministic effects occur when a certain threshold has been exceeded. Consequences are
cataract, sterility, or radiation sickness. Stochastic effects occur with no threshold limits but the risk
of an effect occurring is linearly increasing with the dose. Consequences are cancer or hereditary
defects. During an examination the patient is exposed to a radiation dose corresponding to the dose
of background radiation during a couple of days to a couple of years [28]. A principle of radiation
protection is specified as the As Low As Reasonable Possible-principle, which implies that the risks and
benefits of radiography should be compromised. Therefore, a method to avoid unnecessary exposure
should be identified.
Fig. 2.2: Abdominal Radiography. The air bubbles in the colon appears as black spots [8].
8 of 50
2.2. MAGNETIC RESONANCE IMAGING
2.2
Magnetic Resonance Imaging
MRI utilise magnetism and radio frequencies (RF) to generate sectional images of the body. Hence, no
exposure to radiation for the patient. The nucleus of an atom has a mass and an angular momentum
due to its spin. The mass of the nucleus gives rise to a moment of inertia. This inhibits the nucleus
from aligning in a magnetic field. This factor, the ratio of the magnetic moment to the moment of
inertia, is termed the gyromagnetic ratio. Specific gyromagnetic ratios for different isotopes allow a
differentiation during MRI. [23, 24]
MRI is accomplished through various measurements of the rotational precession of the nucleus and
the precession of the magnetic field caused by the nucleus. The precessions are shown in Figure 2.3.
Because the gyromagnetic ratio is fixed for all elements, the precessional frequency is fixed for a given
strength of a magnetic field.
Fig. 2.3: Precession: A rotation of an axis on rotation [9].
The nuclear magnetic resonance signals are produced by elements having nuclear angular momentum.
Hydrogen atoms are of major importance due to their high gyromagnetic ratio. The magnetic resonance emerges from protons alternating from parallel to antiparallel in proportion to a static field.
The alternation is caused by an electromagnetic RF alternating at the Larmor frequency being applied
to the static field. Excessive energy is emitted from protons antiparallel to the static fields, which is
emitted from the protons as electromagnetic radiation.
The precessional frequency is also known as the Larmor frequency. The Larmor equation below expresses what RF will be emitted by a specific element. ω is the rate of precession, γ is the gyromagnetic
ratio, and B0 is the strength of the static magnetic field. When this RF is detected by a radio receiver
turned to this frequency, MRI can occur.
ω = γ · B0
Summed up, applying a Larmor frequency RF and then detecting the RF emissions from the protons,
produce the magnetic resonance signal. Afterwards techniques for digital reconstructions can create
sectional images. [23, 24]
The appropriate pulse sequence and spatial localisation within the subject are obtained by means of
three sets of gradient coils in the MR scanner, one for each of the x, y, and z axes. Radiofrequency
coils serve as both the transmitter and the receiver of the magnetic resonance signals. Preferably these
coils are placed as close to the requested body part, and that is why they are shaped for their purpose.
Figure 2.4 shows the gradient coils and a radiofrequency coil often used for abdominal MRI.
9 of 50
CHAPTER 2. MEDICAL IMAGING OF THE COLON
Fig. 2.4: Left: Gradient coils creating the static fields [10]. Right: RF coil used for abdominal MRI [11].
2.2.1
MRI Parameters
The pulse sequence being used determines the parameters controlling the MRI process. The following
section briefly outlines the parameters starting with three intrinsic parameters followed by four parameters controllable by the operator.
Proton spin density (ρ) is the quantity of resonating spins in a tissue. It determines the magnetic
resonance strength corresponding to the brightness in the image. When an image is proton density
weighted, a long repetition time (TR) between the pulses has been used to detect RF signal from a
larger number of protons. Due to small differences between the percentages of hydrogen in various
tissues, these images have relatively low contrast. Longitudinal relaxation time (T1) is a measure
of the time required for the protons to realign with the external magnetic field to 63% of the maximum
possible strength. The T1 varies significantly for hydrogen in different tissues. The molecular motion
in fat nuclei is rather slow, the T1 is equivalently short resulting in a rapidly regained longitudinal
magnetisation. In water it takes longer to regain the longitudinal magnetisation due to the high mobility of water molecules. The water nuclei are more adversative to give up energy to surrounding
tissues. In a T1 weighted image tissues with short T1 appear bright, e.g. fat and protein-containing
fluid. Transverse relaxation time (T2) is a measure of the time required for the protons to lose
63% of their alignment with each other (spin-spin) after a RF pulse has cause them to align at an angle
to the external magnetic field. When the transverse magnetisation of the nuclei occurs, a transient
MR signal is produced. This signal will decay to zero during T2. This phenomenon is called the free
induction decay (FID). Fat is efficient at exchanging energy resulting in a short T2. The converse is
true of water, which is less efficient at exchanging energy resulting in a long T2. Since, the efficiency
of energy exchange is a major factor, T2 is temperature dependent. In a T2 weighted image water
appears brighter than fat. Repetition time (TR) is the time interval between successive pulses
sequences. Variations in this value have great impact on the characteristics of the image contrast.
Figure 2.5 shows differences in signal intensity with varying TR. Echo time (TE) is the time interval
between a 90 ◦ pulse and the peak of the echo signal in a specific pulse sequence, the spin-echo pulse
sequence. Inversion time (TI) is the time between a 180 ◦ pulse and a 90 ◦ inversion pulse in a so
called inversion recovery pulse sequence. Controlling this factor allows water or fat suppression. Flip
angle α is used in echo imaging. The longitudinal magnetisation is maintained by using a flip angle.
[23, 24]
10 of 50
2.3. PROBLEM STATEMENT
Fig. 2.5: Differences in signal intensity between water and fat with varying repetition time, TR [12].
2.2.2
MRI Artifacts
To secure a high image quality and optimal use of the image, knowledge about artifacts is important.
Image artifacts are structures not normally seen, but they occur at images due to different sources of
limitations. Two major groups of artifacts exist. The first group of artifacts is caused by limitations
of the MRI devices (hardware and software). For instance, the coils can produce artifacts which are
seen as ringing in the images. The second group of artifacts is caused by the physiology of the subject.
This group can be subdivided into motion artifacts and inherent physics. Respiratory artifacts is a
type of motion artifacts impairing the image quality due to long scan times. Motion artifacts can be
reduced by reducing scan time, cardiac gating, or use of a breath-holding technique. An example of
artifacts caused by inherent physics is chemical shifts. Chemical shift artifacts are caused by a slight
difference (parts per million) between resonating frequencies of similar protons. E.g. the resonance
frequency of hydrogen in fat is slightly higher than that of hydrogen in water. This can result in a
bright band next to a dark band at fat-water interfaces. This kind of artifacts can be reduced by
using higher gradient, even though this also can contribute to a higher amount of noise in the image.
[23, 24]
2.3
Problem Statement
From the above sections it is evident that the state of the art to assess the degree of constipation is
inadequate. Plain AR is highly subjective and has a high degree of variability in the readings. In
addition to this the patient is exposed to ionising radiation, which can induce some undesirable effects.
MRI stands out among imaging techniques for soft tissues. In order to assess the degree of constipation
in terms of eliminating subjectivity, an algorithm could be developed to quantify constipation from
11 of 50
CHAPTER 2. MEDICAL IMAGING OF THE COLON
MRI. This gives rise to the problem statement below:
How can an algorithm for automatic segmentation of the colon content from MRI be
developed?
The objective is to develop an image analysis system that comprises an algorithm for semi-automatic
segmentation of the colon content. The output of this image analysis system are parameters that
could support the physician in assessing the degree of constipation.
The scope of the project applies to physicians responsible for the management of assessing constipation. A compromise between a perfect segmentation and a minimum of operator intervention is
assigned a high priority. The limitations in this projects are mainly due to the limited amount of time
and limited amount of data.
12 of 50
3
Image Analysis System
This section serves to create an overview of the algorithm. It consists of an input, three blocks of
image processing, and finally an output. A graphical representation of the overall algorithm is shown
in Figure 3.1.
Fig. 3.1: Traditional division of an image analysis system [13]. The components are processed in the following
chapters.
The input to this image analysis system is a volume data set described in Section 3.1. The first block
is processed in Chapter 4. The segmentation block is processed in Chapter 5. The output of the image
analysis system is a set of parameters that are calculated in the rightmost component in Figure 3.1.
This process is described in Chapter 6. The image analysis system is tested in Chapter 7.
3.1
Data Set
Five sets of MRI were supplied acquired at the Department of Radiology, Aalborg Hospital. The
images were provided in the Digital Imaging and Communications in Medicine (DICOM) format.
Conformance to this standard claims amongst more a protocol to be followed for network communications, including syntax and semantics when using this protocol, and also a set of media storage
services to be followed for instance file format and medical directory structure.
Different series were investigated in order to visually extract features for the segmentation. All series
were from a GE Medical System scanner using a magnetic field strength of 1.5 T. Series with Liver
Acquisition with Volume Acquisition (LAVA) sequence applied were chosen. This sequence is derived
from a 3D spoiled gradient echo pulse sequence. The gradient echo is achieved by dephasing and then
rephasing the spins with negatively and positively pulsed gradients, respectively, creating the echo.
The first RF pulse in the sequence tilts the magnetization by a flip angle α. The contrast in the image
is controlled by α. For the data sets, TR = 3.784-3.852 ms, TE = 1.776-1.804 ms, TI = 7 ms, α =
12 ◦ . The rather small value of α indicates that the longitudinal magnetization is slightly decreased
resulting in a faster recovery time. This allows for the short TR. As shown in the Figure 2.5 this short
TR time results in a preferable contrast.
Additional, prospective scans of two subjects were performed using a field strength of 3 T. For these
scans, TR = 4.508-4.512 ms, TE = 2.1128-2.128 ms, TI = 5 ms, α = 12 ◦ . These scans were used to
test for reliability.
Short TE and TR enhances T1 in the data set. The colon content is fatty and high in protein, which
results in this tissue appearing bright compared to surrounding tissue. This indicates that an untrained eye could distinguish colon content from other tissues explicitly, which appears to advantage
13
CHAPTER 3. IMAGE ANALYSIS SYSTEM
when the segmentation algorithm is to be continuously validated.
Figure 3.2 shows two images from a LAVA series:
Fig. 3.2: Two coronal MR images from a LAVA series. The colon is labelled with green arrows while artifacts
are labelled with red arrows.
The images in Figure 3.2 demonstrate more features:
• Colon content is distinct/high in intensity in images from a LAVA series
• Colon has a wall and the edge appears quite clear
• The colon is apparent in more slices
• Artifacts are as prominent as colon content which necessitates a technique for removal of artifacts
Figure 3.3 shows an outline of an axial view of the data set.
Fig. 3.3: Axial representation of the data set.
The terms slice thickness (ST) and spacing between slices (SBS) are explained in figure 6.8. The total
volume (TV) is dependent of the number of slices (NS), by:
TV = (ST + SBS) · NS
For the given data set, ST = 3.60 mm and SBS = 1.80 mm.
14 of 50
4
Initial Segmentation
4.1
Purpose
To select and extract features for the further segmentation. The selection and extraction will result
in an initial segmentation. Based on unfortunate features in the data set, methods for suppression of
or removal of artifacts are desired.
4.2
Different Methods for Initial Segmentation
Many ways of extracting features from an image exist [13]. The given data set and applications are
decisive factors. In broad outlines the features can be based on pixel intensities or on calculations of
these. The selection and extraction of features are emphasised since the purpose is to maximise tissue
contrast differentiation. The following section contains considerations regarding selection of method.
4.2.1
Filtering
To filter out noisy single high-intensity pixels and to be able to retrieve statistic measurements on
single connected components for later determination of region growing criteria, all small connected
components (< 50 voxels) were removed.
4.2.2
Region Filling
The symptoms of constipation are due to increased content in the colon. The factors contributing to
the symptoms are both solids and air bubbles. Hence, the segmentation should include both solids
and air bubbles inside solids. In the implementation air bubbles were filled using a flood fill algorithm
in Matlab.
4.2.3
Edge Detection
As stated in section 3.1 the colon content has three distinct features. One feature is concerned with
the colon having a wall. This implies that edge detection could be useful as a method for initial segmentation. An edge in an image is considered as a gradient within the image. Hence, edge detection
extracts the features of an image where an abrupt change in the intensities is present. To detect these
changes edge detection makes use of differential operators. Different edge-detector filter kernels have
been proposed to approximate derivatives of the image gradients. However, using edge-based features
can result in incorrect edges due to noise, or in an over- or under-segmentation.
The following Figure 4.1 shows an image from the data set with different kernels applied.
15
CHAPTER 4. INITIAL SEGMENTATION
Fig. 4.1: Left: Edge detection using a Sobel kernel. Right: Edge detection using a Gaussian smoothing filter
and a Laplacian filter kernel.
The image to the left in Figure 4.1 uses the Sobel kernel, which performs a 2D spatial gradient measurement, the first spatial derivative. The Sobel kernel consists of a pair of 3x3 kernels. The kernels
are identical except that the second is rotated by 90 ◦ . This design assures that the edges horizontally
and vertically relative to the pixel grid are measured. The Sobel kernel is not as sensitive to noise as
other first derivative edge-detectors, since it differentiates in one direction and implements a Gaussian
average filter in the other direction. This reduces the impact from noisy pixel. A first order derivate
filter gives the strength of a gradient but does not provide the precise onset of the edge. The latter
can be obtained by using a second order derivative filter. The result of using a Sobel kernel is a broad
but not continuously outline of the colon content.
The image to the right in figure 4.1 uses the Laplacian of Gaussian method. The Laplacian is a measure of the second spatial derivative of gradients. A second order derivative filter identifies the precise
onset of the edge, but leaves only little information on the strength of the edge and also amplifies
noise. Since this method is very sensitive to noise, a Gaussian smoothing filter has been applied to
suppress noise in the image prior to computing the Laplacian, as the title also indicates. It is evident
that the result of using the Laplacian of Gaussian is an over-segmentation.
The Canny edge detector is acknowledged as one of the best edge detection methods. This method
includes a smoothing by use of a Gaussian kernel, calculations of edge strengths and directions,
suppression of pixel not identified as edge pixels, and finally thresholding the remaining pixels. Figure
4.2 shows the result of applying the Canny edge detector. Despite of acknowledgment of the method,
the result includes too many false positives to base an initial segmentation of this method.
16 of 50
4.2. DIFFERENT METHODS FOR INITIAL SEGMENTATION
Fig. 4.2: Edge detection using a more sophisticated method, the Canny edge detector.
Edge detection as a method for initial segmentation was not satisfying, mainly because no clear border
between colon content and surrounding tissue is apparent. Even though three very different kernels
were applied, indications of an over-segmentation were obtained.
4.2.4
Global Thresholding Based on Intensities
Another way of extracting features is intensity thresholding. The basic idea is that pixels having a
value above a given threshold T are assigned to one region. The output b(x,y) from thresholding an
intensity image I(x,y) is a binary image:
(
b(x, y) =
1
0
if I(x, y) > T
otherwise
This method is really applicable when the input image is divisible into a foreground of objects of the
same intensity and and a background having intensities distinct to the foreground. I.e. this method
is limited to rather simple cases.
The threshold can be estimated in many ways, the most simple being empirically. Operator selection
of thresholds has been reported in a study [13]. However, in the task of segmenting the colon content, the need of manual intervention should be minimised in order to stress the physician minimally.
In addition to that inter- and intra variation might influence the selected threshold. An automatic
method for choosing the thresholding is necessitated.
To identify a proper threshold for the data set described in Section 3.1, histograms were made. An
image and its associated histogram from the LAVA series are shown in Figure 4.3:
17 of 50
CHAPTER 4. INITIAL SEGMENTATION
Fig. 4.3: Left: Coronal slice from a LAVA series. Right: The associated histogram.
To extract features characteristic for the colon content, a cropping was done. The result is shown in
figure 4.4:
Fig. 4.4: Left: Cropping colon content in a coronal slice. The cropped image is labelled with white. Right:
The histogram from the cropped image.
Two histograms in Figure 4.3 and Figure 4.4, respectively, indicate that colon content has higher intensities compared to surrounding tissue. Please note that the resolutions on the abscissas are different.
It is assumed that colon content covers the same intensity spectrum for all slices from the same series.
This property yields that globally thresholding based on intensities could result in a valid initial segmentation. Which threshold intensity to use is the next question. Since no training data was used,
an unsupervised method for classification was required.
Unsupervised classification is often termed as cluster analysis. In this analysis a search for similarity
in the data set is conducted in order to determine whether so-called clusters can be formed. Fuzzy cmeans (FCM) is a method of clustering which allows one piece of data to show a degree of membership
to two or more clusters.
18 of 50
4.2. DIFFERENT METHODS FOR INITIAL SEGMENTATION
Clustering Using Fuzzy c-means
To identify a global threshold, clustering using FCM was attempted. This was done to group colon
content in one cluster in order to conduct an initial segmentation. Clustering is a mathematical tool
for achieving natural groupings in a data set. The objects inside each cluster show a certain degree of
similarity. The algorithm for this method is as follows:
1. Choose a number of clusters and place initial cluster centers in the feature space
2. Measure distances from all samples to the centers resulting in an initial classification
3. Compute mean values of the clusters to find new cluster centers
The above algorithm is repeated until the algorithm converges meaning that the new cluster centers
are equal to or very close to the former cluster centers.
In Matlab the function fcm was used to conduct clustering using FCM. This which gave rise to the
following issues:
1. What is the optimal number of clusters?
2. Which are the optimal initial cluster centers?
3. Which cluster is used as colon content/true positives?
Ad 1: The optimal number of classes is a compromise on how many false positives to exclude, how
many true positives to include.
Ad 2: The FCM clustering function in Matlab suggests random values for initial cluster centers. Being
able to choose the initial cluster centers in advance might reduce the number of iterations.
Ad 3: To address the aforementioned issues, different numbers of clusters were used on the data set.
As shown in Figure 4.4 the colon content has high intensity between 0.3 and 0.6, which is among the
highest intensities in the image. This indicates that the clusters with the highest intensities should
be identified as colon content/true positives. The threshold is calculated as the average of the highest
cluster center and the second highest cluster center, as illustrated in figure 4.5.
19 of 50
CHAPTER 4. INITIAL SEGMENTATION
Fig. 4.5: Histogram of one complete series. The lines indicate the cluster centers, and the calculated threshold.
Table 4.1 contains the results from this review.
No of clusters
2
3
4
5
6
7
Threshold intensity
0.1335
0.3165
0.3601
0.4023
0.4805
0.5292
No of iterations
29
52
73
100*
100*
100*
Proportion**
85/85
76/88
51/84
37/67
39/68
26/54
Table 4.1: *) The FCM clustering function in Matlab has a (default) maximum number of iterations. **)
Proportion of false positives to the total number of connected components.
As previously mentioned in this section most of the colon content has intensities higher than 0.3. As
shown in Table 4.1 this excludes the threshold obtained using two clusters. From segmentations using
3 or 4 clusters, some parts of the liver tissue were included. These are considered as false positives. To
identify the number of clusters resulting in a satisfying compromise between true positives and false
positives, a histogram of the liver was made. In figure 4.6 histograms of both original image (coronal
slice from data set) and cropped image of liver are shown. When comparing the latter histogram
with the histogram of colon content in figure 4.4, it is noticed that colon content is slightly higher in
intensity.
20 of 50
4.2. DIFFERENT METHODS FOR INITIAL SEGMENTATION
Fig. 4.6: Above: Histogram of original image (coronal slice of data set). Below: Histogram of liver tissue.
Based on this a FCM using 5, 6, or 7 clusters could be used. It could be desirable to choose 7 clusters
since this one had the fewest false positives.
Fig. 4.7: Left: True positives from FCM using 7 clusters. Right: True positives from FCM using 5 clusters.
However, as shown in figure 4.7, choosing 7 clusters only finds true positives in some parts of the
present colon content. Use of 6 clusters gave similar results. Since a FCM using 5 clusters resulted in
true positives in more components of the colon content, this approach was chosen. However, because
the artifacts and the breast tissue demonstrate intensities as high, or even higher than colon content,
this approach still included a number of false positives. This fact is considered in the next section.
4.2.5
Manual Intervention
In order to narrow down the initial segmented areas to only contain true positives, i.e. connected
components representing colon content, regions of interests were drawn on several images of the series
with the initial segmentation. As shown in Figure 4.8 to the left, FCM clustering and morphological
filtering include both true positives and false positives. Depending on the quality of the manual
segmentation, the resulting areas represent only components with true colon content. Regardless of
21 of 50
CHAPTER 4. INITIAL SEGMENTATION
the quality of this intervention, this step composes the foundation for the region growing described in
chapter 5.
Fig. 4.8: Left: Candidate areas for initial seed point selection based on FCM clustering and morphological
filtering. Right: Areas of colon content inside the manually drawn region of interest. These areas are used for
final seed points for region growing.
22 of 50
5
Segmentation
5.1
Purpose
The purpose of segmentation is to divide an image into desired regions. These regions are then targets
for a classification.
5.2
Texture analysis
Texture can be described as an area, where the variety of intensity has specific features, that is
perceived as uniform. Textured region comprises of a set of pixels, satisfying a given gray-level
property, which is repeated in an image. They are frequent among others in medical images and well
recognised by human observers. Textures help in recognising borders, types of structures, shapes of
objects and their spatial placement [15]. Textures can be divided into categories according to:
• periodicity – strong textures with regular spacing, or rather stochastic textures without expressed
spatial regularity (weak textures)
• size – large or small primitives (groups of pixels), determine whether the texture is perceived as
fine or coarse
• orientation of primitives – isotropic (same in all directions) or directional (with one leading
direction)
5.2.1
Micro-features
Using this method it is assumed, that the texture contains periodic transitions of image brightness
(called micro-features). This means that sets of points have a specific intensity distribution in the
space domain. The extraction of a local feature is performed by a 2D convolution of the image with
a filter mask, that represents a model of the micro-feature pattern. These masks are empirically set
matrixes, usually of a 3x3 size. Nine masks were tested on our image data, here is an example of 2 of
these masks: [1 2 1; 2 4 2; 1 2 1]/36 and [1 0 -1; 2 0 -2; 1 0 -1]/12.
The segmentation method comprised of the following steps. The whole image was convolved with the
nine mentioned masks. The result of convolution of reference pixel (center of a single component after
FCM clustering) was saved in a separate variable. The convolved image was then scanned pixel by
pixel and using the euclidean distance formula:
v
uN
uX
D(a, b) = t (ai − bi )2
i=1
a difference of each pixel from reference pixel (both after convolution) was calculated for each mask.
The distance for all masks in each pixel was then averaged. Finally a simple threshold value (a) was
set. Every pixel with a distance < a was concerned as a member of the texture. The problem arose
with the final step – setting a threshold. The liver appeared to have very similar texture properties
and was very often included as undesired false positive result. Setting a small threshold resulted in
23
CHAPTER 5. SEGMENTATION
segmenting small areas around the center of each component, with missing areas while approaching
the boundary of components (drops in intensity values). This could partially be solved by setting a
multilayer perceptron artificial neural network as a classifier, as proposed in [29].
5.2.2
Co-occurrence matrix
A gray level co-occurrence matrix (GLCM) is a tabulation of frequency of different pixel brightness
value combinations occurring in an image. It is a n x n matrix P , where n represents the number
of gray levels within an image. In a GLCM texture, the relation between two pixels (reference and
neighbor) at a time is considered. P [i, j] accumulates the pixel pairs having intensity values i and j.
Pixel pairs in the matrix are defined by a displacement vector d = (dx, dy), where dx is the number
of pixels moved in horizontal and dy a number of pixels moved in vertical direction.
Fig. 5.1: Co-occurrence matrix, [14]
For quantification purposes various textural features can be calculated, including energy, variance,
correlation etc. Since the measures consider relations between groups of two pixels (in the original
image), the texture calculations are of second order. The main problem using this method, considering
the number of gray scale levels in a DICOM image, is computational time. The number of levels can of
course be reduced, resulting commonly in a 32x32 or 64x64 co-occurrence matrix. Further, for different
numbers of displacement vectors, this matrix grows in third dimension. In our case the following
displacement vectors, forming 4-connected pixels, were used: d1 = (1, 0), d1 = (0, 1), d1 = (−1, 0)
and d1 = (0, −1). This quantisation leads in final to loss of valuable details in image and many false
positives especially among the liver area and parts of the neighboring small intestine.
5.2.3
Difficulties Regarding Texture Analysis
For the co-occurrence matrix a limitation presents the computational time, depending on the predefined number of levels and displacement vectors. The same accounted for micro-features based on the
number of chosen matrixes to convolve with the original image.
The main difficulty in applying texture analysis in these series was with redefining and setting parameters for including voxel as a part of the texture. In other words, it would be necessary to implement
a neural network classifier, or similar complex sorting algorithm to correctly classify voxels in the
image. A basic approach by applying simple thresholds resulted in varying results of over- or undersegmentation. Additional user interaction would be necessary for inputting and fine tuning thresholds
for single image series, which is acceptable in technical practice, however not in clinical application.
Texture analysis could provide satisfying results when applied and adjusted on a series by series basis
which does not conform with the aim of the project.
5.3
Region Growing
Region growing falls into region-based segmentation methods. That is segmentation based on areas
of certain common property, using the concept of homogeneity. The principle of region growing is
24 of 50
5.4. REGION GROWING ALGORITHM
planting a seed in a potential region. The seed has properties typical for that region, usually intensity
or local mean. The surrounding pixels are aggregated to the region when fulfilling certain criteria of
homogeneity. These criteria are termed as growing criteria.
This iterative approach is then followed: For every pixel already recognized as a member of region
(initially the seed), check neighboring pixels (in 2D, 4- or 8- connected) and compare parameters of
these candidates to parameters of the seed or other set condition. If the conditions are met, add the
current pixel to region, if not, mark pixel as not accepted. Algorithm ends, when there are no pixels
to be added – for all pixels belonging to region, their neighbors are either members of the region or
marked as not accepted [15].
Fig. 5.2: 4- and 8-connected neighborhoods for region growing [15].
5.3.1
Common Approaches
Static Criterion
New pixel (pj ) is added to region after comparing with original seed pixel (ps ) in a predefined value
(tolerance) T according to formula: |ps − pj | ≤ T . The initial selected seed largely influences resulting
region. This method is not robust against noise. Another approach to initialize the region is to use
multiple seeds. This is either done by calculating local statistics in a mask around an initial seed, for
instance mean value or standard deviation, or by planting more initial seeds in an intended region to
better span the variance.
Dynamic Criterion
New pixel is added to region after comparing with previously added pixel (pi ), under condition that
|pi − pj | ≤ T . Another option is a comparison to current mean value of region |pj − pa | ≤ T . This
helps to grow on a slowly changing area. The order in which pixels are added to region affect the final
region [15].
5.4
Region Growing Algorithm
The main aim of the region algorithm in this project is to segment all desired areas of colon content
in the image series. This works ideally by letting initial seed points grow and extend their areas to
include every single pixel which is considered as colon content and not include any unwanted non-colon
content pixels. This should ideally produce a segmentation where every single pixel belongs to either
the desired colon content (true) or non-colon content (false). This way the region growing produces
a binary image where all pixels either are true or false e.g. colon content or not. To achieve this,
two conditions must be true: a seed point is placed in every colon content component of the image
and a ideal criterion that makes the region growing include all true colon content pixels and exclude
all non-colon content pixels. This then requires two things: correct placement of seeds and correct
growing criteria. As earlier described, the initial segmentation for selection of seed points relies on
FCM clustering and drawing of manual regions of interest. This ideally produces seed points in all
components of colon content, which then requires the selection of criteria for the region growing to
25 of 50
CHAPTER 5. SEGMENTATION
expand the areas to the desired segmentation.
A common criterion used for region growing is based on the intensities of the image and thereby the
contrast between wanted and unwanted areas. This feature is usable as a criterion for segmentation
of the colon content while the intensity of the colon content generally is brighter than that of colon
lumen and surrounding tissue as seen in Figure 5.3.
Fig. 5.3: Crop of the right flexure at the transition from the ascending to the transverse colon. The contrast
between bright colon content and the colon lumen and surrounding tissue is quite clear.
Two approaches of using the intensity as a criterion for growing were tried out. The initial segmentation
results in a number of connected components that all are possible colon content. In the first approach
the pixel with maximum intensity of each component was used as seed point. The threshold for
inclusion was calculated as the intensity mean of a 3x3x3 mask around the seed point. This would
label each new potential pixel I as either true (colon content) or false (non-colon content) by the
following:
(
I=
true
f alse
if Tlow < I < Thigh
otherwise
where Tlow and Thigh are the lower and upper threshold that the potential new pixel has to satisfy
before being added to the region and is calculated as:
Tlow = µs,3x3x3 − γ
(5.1)
Thigh = µs,3x3x3 + γ
(5.2)
where µs,3x3x3 is the mean of a 3x3x3 mask around the original seed pixel and γ is an empirical chosen
value e.g. 0.2.
This approach however did not take into account the already segmented pixels in the components
from the initial segmentation and the empirical chose value γ did not take the span of pixel intensities
of the colon content into account. Therefore, the second approach used all of the initial segmented
pixels as potential seed points, but checked whether surrounding pixels had already been marked as
colon content from the initial segmentation. Also, the threshold used in this approach were based on
statistics of the initial segmented components. This resulted in the following threshold Tlow and Thigh :
Tlow = µcc − α · σcc
(5.3)
Thigh = µcc + α · σcc
(5.4)
26 of 50
5.4. REGION GROWING ALGORITHM
where µcc is the mean intensity of the connected component used as seeds, α is a empirically chosen
growing parameter, σcc is the standard deviation of the connected component.
This resulted in an improved segmentation by tuning of the α parameter and takes into account the
variation in pixel intensities of the colon content initially segmented.
However, the region growing would still tend to leak in areas when the contrast between colon content
and surrounding tissue was low. Using a too low value of α would result in under-segmentation and
a too high value would result in leakage and over-segmentation. Therefore a second growing criterion
was added to the region growing to give better control of the growing.
To retrieve the second growing criterion, edge detection based on sobel filters in horizontal and vertical
direction was performed on the original image series. They were respectively used to compute the
two-dimensional convolution as estimates of first order derivatives on the x- and y-axis. These two first
order derivatives were further combined into a gradient magnitude. Colon content having a brighter
intensity than surrounding tissue tends to have a surrounding dark area formed by the colonic wall,
making a step in intensity that could be captured using the mentioned edge detection method. Forming
a gradient image with edges allowed performing of a morphological flood-fill operation that fills the
surrounded areas of colon content. The resulting gradient filled image can be seen in Figure 5.4.
Fig. 5.4: Crop of the resulting gradient filled image after applying a sober filter and morphological flood-fill
operation. The colon content of areas of high gradient peaks results in high contrast between content and
surrounding tissue.
This laid the basis for the second region growing criterion. For every new potential pixel of a region,
the mean of a 5x5x5 mask on the pixel position was calculated from the gradient filled image which
was compared to the mean of the all pixels in the connected component of the seed points. This
resulted in the following second growing criterion:
(
I=
true
f alse
if µgradient,5x5x5 > β · µgradient,cc
otherwise
where I is the new potential pixel under investigation, µgradient,5x5x5 is the mean of a 5x5x5 mask on
the gradient image at the pixel position currently being investigated, β is a empirically chosen growing
parameter, and µgradient,cc is the mean of the connected component from the gradient image.
The resulting segmentation using both intensity and gradient features of the image results in a better segmentation with less vulnerability to surrounding low contrast tissue and hence less leakage.
Figure 5.5 shows the result from applying the region growing algorithm.
27 of 50
CHAPTER 5. SEGMENTATION
Fig. 5.5: Segmentation after region growing has been applied. The pink color shows the initial seed points
segmented based on FCM clustering, and the blue color shows the expanded segmentation from region growing.
28 of 50
6
Post Processing
6.1
Purpose
The purpose of the post processing is to add clinical value to the segmentation. By providing the
physician with measures, he or she can relate to a degree of constipation, the level of subjectivity in
assessing constipation might decrease. The proposed algorithm presents an objective measurement of
the colon content.
To add clinical value to the segmentation, some measures should be extracted from the segmentation.
These measures should be objective and allow for a comparison of scans in time and performing of
follow-ups of treatments.
6.2
Estimation of Cross-sectional Areas
The degree of constipation is related to the cross-sectional area of the colon content, given that is a
measure of how loaded the colon is.
Fig. 6.1: Cross-sectional area of a component (yellow). The centroid is marked with a dot, and the direction
of the component is marked with an arrow.
6.2.1
Approach
To obtain an estimate of the cross-sectional area it was assumed that colon content was tubular shaped
or like a pearl necklace consisting of many small components. With this assumption, the largest crosssectional area for each component would be located in the centroid of each component. The centroid
was provided by a Matlab function regionprops and its method Centroid.
To be able to calculate the correct cross-sectional area, the direction of each component of the colon
content was found. This was done by calculating a convex hull of each component was calculated
by means of a Delaunay Triangulation (DT). In 3D the DT is a triangulation of every 3-simplexes
created from points within the component which circumsphere does not contain any other points in
the component. From the DT a principal component analysis (PCA) was conducted. The PCA is an
orthogonal transformation converting a set of observations into linearly uncorrelated variables. This
transformation is based on calculated pixel intensities. PCA uses linear correlation and identifies a
connection between samples and a linear trend line. The variable having the largest variance is assigned as the first principal component. The second principal component has the highest variance
29
CHAPTER 6. POST PROCESSING
possible under the constraint to be orthogonal to the preceding component, etc. Hence, this analysis
transforms data into a new coordinate system wherein the first coordinate is held by the component
having the largest variance.
Hence, the principal components are obtained by computing a covariance matrix for the coordinates
that compose the DT of each component, and then finding the eigenvectors and eigenvalues from this
matrix. The eigenvector having the largest eigenvalue has the direction of the component, assuming
that the largest variation is the largest distance.
To retrieve the location of the content, the eigenvector was evaluated. Depending on the eigenvector
having the largest absolute eigenvalue, each component was classified into belonging to either the
descending/ascending or the transverse colon. When the eigenvalue was larger in the x-direction,
the component was classified into the transverse segment, and when the eigenvalue was larger in the
y-direction, the component was classified into either the ascending or descending segment. In case
having the larger eigenvalue in the z-direction, the second largest eigenvalue was used to indicate of
the direction of the component. Figure 6.2 summarises this.
Fig. 6.2: The direction of each component depends on the components of the eigenvector having the largest
absolute eigenvalue: When the absolute value of the component in the x-direction was largest, the component
was classified into the transverse segment (red). Conversely, when the absolute value of the component in the
y-direction was largest, the component was classified into either the ascending or descending segments (purple).
The cross-sectional area was calculated as the plane through the centroid orthogonal to the direction,
i.e. the eigenvector.


a


N =E= b 
c
The coordinates for the centroid are cf. item 1:


x0


center =  y0 
z0
30 of 50
6.2. ESTIMATION OF CROSS-SECTIONAL AREAS
These properties are utilised in the following function for a plane:
plane(x, y, z) = b · (x − x0 ) + a · (y − y0 ) + c · (z − z0 )
The coordinates for each component were evaluated in the equation for the plane. The sets of coordinates resulting in the equation being approximately equal to zero are coordinates on the disk.
Fig. 6.3: Segmented colon content. The red lines indicate the plane through each component.
The area of this plane is one measure of how constipated a patient is. To compute the area the slice
thickness, spacing between slices and pixel spacing are read from the dicomHeader.
In general the area of the plane is calculated as the area in the slices added to the area between slices:
oneSliceArea = areaInSlice + areaBetweenSlices
The area in slice is computed from the slice thickness and the length of the plane in the current slice.
The area between slices is estimated by the approach shown in Figure 6.8. The length of the disks
in the surrounding slices, slice(n) and slice(n+1) are multiplied by half the spacing between slices,
respectively, and then summed. Using this approach it is assumed that the colon does not change in
shape inbetween the slices. This is a fair assumption since the spacing between slices is rather small
compared to the slice thickness.
31 of 50
CHAPTER 6. POST PROCESSING
Fig. 6.4: Axial view of a disk (blue) in a component (not shown). The interslice area is calculated assuming
that the colon does not change in shape in between the slices.
The area is then calculated as the length of the plane in each slice (planeLength) multiplied by the
pixel spacing, this is then multiplied with the slice thickness (ST) multiplied by the spacing between
slices (SBS):
totalArea =(planeLength(1) · pixelSpacing) · (ST + SBS · 0.5)
+ planeLength(n) · pixelSpacing) · (ST + SBS · 0.5)
+
n−1
X
planeLength(i) · pixelSpacing) · (ST + SBS)
i=2
with n being the number of slices the area appears in.
6.2.2
Difficulties Regarding Cross-sectional Areas
Since the degree of constipation depends on the volume of the colon content, mainly, initial attempts
on adding clinical value to the segmentation were based on geometrical parameters. However, the
segmentation resulted in components having more than one direction, e.g. belonging to both the
ascending and transverse part. From this incorrect direction vectors would be calculated which gives
rise to incorrect cross-sectional areas.
6.3
Diameter of Inner Joined Spheres
The diameter of the colon content could be a measure of the degree of constipation. The shape of
the colon content can range from simple and circular to very complex and not ideal in any way. That
complicates the task of estimating the diameter. However by assuming that the content in some degree
is tubular it could be possible to estimate a diameter.
6.3.1
Approach
A distance map provides the distances from all pixels (ones) in one component to the nearest zero
pixel. The peak points of this distance map would then compose radii of inner joined spheres. The
amount of peak points in a component depends on the complexity of the object. The more irregular
the more peak points are expected. Figure 6.5 below shows a distance map of the segmentation from
one series.
32 of 50
6.3. DIAMETER OF INNER JOINED SPHERES
Fig. 6.5: Left: Distance map of the segmented colon content. Right: Enlarged segment of the distance map.
Peak points are brightest resulting in a rough skeleton of the shapes. The values of the distance map (not
shown) are radii of inner joined spheres.
6.3.2
Difficulties Regarding Diameter of Inner Joined Spheres
The quality of this estimate truly depends on the shape of the colon content. Colon content is tubular
to some degree. The smooth muscle layer in the colon wall, Tunica muscularis, has an inner circular
layer. The surface of this layer is the myentheric plexus, which contributes to peristalsis, as described
in Section 1.1. This feature of the colon wall could contribute to the colon being predominantly circular. However, the segmented volumes can be far from circular in that the concavities are dominant.
This is shown in Figure 7.5, where the algorithm nicely segments the bright colon content.
This method tends to produce small values as estimates of the diameter, especially when the shape is
complex. The larger the value, the more reliable the estimate is. The Figure 6.6 shows four examples
of shapes of colon content and associated inner spheres. It is seen that the more circular the shape of
the colon content is, the more reliable the estimate is.
Fig. 6.6: From left to right: Different estimates of the diameter of the colon content. 1) A reliable estimate.
The diameter of the inner sphere is close to that of the colon content. 2) Less reliable estimates. The irregular
shape limits the size of the diameter of the inner sphere. 3) Not reliable estimates. 4) Less reliable estimates.
If this is the situation in one of the flexures, the reliability might increase.
In terms of adding clinical value to the algorithm, however, the δ largest values from distance maps
of the entire segmentation could indicate a distribution of the colon content. In some cases every of
the δ points will be located in the flexures. This will not add any information since the flexures, after
segmentation, often will be composed by both the ascending/descending and transverse part, hence
this will give large distances in the flexures.
33 of 50
CHAPTER 6. POST PROCESSING
Fig. 6.7: Since the components of the flexures can be composed of both ascending/descending and transverse
part, the majority of the δ points will be located in the flexures. In these cases the distance map will lack
reliability as an estimation of the diameter of the colon content.
6.4
Model-based Approach
A very simple model of the colon was made. It simply divided the colon content into three parts
corresponding to ascending, transverse and descending segments.
6.4.1
Approach
The approach is shown in Figure 6.8. By clicking in the two flexures of the colon three lines were
drawn, estimating approximately shape of the colon. A line connecting two selected points represented
the transverse part and two lines drawn vertically in the inferior direction of the image represented
the ascending and descending part of the colon. Angles formed in both of the colon corners were split
by axial lines, these two lines (division lines) intersected roughly in the center of the abdomen creating
a central point from which another division line was drawn vertically in the inferior direction. The
division lines now created three areas in the image plane, which could be used to cluster segmented
pixels into three parts of the colon.
Fig. 6.8: Estimating a rough distribution of the colon using lines created by choosing flexure points.
34 of 50
6.4. MODEL-BASED APPROACH
The division of colon content volume in three anatomical colon segments - ascending, transverse,
descending, could provide valuable information. Further work could include showing the areas with
largest accumulation of colon content. Knowing, that largest masses of content tend to cumulate in
the right side of the colon, especially in the connection of terminal ileum and cecum, could help in
evaluating results from segmentation.
However, due to heterogeneous manifestations of the locations and shapes of the colon content, this
model-based approach showed some difficulties.
6.4.2
Difficulties Regarding the Model Based Approach
A drawback to this approach was its sensitivity. According to the locations of the clicks of the mouse,
the total volume was divided into three partial volumes. To easily repeat a test of sensitivity, the two
points were always placed in the inner side of the colon flexure. This way the impact of placing flexure
points on the distribution of volume among the three main parts of the colon (ascending, transverse,
descending) was calculated. The test was repeated 10 times, with calculated results of percent of total
volume (TV) in each part of the colon in Table 6.1.
# measurement
1
2
3
4
5
6
7
8
9
10
µ
σ
% of TV in ascending
26.22
32.40
29.95
30.66
33.83
29.94
31.13
34.54
32.33
30.31
31.13
2.35
% of TV in transverse
54.27
50.41
51.05
48.80
47.85
49.13
49.27
46.78
49.35
50.77
49.77
2.05
% of TV in descending
19.51
17.19
19.00
20.54
18.31
20.94
19.60
18.68
18.32
18.92
19.10
1.10
Table 6.1: Percent of total volume in each part of colon. µ is the mean, and σ is the standard deviation.
Figure 6.9 shows a slice with colon content divided into three parts, differentiated by color. Mouseclicks in corner points are marked for each image. Since it is a 2D projection, it may be hard to
distinguish exactly the position of the flexure, while colon content tends to form larger masses in the
connection of the terminal ileum and cecum and also the flexures.
Fig. 6.9: Two dividings of colon according to different placements of corner points.
35 of 50
CHAPTER 6. POST PROCESSING
The image on the left results in dividing the colon into percentage of the total volume in ascending
33.18%, transverse 48.60% and descending 18.22%. The right image gives a result of 32.73% in
ascending, 45.21% in transverse and 22.06% in descending part of colon. The human colon is very
individual with various lengths and location of the single colon parts in the abdomen. In the example
shown in Figure 6.10 the transverse part of colon and its location would completely disable an accurate
use of this method, due to the placement of flexures above this segment. In this case, most volume of
the transverse part would be estimated as either ascending or descending.
Fig. 6.10: Problematic position of transverse segment for use of this method.
6.5
Downgrading the Complexity of the Measurements
The measurements requested by the physician seemed simple. However, in this case these rather
simple measurements required complex assumptions due to the nature of colon content. First of all
anatomically all individuals are different. Secondly, different eating habits and levels of exercise lead
to diverse manifestations of colon content. This complicates the statement of valid assumptions. The
following assumptions could ease the extraction of parameters for the physician, but as will follow,
the assumptions are invalid.
Colon Content is Present in All Segments
In some patients this assumption could be valid. However, even very constipated patients might have
the majority of colon content in the right half and only little in the left half, as seen in the results
Section 7.2. However, if such an assumption were to be made, more solutions to get the requested
measurements arise. Fitting to a simple model as suggested in section 6.4 could provide a valid
direction of the components leading to a possible calculation of the true cross-sectional area. On the
other hand many variations of the shape and location of colon exist, as shown in Figure 1.2.
Colon Content is Tubular
Colon content is tubular to some degree. The smooth muscle layer in the colon wall, Tunica muscularis, has an inner circular layer. The surface of this layer is the myentheric plexus, which contributes
to peristalsis, as described in Section 1.1. This feature of the colon wall could contribute to the colon
being predominantly circular. However, the segmented volumes can be far from circular in that the
concavities are dominant. This is shown in Figure 7.5, where the algorithm nicely segments the bright
36 of 50
6.6. ESTIMATION OF VOLUMES
colon content. Assuming that colon content is tubular would allow for valid estimates of the diameter
of the components based on a distance map or skeletonisation. The distance map could provide a
distribution of the diameters along a component. Skeletonisation, on the other hand, could provide
both a measure of lengths and a distribution of cross-sectional areas. The approach in section 6.2 was
not successful in the flexures, but using the skeletonisation might overcome that issue.
However, due to the diverse appearance of colon content, even simpler measures were aimed at. The
proposed algorithm calculates the total and two partial volumes. This description of volumes can
indicate a rough distribution, which might add clinical and robust value to the algorithm.
6.6
Estimation of Volumes
Previous attempts to add clinical value to the segmentation demonstrated limitations.
6.6.1
Approach
To calculate the total volume (TV) and some partial volumes, the following is needed.
1. Total voxels of each component
2. Volume of each component
The total voxels of each component is provided by the Matlab function regionprops and its method
Area.
When provided with the total number of voxels of each component the total volume of the component
is computed by the slice thickness, spacing between slices and the pixel spacing, which are read from
the dicomHeader.
In general the volume is calculated as the volume in the slices added to the volume between slices:
oneSliceVolume = volumeInSlice + volumeBetweenSlices
The volume in slice is computed from the slice thickness and the volume of the component in the slice,
i.e. the total amount of pixels multiplied with the pixel spacing.
To calculate the volume between slices the same assumptions are made as when calculating the area
between slices, thus the volume of the whole component is calculated as:
totalV olume =(numberOf P ixels(1) · pixelSpacing 2 ) · (sliceT hickness + spacingBetweenSlices · 0.5)
+ numberOf P ixels(n) · picelSpacing 2 ) · (sliceT hickness + spacingBetweenSlices · 0.5)
+
n−1
X
numberOf P ixels(i) · pixelSpacing 2 ) · (sliceT hickness + spacingBetweenSlices)
i=2
To find a distribution of the colon content the segmented content is divided into a left and a right
part by manually outlining the middle of the images. Afterwards a right and left volume is calculated
the same way as the total volume.
37 of 50
7
Results and Validation
To assess the performance of the proposed algorithm for segmentation of colon content, the accuracy
and reliability were validated. The accuracy was evaluated both quantitatively and qualitatively, while
the reliability was evaluated mainly quantitatively.
7.1
Results from Segmentation
Accurate classification of segmented tissues is important due to correct readings of the segmentation.
To assess the accuracy of a segmentation, knowledge about a ground truth, a voxel-wise gold standard,
can be useful. A physician was asked to manually outline regions of colon content on two MRI series
using the programme ImageJ (ROI Manager) as seen in Figure 7.1.
Fig. 7.1: Manual segmentation in ImageJ indicated by yellow.
These segmented regions were then imported to Matlab and compared to results from our automatic
segmentation using the Dice Similarity Coefficient and a qualitatively review of the segmentations as
described in the following section [30].
7.1.1
Test of Accuracy: Comparison to ground truth using Dice Similarity
Coefficients and a Qualitatively Validation
Dice similarity coefficients s were calculated based on the results of the algorithm and results from
two manual segmentations. With ground truth available, the s value is a measure of the accuracy
of a system, since it is a simple measure of spatial overlap between two sets. The value of s ranges
39
CHAPTER 7. RESULTS AND VALIDATION
from 0-1, with 0 indicating no spatial overlap and 1 indicating complete overlap between two sets of
segmentation results (X and Y ).
s=
2|X ∩ Y|
|X| + |Y|
(7.1)
s values were calculated based on the results of the algorithm and results from two manual segmentations (ground truth). Dice coefficient for the first series was s = 0.40 and s = 0.63 for the second.
Total volume segmented by manual outlining was 1 l in first series and 1.4 l in second series. Volume
estimation using the proposed algorithm gave results 0.8799 l and 0.6772 l respectively. When comparing s with estimated volumes, it is evident that in the first series (small s value, similar volume)
the algorithm segments a lot of false positives. In the second case (larger s value, greater difference in
volume) by visualising results from both segmentations, it is seen that the automatic has fewer false
positives, but also does not extract the dark areas as in the manual segmentation. Figures 7.2 and
7.3 show how the changes in s values slice by slice in both segmented series. In series 4 a rapid drop
occurs when manual segmentation outlined very dark areas and borders of colon content, see fig. 7.6.
Results from series 5 are more consistent, s ranging from 0.6–0.8 in the central slices, with drops in
the few first and last slices.
1
0.9
Dice similarity coefficient
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
10
20
30
Slice number
40
50
60
Fig. 7.2: s values in each slice, series 4
0.9
0.8
Dice similarity coefficient
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
5
10
15
20
25
30
Slice number
35
40
45
Fig. 7.3: s values in each slice, series 5
40 of 50
50
7.1. RESULTS FROM SEGMENTATION
In both series, ground truth contained also areas of dark intensities, where presence of colon content is
presumed by clinical experience of physician. Contours of the colon content were outlined, while it is
very hard and time demanding, requiring zooming in, to outline the gaps and holes of the content itself.
This is of course very challenging for any automatic algorithm. Colon content is presented as bright
intensities, with drops in the boundaries towards surrounding tissue. Including a wide span of pixel
intensities would at the same time cause a certain amount of leakage. In one series the region growing
algorithm segmented also areas of the terminal ileum, which had similar (and higher) intensities. In
a few slices the ileum (marked by yellow arrow) anatomically connects with the cecum, thus forming
a large connected component, as shown in Figure 7.4, recognised by the algorithm as colon content.
This creates a large false positive area, rapidly decreasing the s value.
Fig. 7.4: Bright intensity terminal ileum segmented as part of colon content
s-values can be calculated for each slice separately, allowing us to view specific parts of the image,
where results from automatic segmentation best agree with those from manual outlining. Figure 7.5
and 7.6 show the outlines of regions segmented manually (green) and using region growing algorithm
(red). Clearly can be seen, that red lines follow the bright edges of regions, while the green lines take
into account the whole bright region including darker intervening areas forming a rounder contour.
Fig. 7.5: Manual (green) and automatic (red) segmentation resulting in s = 0, 74
In the first Figure 7.5, the segmented region from automatic segmentation is similar to the manual,
in this case is s = 0, 74. The image nicely represents the difference between contours of colon content
and the visible parts represented by higher intensity voxels.
41 of 50
CHAPTER 7. RESULTS AND VALIDATION
Fig. 7.6: Manual (green) and automatic (red) segmentation resulting in s = 0, 39
In the second Figure 7.6 on the other hand a large difference between the algorithm result and result
from what the physician outlined as part of colon content from clinical experience.
Summed up, during this qualitative validation three major areas were found in which the segmented
areas differed from ground truth:
1. The proposed algorithm causes leakage in some series (false positives)
2. Ground truth includes a wider span of intensities.
3. Ground truth has a inconsistently tendency in some slices.
7.2
Results from Volume Estimation
Results from calculation of total volume of colon content, further divided into left and right according to
a manually selected center point of the patient are shown in Table 7.1. Volumes from both automatic
and manual segmentation are mentioned. As mentioned in previous section, manual segmentation
included whole contours of colon content with dark surrounding areas which the algorithm can not
segment. Therefore larger volumes in both manually segmented series (4* and 5*) were estimated
compared to results from automatic segmentation.
Data Set
#1
#2
#3
#4
# 4*
#5
# 5*
Total Volume [l]
0.8667
0.0775
0.6542
0.8799
1.0089
0.6772
1.4
Right Volume [l]
0.5091
0.0182
0.6223
0.5881
0.7935
0.4493
0.8835
Left Volume [l]
0.3575
0.0592
0.0319
0.2918
0.2154
0.2279
0.5156
Table 7.1: Estimation of volumes. *) Manual segmentation, ground truth.
Difference in ground truth to results from algorithm are also reflected in the calculated s values for
both series s = 0.40 and s = 0.63, respectively. Results also show, that the right side tends to contain
a larger volume of colon content. This could be due to lack of motility in the colon of constipated
patients and the fact, that content tends to accumulate in the cecum (where content enters the colon
42 of 50
7.2. RESULTS FROM VOLUME ESTIMATION
after leaving the terminal ileum) and flexures of the colon. An illustration of this situation is shown
in Figure 7.7. Long term constipation on the other hand involves the entire length of the colon.
Fig. 7.7: Total volume segmented from data set # 5
7.2.1
Test of Reliability: Comparing Consecutive Scans of one Subject
Reliability refers to the reproducibility of the results from a test [31], which has been reported as an
important measure of confidence in segmentations [13]. Better reliability implies better detection of
changes in measurements. The standard deviation is a simple measure of variation from scan to scan.
To obtain a trustworthy estimate of the reliability many scans of many subjects are required. However,
due to practicalities, only 3 scans of 2 subjects were available. Unfortunately only the scans from one
subject were useful, since the scans from the other subject barely had any visible colon content.
Due to a different magnetic field strength, the algorithm was adjusted concerning two points. Firstly,
the global thresholds for identifying candidate initial seed points were manually selected for each
series. The final initial seed points were chosen as described in Section 2.2. Secondly, the α threshold
parameter was identical with that from the 1.5 T scanner, while the β threshold parameter was tuned
to the 3 T scanner. Hence, the following parameters were used: α = 4 and β = 3
Results from the 3 volume estimations are shown in Table 7.2 with mean (µ) and standard deviation
(σ) for the total volume and left and right volumes.
Total Volume [l]
Right Volume [l]
Left Volume [l]
1
0.5625
0.5460
0.0164
2
0.5400
0.5282
0.0118
3
0.5157
0.4893
0.0264
µ
0.5394
0.5212
0.0132
σ
0.0234
0.0290
0.0075
Table 7.2: Test of reliability.
Visible consistency in the resulting volumes can conclude and indicate a certain level of reliability
while segmenting consequent scans from the same subject. The volume ratios are preserved for both
left and right volume. Figure 7.8 shows regions where region growing leaked into surrounding tissue,
forming a false positive result, which added up to the total volume and subsequently either left or
right volume.
43 of 50
CHAPTER 7. RESULTS AND VALIDATION
Fig. 7.8: A leak in the right side of the patient was seen in all 3 series used to test the reliability.
7.2.2
Computational Times
What is most important, accuracy or speed? In surgical planning, accuracy might come first. However,
regarding assessment of colon content or treatment follow-ups both speed and accuracy are important,
and the crucial factor could be the level of user intervention. The physician might want a accurate
result, but he or she does not have time for higher levels of intervention. Balancing accuracy with the
level of intervention is really the task. The proposed algorithm is semi-automatic, and all the physician
has to do is to review if only true positives are present. In the following computational times for the
algorithm are presented. A comparison with the time taken to do a manual segmentation is made.
The total computational time for each series and the average computational time are shown in Table
7.3.
Time [min:s]
1
35:24
2
24:17
3
46:05
4
36:03
5
30:23
µ
34:26
Table 7.3: Total computational times for each series.
The two manual segmentations took 41:20 and 54:27 min, respectively, with a average of µman = 47:54
min.
In the automatic segmentation attention from the physician is only needed while drawing the regions
of interest for selection of initial seed points.
44 of 50
8
Discussion
This paper proposes a novel method for semi-automatic segmentation and quantification of colon
content from MRI based on region growing. The quantification of colon content relies on a well
executed segmentation of the colon content. The algorithm has proved successful in segmenting bright
colon content allowing for a robust quantified representation of the colon content.
8.1
Segmentation
Semi-automatic quantification of colon content from MRI is proposed using a region growing algorithm. The challenge is in segmenting disconnected components with irregular shapes having no or
limited knowledge about their position in the abdomen. Series also differ in intensity representation of
colon content, making it necessary to adjust parameters of our segmentation algorithm. The proposed
algorithm for segmentation of colon content from MRI, is dependent on a minor user intervention.
The operator has to draw a region of interest for planting of seeds. The level of user intervention
compromises the accuracy of the segmentation. In some scans artifacts and unwanted tissues have
similar or even higher intensities than that of the colon content. When drawing regions of interest the
unwanted areas are excluded resulting in an increased number of false positives in the result.
The proposed algorithm utilises region growing to expand the initially segmented areas used as seed
points. The growing of these regions applies intensity differences and gradient differences between
colon content and surrounding tissue including air in the colon lumen. Hence, the segmentation
requires colon content to be distinguishable from neighboring tissue. This holds true for both the
placement of seed points which is based on a fuzzy c-means clustering and the first of two growing
criteria in the region growing. The results show that the algorithm is successful in segmenting the
desired colon content. However, in MRI series with a low contrast, the algorithm might leak. This
will result in an incorrect classification of tissues, and hence an over-segmentation. To hinder this a
second growing criterion is formed. This relies on the gradient magnitudes of the images. The colon
content tends to consist of areas of several high gradient peak areas. This property is used to control
the region growing by restricting the growing into these areas. Still, some over-segmentation can occur
due to closely located tissues that are similar to colon content in intensities.
The algorithm relies on universal segmentation criteria. This however results in some parts of the
segmentation leaking to surrounding tissue areas. The present results demonstrate that region growing
is capable of segmenting the bright colon content. Differences occur between automatic segmentation
and the ground truth obtained by manual segmentation performed by a physician. Differences are due
to the core of the region growing algorithm being able to segment bright colon content from darker
surrounding matter. Manual segmentation on the other hand takes into account prior anatomical
knowledge of the human colon, therefore outlining contours of the colon content, disregarding the
concavities at the edges of the bright matter. In general, manual outlining segmented a larger volume
of colon content. Automatic segmentation includes false positive bright tissue (little intestine) as colon
content and demonstrates areas of leakage towards the abdomen wall and other soft tissues.
45
CHAPTER 8. DISCUSSION
8.2
Automatic Versus Manual Segmentation
The developed algorithm for segmentation of colon content from MRI, is dependent on a minor user
intervention. The operator has to draw a region of interest, in which the algorithm plants seeds for
region growing. This was chosen in order to increase the amount of false positives in the final segmentation. Without this implementation parts of the liver and artifacts were included. Wyatt et al.
have developed a method for automatic segmentation of the colon for virtual colonoscopy [32]. They
took advantage of the high contrast in the gas/tissue interface in helical computed tomography (CT)
images. The gas-filled lumen appears dark compared to the neighboring tissue. Prior knowledge of
colon geometry was used to localise seed points. They achieved a complete but not highly accurate
representation of the colon.
Backmana et al. suggest that segmentation algorithm for use in medical domains should only have a
minor user intervention [33]. The operators which skills are specialised into the medical field should not
be burdened with modifying parameters, selecting thresholds, or other tasks related to prior knowledge
of the algorithm and application field. They also state that current supervised segmentation algorithm
are superior to the unsupervised algorithms in accuracy, which indicates that the degree of operator
intervention and accuracy of segmentation should be compromised. The algorithm in the work of
Backmana et al. requires one operator intervention, the selection of a single seed voxel. Mancas et al.
developed a semi-automatic region-growing segmentation method in which the operator also has to
select a seed inside the region of interest [34]. Their method is based on both a topology distance map
and a probability map, and when the structure in the image is homogeneous according to a criterion,
they are able to segment the spinal cord.
From this review, the minor operator intervention required in present project is considered as acceptable.
8.3
Sensitivity to Selection of Seed Point
As for now the algorithm is sensitive to the physician picking the right seed points for the region
growing. When drawing a region of interest around a component, the whole 3D component is selected
as seed points, that is choosing a component that are not only true positive will lead to an oversegmentation. A way to overcome this problem could be to draw region of interest around a component,
that both are colon content and not, when it appears on a slice as only true positive, letting this 2D
component be selected as seed points. This will make it possible to choose more seed points in colon
content. The second threshold criterion will hinder leakage into the tissue that was connected with
the outlined 2D component from FCM. However, utilising a manual intervention will always indicate
a certain degree of sensitivity and subjectivity.
8.4
Optimising the Algorithm for Clinical Use
Even though the test of computational times revealed quite good results compared to a manual segmentation. Not only do the times taken to do the segmentations differ, hence µaut = 34:26 min and
µman = 47:54 min., the level of attention from the user also differs widely. A few initiatives might
decrease the computational time even more. Matlab does not fully compile its for-loop, which results
in for-loops tending to be slow. By vectorising the code, the computational time might decrease. Another issue heavily affecting the computational time is the FCM clustering. To locate 5 cluster centers
demands more than 100 iterations. By integrating prior knowledge the number of needed iterations
might decrease resulting in a faster algorithm.
The perspective of the proposed algorithm is to provide an objective measurement for quantification of
colon content. This could support the clinical assessment of constipation. The method and its results
are very promising as new objective measurements for physicians. Future work targets an improved
segmentation, possibly by including additional image features, fine tuning of parameters, and reducing
46 of 50
8.4. OPTIMISING THE ALGORITHM FOR CLINICAL USE
the sensitivity of the manual intervention. Also, further investigations of limitations of the proposed
method and further comparison with current clinical practice should be conducted.
47 of 50
Bibliography
[1]
Committee for Medicinal Products for Human Use. Concept paper on the need of a guideline for
clinical investigation of medicinal products for the treatment of chronic constipation. European
Medicines Agency. 2012;.
[2]
World Gastroenterology Organisation. World Gastroenterology Organisation Practice Guidelines:
Constipation. 2007;.
[3]
Lunn S. The Colon; 2012. Available from: http://www.sharonlunn.co.uk/.
[4]
Cummings JH, Hill MJ, Jenkins DJA, Pearson JR, Wiggins HS. Changes in fecal compostion and
colonic function due to cereal fiber. The American Journal of Clinical Nutrition. 1976;29:1468–
1473.
[5]
Cummings JH, Macfarlane GT. The control and consequences of bacterial fermentation in the
human colon. Journal of Applied Bacteriology. 1991;70:443–459.
[6]
Geneser F. Histologi. 1st ed.; 2008.
[7]
Harris T. The X-Ray Machine. 2012;Available from: http://science.howstuffworks.com/
x-ray2.htm.
[8]
Alessandro DMD. What Causes Recurrent Abdominal Pain?; 2012. Available from: http://www.
pediatriceducation.org/2004/12/27/.
[9]
Stern DP. Educational Web Sites on Astronomy, Physics, Spaceflight and the Earth’s Magnetism;
2012. Available from: http://www.phy6.org.
[10] Rinck P. Magnetic Resonance in Medicine. The Basic Textbook of the European Magnetic
Resonance Forum.; 2012. Available from: http://www.magnetic-resonance.org/MagRes%
20Chapters/03_10.htm.
[11] MedicalExpo. The Virtual Medical Exhibition; 2012. Available from: http://www.medicalexpo.
com/product-manufacturer/toshiba-mri-coil-759-22.html.
[12] MRI Sequences;Available from: http://www.bioc.aecom.yu.edu/labs/girvlab/nmr/course/
COURSE_2012/chapter5.pdf.
[13] Clarke LP, Velthuizen RP, Camacho MA, Heine JJ, Vaidyanathan M. MRI Segmentation: Methods and Applications. Magnetic Resonance Imaging. 1995;13.
[14] Gray level co-occurrence matrix;Available from: http://www.mathworks.se/help/images/ref/
graycomatrix.html.
[15] Jan J. Medical Image Processing, Reconstruction and Restoration; 2005.
[16] Connell AM, Hilton C, Irvine G, Lennard-Jones JE, Misiewicz JJ. Variation of Bowel Habit in
Two Population Samples. British Journal of Medicine and Medical Research. 1965;.
49
BIBLIOGRAPHY
[17] Johanson JF.
Definitions and Epidemiology of Constipation. 2007;Available from:
http://www.researchgate.net/publication/226057621_Definitions_and_Epidemiology_
of_Constipation.
[18] Rostgaard J, Tranum-Jensen J, Qvortrup K, Holm-Nielsen P. Hovedets, halsens og de indre
organers anatomi. 10th ed.; 2006.
[19] Martini FH, Nath JL. Fundamentals of anatomy & physiology; 2009.
[20] Panchal SJ, Müller-Schwefe P, Wurzelmann JI. Opioid-induced bowel dysfunction: prevalence,
pathophysiology and burden. International Journal of Clinical Practice. 2007;61, 7:1181–1187.
[21] Guyton A, Hall J. Textbook of Medical Physiology. 11th ed.; 2006.
[22] Thompson AC, Attwood DT, Gullikson EM, Howells MR, Kortright JB, Robinson AL, et al.
X-Ray Data Booklet. 2009;Available from: http://xdb.lbl.gov/.
[23] Bushberg JT, Seibert JA, Edwin M Leidholdt j, Boone JM. The Essential Physics of Medical
Imaging. 2nd ed.; 2002.
[24] Carlton RR, Adler AM. Principles of Radiographic Imaging. 4th ed.; 2006.
[25] Reuchlin-Vroklage M, Bierma-Zeinstra S, Benninga MA, Berger MY. Diagnostic Value of Abdominal Radiography in Constipated Children. Archives of Pediatrics & Adolescent Medicine.
2005;159.
[26] Moylan S, Armstrong J, Diaz-Saldano D, Saker M, Yerkes EB, Lindgren BW. Are Abdominal
X-Rays a Reliable Way to Assess for Constipation? The Journal Of Urology. 2010;184.
[27] Bongers MEJ, Voskuijl WP, van Rijn RR, Benninga MA. The value of the abdominal radiograph
in children with functional gastrointestinal disorders. European Journal of Radiology. 2006;59.
[28] Health Physics Society. Radiation Exposure From Medical Diagnostic Imaging Procedures. Health
Physics Society Fact Sheet. 2010;.
[29] Andrzej Materka MS. Texture Analysis Methods - A Review. 1998;.
[30] Dice LR. Measures of the Amount of Ecologic Association Between Species. Ecology. 1945;26:297–
302.
[31] Hopkins WG. Measures of Reliability in Sports Medicine and Science. The American Journal of
Sports Medicine. 2000;30.
[32] Wyatt CL, Ge Y, Vining DJ. Automatic segmentation of the colon for virtual colonoscopy.
Computerized Medical Imaging and Graphics. 2000;24.
[33] Backmana NJ, Whitney BW, Furstc JD, Raicuc DS. A prioritized and adaptive approach to
volumetric seeded region growing using texture descriptors. Electronic Imaging. 2006;6065.
[34] Mancas M, Gosselin B, Macq B. Segmentation using a region-growing thresholding. In: Image
Processing: Algorithms and Systems; 2006. p. 388–398.
50 of 50