INVESTIGATION OF SURFACE EVOLUTION FOR
ABDOMINAL AORTIC ANEURYSMS INTERACTING WITH SURROUNDING TISSUES
By
Alexander C. Dupay

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
MASTER OF SCIENCE
Engineering Mechanics
2012

ABSTRACT
INVESTIGATION OF SURFACE EVOLUTION FOR
ABDOMINAL AORTIC ANEURYSMS INTERACTING WITH SURROUNDING TISSUES
By
Alexander C. Dupay
Three-dimensional patient-specific models of abdominal aortic aneurysms (AAAs) were
reconstructed from CT scans and vessel migration, eccentricity, and cross-section line curvature
were identified as morphological features which are influenced by surrounding tissues during
disease progression. A novel centerline algorithm was developed in Matlab and applied capable
of qualitatively estimating vessel migration away from the vertebrae during disease progression
and automatically parameterizing patient-specific surface models with respect to the generated
centerline.

The parameterized surface is then approximated using shape functions and

reconstructed appropriately.

Image analysis and reconstruction was performed through a

commercial biomedical imaging software suite for a high-resolution CT dataset consisting of 51
scans across 11 patients which were obtained through collaboration with Dr. Whal Lee at Seoul
National University Hospital (SNUH). A longitudinal dataset, multiple scans taken at multiple
times per patient, presents additional challenges such as the need for accurate model registration
and the unique opportunity of increasing our understanding of geometrical evolution during the
progression of AAA. This study discovers that a rich and complex biomechanical relationship
exists between an AAA and the spine during AAA disease progression. The most significant
implication of these results highlights the necessity of incorporating surrounding tissue, whether
patient-specific or idealized, for computational biomechanics stress simulations of AAAs to
improve their accuracy and prediction of rupture potential. Incorporating such interactions
would also considerably enrich documented analyses.

ACKNOWLEDGEMENTS

I would like to thank Dr. Seungik Baek for his mentorship and the vast amount of knowledge and
support he has imparted on me throughout our work together.

I would like to thank my

committee members Dr. Tamara Reid Bush and Dr. Thomas Pence for their time and insight
with regards to this work. The same can be said for the colleagues I’ve had the pleasure of
working with in my lab: Dr. Shahrokh Zeinali-Davarani, Dr. Jungsil Kim, Byron Zambrano R,
Paul Snyder, and Jared Kavinsky. They all helped me further understand our field of study and
helped shape my own studies as we worked together. Upon shaping my studies our collaboration
with Dr. Lee was the invaluable basis of this thesis, and for that I am grateful. I had the pleasure
of working with Drew Kim and Angela Kolonich during the Summer semester of my studies
through the Research Education for Teachers (RET) program, and would like to acknowledge the
support given to me by Drew Kim and the research help rendered by Angie. Outside the lab I
would like to thank my friends and family and foremost my wife Cori, who have been there to
support me every step of the way.

iii

TABLE OF CONTENTS

LIST OF TABLES ......................................................................................................................... vi
CHAPTER 1. Introduction.............................................................................................................. 1
1.1 Motivation ........................................................................................................................ 1
1.1.1 Clinical Review of AAAs .............................................................................................. 1
1.1.2 Biomechanics Research of AAAs .................................................................................. 5
1.1.3 Surrounding Tissues..................................................................................................... 11
1.2 Objectives ........................................................................................................................... 14
1.3 Contributions....................................................................................................................... 16
CHAPTER 2. Biomedical Image Processing: Segmentation and Registration ............................ 19
2.1 Introduction ......................................................................................................................... 19
2.1.1 Biomedical Image Processing Background ................................................................. 19
2.1.2 Biomedical Imaging Modalities ................................................................................... 20
2.2 Methods............................................................................................................................... 25
2.2.1 Data .............................................................................................................................. 25
2.2.2 Segmentation................................................................................................................ 30
2.2.3 3-D Modeling ............................................................................................................... 34
2.3 Results ................................................................................................................................. 38
2.4 Discussion ........................................................................................................................... 47
CHAPTER 3. Morphological Features Influenced by Surrounding Tissue ................................. 50
3.1 Introduction ......................................................................................................................... 50
3.2 Method ................................................................................................................................ 52
3.3 Results ................................................................................................................................. 54
3.4 Discussion ........................................................................................................................... 61
CHAPTER 4. Parameterization of Luminal and AAA Surfaces .................................................. 64
4.1 Introduction ......................................................................................................................... 64
4.2 Methods............................................................................................................................... 65
4.2.1 Importing the data ........................................................................................................ 65
4.2.2 Centerline Generation Algorithm ................................................................................ 67
4.2.3 Smooth Centerline Approximation .............................................................................. 72
4.2.4 Longitudinal Acquisition Plane Generation................................................................. 76
4.2.5 Parameterization of Cross Sections into r(s, θ) ............................................................ 80
4.2.6 Coordinate transformation from r(s, θ) to (x, y, z) ...................................................... 82
4.2.7 Surface Approximation ................................................................................................ 84
4.2.8 Eccentricity and Polyline Curvature ............................................................................ 88
4.3 Results ................................................................................................................................. 90
4.3.1 Centerline generation ................................................................................................... 90
4.3.2 Eccentricity and Curvature .......................................................................................... 92
iv

4.4 Discussion ........................................................................................................................... 94
4.4.1 Centerline Generation Algorithm ................................................................................ 94
4.4.2 Parameterization and Surface Approximation ............................................................. 97
4.4.3 Morphological Parameters ........................................................................................... 99
4.4.4 Future Work ............................................................................................................... 102
APPENDICES ............................................................................................................................ 105
APPENDIX A. Detailed information and measurements of CT scan data............................. 105
APPENDIX B. MatLab code developed for automatic parameterization of AAAs............... 111
APPENDIX C. Results of centerline and surface approximation algorithm performance ..... 148
REFERENCES ........................................................................................................................... 159

v

LIST OF TABLES

Table 1. Patient age, number of scans, and sex after discounting irrelevant scans. ..................... 27
Table 2. Statistical summary of the CT dataset discounting irrelevant scans. The surveillance
interval refers to the time between individual scans, such as between J1 and J2 or J4 and J5. The
total study length refers to the total time between the first and last scan of a patient, such as
between J1 and J5. ........................................................................................................................ 27
Table 3. Optimal STL model generation parameters used for every STL model generated in the
study. Contour interpolation results in a smoother model with less gaps and is recommended for
medical CT applications [59]. Matrix Reduction is set to 1, equivalent to no matrix reduction,
which inherently nullifies the need to select a matrix reduction algorithm as well (not shown). 36
Table 4. The number of employed smoothing iterations at a constant smoothing factor of 0.7 on
each STL model after optimal model generation. The smoothing was limited to a maximum of
10 iterations to avoid significant loss of detail. ............................................................................ 36
Table 5. Summary of the definition of techniques employed over the survey of 56 studies by
approximate percentage [51–53]. Some groups reported multiple methods of measurement to be
in use rather than a single method, resulting in total percentages not equal to 100%. ................. 52
Table 6. Linear approximation of annual transverse maximum AAA tissue diameter growth and
R2 fit value. ................................................................................................................................... 56
Table 7. Linear approximation of annual transverse maximum AAA lumen diameter growth and
R2 fit value. ................................................................................................................................... 59
Table 8. Conditional table for shape function assignment of smooth line approximation. .......... 74
Table 9. Conditional table for calculation of θ from a given rh vector. First it is determined
whether θ lies on the N2, N3, -N2, or -N3 axis. If it does not lie on an axis than the region it lies in
is determined and θ is calculated appropriately. ........................................................................... 81
Table A.1. Detailed scan resolution information for the CT data set studied. ............................ 106
Table A.2. Detailed scanner manufacturers and models, date of scans, and total time elapsed
from initial scanning. .................................................................................................................. 107
Table A.3. Tissue Measurements from individual patient CT scans. ......................................... 109
Table A.4. Lumen Measurements from individual patient CT scans. ........................................ 110

vi

Table C.1. SSE values from surface approximations sorted first by Nη and then by Nθ ............. 152
Table C.2. SSE approximations sorted based on SSE values ..................................................... 155

vii

LIST OF FIGURES

Figure 1. A sample CT scan of a patient diagnosed with an AAA highlighting a coronal slice
(left, upper-left), a sagittal slice (left, upper-right), a transverse slice (left, bottom), and the
intersection of all three (right). White arrows are used to identify the AAA. For interpretation of
the references to color in this and all other figures, the reader is referred to the electronic version
of this thesis The text is not meant to be readable this figure and other CT reference figures, and
is for visual reference only. ............................................................................................................. 2
Figure 2. Transverse slice of a patient CT scan in the abdominal region (left), and a zoomed
section of a transverse slice in the abdominal region (right) focusing on AAA transverse crosssections exhibiting varying thickness of the arterial wall due to ILT presence shown by the
arrows. The intensity scale is adjusted on the right to show clear contrast between lumen, tissue,
and intraabdominal space. ............................................................................................................... 5
Figure 3. Near the lumbar vertebrae (white) the inferior vena cava (blue) runs opposite to the
aorta (red) on its right lateral side. The inferior vena cava and the aorta are connected, most
notably for AAAs, to the renal system (purple). The left renal vein crosses anteriorly over the
aorta inferior to the superior mesenteric artery. ............................................................................ 12
Figure 4. A process diagram highlighting key elements of the biomedical imaging process from
data obtainment to generating technical computing models. ........................................................ 20
Figure 5. A 3-D left side view of F1 after initial processing. The lower scan boundary is
equivalent to the bottom of the image which clearly does not contain the abdominal cavity nor
the full presence of an AAA. The scan F1 was therefore discounted from this study on the basis
of insufficient data. ....................................................................................................................... 28
Figure 6. A 3-D transparent frontal view of the tissue surface (green, left) and the ELG (blue,
left), and the luminal surface (red, right) and the ELG (blue, right) of scan H8. The presence of
an ELG indicates this to be data from a post-EVAR scan, which is not the focus of this study.
The scans H8 – H9 were therefore discounted from this study on the basis of irrelevant data. ... 29
Figure 7. a 3-D view overlaid on a respective CT scan coronal plane for scan K4 (left) and for
scan K6 (right). The shape of the AAA from scan K6 does not represent natural AAA disease
progression; it instead more closely represents a post-EVAR geometry. The focus of this study
is pre-EVAR AAAs and K6 is therefore fully discounted. ........................................................... 29
Figure 8. A 3-D view overlaid on a respective CT scan coronal plane for scan E1. A thoracic
aneurysm is present, but no AAA is present. The focus of this study is AAAs, so the scans are
discounted. .................................................................................................................................... 30

viii

Figure 9. A segmented mask for the arterial lumen (red) and arterial tissue (green) of scan B2
shown for a coronal slice (left) and transverse slice (right). The transverse slice (right) highlights
the rectangular nature of voxels which occupy a segmented mask. The arterial tissue is seen to
contain a wide region on the upper-left, or left anterior, side indicating the presence of an ILT. 33
Figure 10. A 3-D AAA lumen STL reconstructed from scan B2. STL format describes a
triangularly meshed surface, as can be clearly seen. The data exported from the 3-D models for
Chapters 3 & 4 uses point cloud format. This format will only contain the vertex information
from STL models. ......................................................................................................................... 35
Figure 11. Front coronal plane view of 3-D luminal STL models generated for all 11 patients.
From top-left to bottom-right in reading order the STLs originate from scans: A2, B3, C2, D3,
E2, F7, G6, H7, I6, J6, K5. ........................................................................................................... 39
Figure 12. Right sagittal plane view of 3-D luminal STL models generated for all 11 patients.
From top-left to bottom-right in reading order the STLs originate from scans: A2, B3, C2, D3,
E2, F7, G6, H7, I6, J6, K5. ........................................................................................................... 40
Figure 13. A 3-D STL model of a reconstructed infrarenal AAA luminal surface as seen from the
left sagittal plane (left), front coronal plane (center), and right sagittal plane (right). A bounding
box indicates the infrarenal region specifically to be exported to a technical computing software
package. ........................................................................................................................................ 41
Figure 14. An example of two registered vertebral models from scans I1 & I3 as viewed from the
left (left), front (center), and right (right). Each model is the same color, and has been registered
to match the position of its pair. Near the top and bottom of the models clipping can sometimes
be observed, which highlights the boundaries of a model (red). In general, however, it is clear
that the models’ positions match extremely closely given that they are independently generated
from different CT scans. The associated lumen and tissue models are also automatically
registered (not shown) according to the transformation matrix generated by minimizing the
distance field between vertebral models. ...................................................................................... 42
Figure 15. Transparent and registered front coronal plane view (left) and right sagittal plane view
(second from left) of 3-D luminal STL models for B1 and B3. Transparent and registered front
coronal plane view (second from right) and right sagittal plane view (right) of 3-D tissue STL
models for B1 and B3. The lumen and tissue models presented are not to scale, as the lumen
models are inherently always smaller than the tissue models....................................................... 43
Figure 16. Transparent and registered front coronal plane view (left) and right sagittal plane view
(right) of 3-D luminal and tissue STL models for B3. .................................................................. 44
Figure 17. A transverse slice un-zoomed (upper-left) and zoomed (upper-right) from scan A2. A
3-D reconstruction of the AAA tissue geometry for scan A2 from the right (bottom-left) and an
opened version of the AAA tissue geometry for scan A2 from the right (bottom-right). Red
circles indicate linked areas between images. Within the red circles it is clear that the AAA
tissue is interacting with the spine and that the local geometry of the tissue is affected. ............. 45
ix

Figure 18. A transverse slice un-zoomed (upper-left) and zoomed (upper-right) from scan A2. A
3-D reconstruction of the AAA tissue geometry for scan A2 from the left (bottom-left) and an
opened version of the AAA tissue geometry for scan A2 from the left (bottom-right). Red circles
indicate linked areas between images. Within the red circles it is clear that the AAA tissue is
interacting with the spine and that the local geometry of the tissue is affected............................ 46
Figure 19. Reconstructed 3-D models for the vertebrae (white), AAA lumen (red), and inferior
vena cava tissue/lumen (blue) as seen from the left (left), front (center), and right (right) for scan
B2. ................................................................................................................................................. 47
Figure 20. A sagittal slice (left) of scan J5 and a transverse slice (right) of scan J5 with the
segmented lumen/tissue (red/green). The blue line indicates an anteroposterior maximum AAA
tissue diameter measurement axis in the transverse plane (left, right), the solid red line indicates
an anteroposterior maximum AAA tissue diameter measurement in the longitudinal plane (left),
the dashed red line indicates a crude longitudinal axis approximation (left), and the orange line is
simply the AAA tissue diameter measurement at 90° in the transverse plane to the
anteroposterior measurement axis. The caliper positioning for this case was external/tissue. .... 53
Figure 21. A transverse slice of scan D2 with the segmented lumen (red), maximum
internal/lumen diameter axis (blue) and axis 90° to the maximum internal/lumen diameter (red).
It is clear that anteroposterior measurements will not yield maximal results in the presence of an
eccentric ILT. This is the reason that the maximum diameter axis was used and anteroposterior
was not strictly used in all cases. .................................................................................................. 54
Figure 22. 3D CT reconstruction of patient C as viewed from the right. Two aneurysmal sacs are
present and the primary is referred to as Patient C, while the secondary is referred to as Patient C
(sac). .............................................................................................................................................. 55
Figure 23. The maximum transverse diameter of the arterial tissue vs. scanning date for each
patient. ........................................................................................................................................... 56
Figure 24. The annual expansion rate of maximum AAA tissue diameter growth calculated on a
successive scan basis for each patient. .......................................................................................... 57
Figure 25. The tissue eccentricity based on eccentricity = (Dmax) / (D @ 90° to Dmax) for each
patient scan.................................................................................................................................... 58
Figure 26. The maximum transverse diameter of the arterial lumen vs. scanning date for each
patient. ........................................................................................................................................... 59
Figure 27. The annual expansion rate of maximum AAA lumen diameter growth calculated on a
successive scan basis for each patient. .......................................................................................... 60
Figure 28. The lumen eccentricity based on eccentricity = (Dmax) / (D @ 90° to Dmax) for each
patient scan.................................................................................................................................... 61
x

Figure 29. The coordinate planes inherited from the anatomical planes where the coordinate axis
is such that the x-y plane is parallel to the transverse plane (green), the y-z is parallel to the
sagittal plane (red), and the x-z plane is parallel to coronal plane (blue). .................................... 66
Figure 30. The arterial surface (blue) with two maximally inscribed spheres (red) and their
respective constraint planes (green) shown. The centerline (black) is shown within the lumen
model for scan B2 as seen from an isometric point of view (left), parallel to the x-z or coronal
plane (center), and the y-z or sagittal plane (right). ...................................................................... 67
Figure 31. The lumen point cloud model for scan B2 colored by z-coordinate values with a red
sphere at approximately (155, -180, -163) (left), and the sphere itself colored by z-coordinate
values with vectors indicating the orientations of vectors A (blue), B (black), and C (red) for a
sample center point (right). ........................................................................................................... 70
Figure 32. The final collection of centerline points plotted within the lumen point cloud model
for scan B2 colored by z-coordinate value as seen from an isometric point of view (left), parallel
to the x-z or coronal plane (center), and the y-z or sagittal plane (right). .................................... 71
Figure 33. The meaning of p, d, and l can be hard to ascertain from equations alone. The
variable d represents the point to point distances between p values. The variable l represents the
sum of d values thus far at a given l.............................................................................................. 72
Figure 34. The final collection of centerline points plotted within the point-cloud lumen model of
scan B2 colored by z-coordinate value and the fitted centerline (black) as seen from an isometric
point of view (left), parallel to the x-z or coronal plane (center), and the y-z or sagittal plane
(right). ........................................................................................................................................... 76
Figure 35. The vectors N1 (red), N2 (blue) and N3 (green) are plotted for each s originating from
the points Xk for k = [1, m] for the lumen model of scan B2. ....................................................... 77
Figure 36. The imported point-cloud lumen model composed of points qN for scan B2 colored by
z-coordinate value with points on the smooth approximated centerline shown (black, left), and
the subset of points vh colored by z-coordinate value from the point-cloud lumen model for scan
B2 which lie on longitudinal planes defined by (N1, N2, N3) (right). For example, the vectors in
rh from Xk to the associated points on the surface vh at the plane defined by (N1, N2, N3) for k = 40
are shown in green (right). ............................................................................................................ 78
Figure 37. The radii values are calculated for all h by ||rh||. The mean is then calculated at each
longitudinal plane for the lumen of scan B2 (left), and the circles of mean longitudinal radii at
each longitudinal plane can be visualized for the lumen of scan B2 (right). ................................ 80
Figure 38. Schematic representation for the calculation of θ from a given rh vector. First it is
determined whether θ lies on the N2, N3, -N2, or –N3 axis (left), and if it does not lie on an axis
then the region it lies in is determined and θ is calculated appropriately (right). ......................... 82

xi

Figure 39. Parameterized surface data (red) from scan B2 shown as r(s, θ) surface. ................... 83
Figure 40. Subset of the original point cloud data which lie on longitudinal planes (left) and the
data reconstructed from the r(s, θ) data (right). ............................................................................ 84
Figure 41. A sample arterial lumen surface from scan C2 (left) is shown to undergo the steps of
parameterization (center) followed by surface approximation (right) with the optional
reconstruction to Cartesian coordinates shown for clarity in both cases. ..................................... 85
Figure 42. Parameterized surface data (red) and the respective points from the best fit surface
(blue) from scan B2 shown as r(s, θ) surfaces. ............................................................................. 87
Figure 43. Plot of the residual error between the parameterized surface data and the respective
points from the best fit surface from scan B2 shown as an r(s, θ) surface. .................................. 87
Figure 44. Parameterized surface data (red) and the respective cross-sectional lines from the best
fit surface (blue) from scan B2 shown as r(s, θ) surfaces. ............................................................ 88
Figure 45. Major axis (left, red) and minor axis (left, green) of a cross-section (blue) used for
eccentricity determination. Terminology of the 3-point polyline curvature method is detailed
(right). ........................................................................................................................................... 89
Figure 46. Isometric view of the centerlines (left) and a transverse view (right) for scans H1 – H5
where black, blue, green, red, cyan, magenta represent H1, H2, H3, H4, H5, H6 respectively. In
this case the general region of the spine is the x-z plane at y = 100 such that the arrow indicates
migration away from the spine. .................................................................................................... 91
Figure 47. x-z/frontal/coronal view of the centerlines (left), and y-z/sagittal view of the
centerlines (right) for scans H1 – H5 where black, blue, green, red, cyan, magenta represent H1,
H2, H3, H4, H5, H6 respectively. In the x-z plane (left) migration of the vessel centerline is
taking place towards the left, and in the y-z plane (right) migration is also taking place towards
the left. In this case the general region of the spine is the x-z plane at y = 100 such that the arrow
indicates migration away from the spine. ..................................................................................... 92
Figure 48. Eccentricity parameter values displayed for B2 lumen before surface approximation
with a cutoff of 1.4 (left) and after surface approximation (right). ............................................... 93
Figure 49. Curvature parameter values displayed for B2 lumen after surface approximation with
a cutoff of mean + 1 standard deviation. ...................................................................................... 94
Figure 50. SSE values curve sorted by Nη and Nθ for 100 cases of shape function surface fitting.
....................................................................................................................................................... 99
Figure 51. SSE values curved sorted by SSE value for 100 case of shape function surface fitting.
....................................................................................................................................................... 99

xii

Figure 52. Eccentricity of the arterial lumen surface of B2 pre-fit (blue) and post-fit (red). ..... 101
Figure C.1. A 3-D model of the lumen for D3 as seen in the x-z plane (upper-left), a coronal, or
x-z, slice from scan D3 (upper-right), an x-z (bottom-left), and a y-z (bottom-right) view of the
centerline algorithm results for scan D3. It is clear that scan D3 features a geometric anomaly.
This anomaly is highly unusual and tends to force the centerline generation to follow the
observed path regardless of input parameters. ............................................................................ 149
Figure C.2. An x-z (left) and y-z (right) view of the centerline algorithm results from scan F5. If
the final approximation of a centerline point differs greatly from the initial guess, the vector
between the previous centerline point and current centerline point becomes more extreme. In the
case when the vector points directly at a vessel boundary the algorithm will tend to fail. This
failure occurs because the minimum distance vector is nearly parallel to the axis normal to the
projection plane used for translations, resulting in near-zero translational adjustment vectors.
The algorithm assumes that at any given two points the vector between them is relatively parallel
to the centerline of the vessel. ..................................................................................................... 150
Figure C.3. An x-z (left) and y-z (right) view of the centerline algorithm results from scan H5. If
the final approximation of a centerline point differs greatly from the initial guess, the vector
between the previous centerline point and current centerline point becomes more extreme. In the
case when the vector points directly at a vessel boundary the algorithm will tend to fail. This
failure occurs because the minimum distance vector is nearly parallel to the axis normal to the
projection plane used for translations, resulting in near-zero translational adjustment vectors.
The algorithm assumes that at any given two points the vector between them is relatively parallel
to the centerline of the vessel. For this case the algorithm was able to recover, but it resulted in a
large erroneous region before recovery. ..................................................................................... 151
Figure C.4. An x-z (left) and y-z (right) view of the centerline algorithm results from scan K3. It
is unclear why the centerline fluctuated from its expected center, but it is clearly erroneous and
requires adjustment of operating parameters. ............................................................................. 151

xiii

KEY TO ABBREVIATIONS

AAA – Abdominal Aortic Aneurysm
COPD – Chronic Obstructive Pulmonary Disease
CT – X-ray Computed Tomography
DICOM – Digital Imaging and Communications in Medicine
ELG – Endoluminal Graft
EVAR – Endovascular Aneurysm Repair
FEA – Finite Element Analysis
FSI – Fluid Structure Interaction
G&R – Growth and Remodeling
HU – Hounsfield Unit
ILT – Intraluminal Thrombus
MPR – Multi Planar Reconstruction
MRI – Magnetic Resonance Imaging
NMR – Nuclear Magnetic Resonance
PWS – Peak Wall Stress
RF – Radio Frequency
SSE – Sum of the Square Error
STL – Standard Tessellation Language
SNUH – Seoul National University Hospital
US - Ultrasound

xiv

CHAPTER 1. Introduction
1.1 Motivation
1.1.1 Clinical Review of AAAs

The aorta is the major artery which blood circulates through from the heart. It can be
segmented into different sections such as the aortic arch, thoracic aorta, and abdominal aorta, or
ascending and descending aorta. The infrarenal aorta is the section of the abdominal aorta which
lies between the renal branches and the iliac bifurcation. The normal diameter of the infrarenal
aorta is approximately 2 cm. An aortic aneurysm is identified as an enlargement of the aorta
greater than 50% of the normal diameter, or in the case of the infrarenal aorta a diameter of
greater than 3 cm (Figure 1). Aneurysms affecting part of the vessel’s circumference are known
as saccular, while those affecting an area spanning the entire circumference are known as
fusiform. The vast majority of aortic aneurysms are AAAs (abdominal aortic aneurysms), and
over 90% of AAAs occur specifically within the infrarenal aorta [1–3]. AAAs are a serious
medical condition that when left untreated ends in vessel rupture with patient mortality rates up
to 95% [4], [5]. AAAs are responsible for at least 10,000 deaths annually in the U.S. alone [1],
[6], [7].

1

Figure 1. A sample CT scan of a patient diagnosed with an AAA highlighting a coronal slice
(left, upper-left), a sagittal slice (left, upper-right), a transverse slice (left, bottom), and the
intersection of all three (right). White arrows are used to identify the AAA. For interpretation of
the references to color in this and all other figures, the reader is referred to the electronic version
of this thesis The text is not meant to be readable this figure and other CT reference figures, and
is for visual reference only.

The highest risk group is men over 65 years of age, and AAAs are less common for
women than men and for individuals of black race/ethnicity [8]. The reasoning for this is not
clear, and studies are underway suggesting things such as hormones and different anatomically
influenced gender-based hemodynamics in the abdomen due to reproductive organs as possible
reasons for gender differences while differences between race/ethnicity are harder to discern.
Other risk factors include smoking, chronic obstructive pulmonary disease (COPD), high
cholesterol, atherosclerosis, and hypertension; approximately 50% of individuals with AAAs
also have hypertension [1], [2], [9]. Arterial tissue consists of three layers: the intima, media,

2

and adventitia. Loss of medial tissue and atherosclerosis, the stiffening of vessels due to internal
plaque buildup commonly from high cholesterol, are the most common causes of AAAs. AAAs
are usually asymptomatic until rupture, but can be associated with lower back and abdominal
pain as they progress; as an AAA’s size increases it may become palpable and/or tender [1], [10].
Rupture risk is the ultimate complication which clinicians consider for surgical decisions; rupture
risk for aortic aneurysms with diameters larger than 5 cm, 6 cm, and 7 cm are approximately
20%, 40%, and 50% respectively based on documented diameters immediately following rupture
[4], [9].

Currently medical professionals consider one geometrical parameter as the ultimate
indicator of rupture risk for a patient diagnosed with an AAA: the maximum cross-sectional
diameter. If the maximum cross-sectional diameter exceeds 5.0 cm or 5.5 cm an elective
surgical intervention, either open surgical repair or EVAR (endovascular aneurysm repair), is
suggested. Aortic replacement with a prosthetic graft was first introduced in 1952 by DuBost and
has been the treatment of choice until the introduction of transfemoral intraluminal graft
implantation in 1991 by Parodi [11], [12]. Implanting an ELG (endoluminal graft) promotes
beneficial tissue remodeling to decrease the aneurysm sac size by redirecting blood to flow
through the ELG to relieve pressure from the weakened aortic wall [13–15]. AAA research is
usually divided into two sub-fields: the study of the natural progression of AAAs pre-EVAR, and
the healing/recovery process of AAAs post-EVAR. This thesis will focus on pre-EVAR AAA
pathology and biomechanics.

3

The 5.0 cm to 5.5 cm maximum cross-sectional diameter criteria was developed to set a
clear boundary where the risk of rupture becomes appreciably greater than the risks involved in
surgery, which has mortality rates of up to 5.8% - 6% [16], [17]. Statistically significant studies
suggest continued monitoring for patients in lieu of early elective surgery carries a lower rate of
mortality [7], [17], [18].

This diagnosis can best be made through the attainment of an

Ultrasound (US), X-Ray Computed Tomography (CT), or Magnetic Resonance Imaging (MRI)
screening of the region of interest for an individual; it is estimated that physical examination
alone can miss more than one third of diagnoses [15], [19]. After the initial diagnosis is made,
patient monitoring is required for an expansion rate to be determined and for general AAA
stability to be assessed. From such scans parameters can be defined and analyzed based on
morphological features to develop more sophisticated measures and attainment algorithms that
account for more factors inherent to patient variance. Doing so will allow clinicians and patients
to make more informed treatment decisions and allow for the more accurate estimation of rupture
risk. Perhaps the least understood part of AAA pathology and rupture risk is the influence of an
intraluminal thrombus (ILT) though clinical studies have shown that 75% of AAAs contain ILT
[20]. An ILT is essentially a blood clot within the vessel which adheres to the intraluminal
surface (Figure 2). Formation of an ILT is influenced by hemodynamic and biochemical factors.

4

Figure 2. Transverse slice of a patient CT scan in the abdominal region (left), and a zoomed
section of a transverse slice in the abdominal region (right) focusing on AAA transverse crosssections exhibiting varying thickness of the arterial wall due to ILT presence shown by the
arrows. The intensity scale is adjusted on the right to show clear contrast between lumen, tissue,
and intraabdominal space.
1.1.2 Biomechanics Research of AAAs
1.1.2.1 Introduction

The primary goals of current research of small pre-EVAR AAAs are the accurate
quantification of morphological features and the accurate estimation of mechanical stress
experienced by the arterial wall. Doing so builds a greater understanding of AAA disease
progression on an individual level such that the estimation of rupture risk can shift from general
statistic probabilities to estimates tailored to patient specific concerns.

There exists much

controversy on the merits of surgery for small AAAs less than 5-cm in diameter. In general it is
assumed the risks of surgery outweigh the risk of rupture for such cases, but rupture of small
AAAs is still responsible for high mortality rates [21], [22]. Understanding patient specific
biomechanics is most beneficial for treatment decision making with regards to small AAAs

5

where rupture risk may exceed surgical risks but the AAA maximum diameter criterion suggests
otherwise.
1.1.2.2 Arterial Wall Mechanics

Technical studies tend to focus on the geometrical analysis and finite element simulation
of 3-D models generated from the segmentation of 2-D multi planar reconstructions (MPRs) of
CT scans, which introduces the additional step of 3-D model generation. Many research groups
have reconstructed 3-D AAA geometries from CT scans with the aim of evaluating wall stress
through conventional finite element analysis (FEA). Giannoglou et al. evaluated peak wall stress
(PWS) and found good correlation between PWS and mean centerline curvature, maximum
centerline curvature, and maximum centerline torsion, but not with maximum diameter, the
standard clinical correlation [23]. They based the surrounding tissue model on a uniform wall
thickness derived from Raghavan and Vorp, whom are leading experts on AAA tissue wall
material mechanics [24].

There have been significant developments for computational biomechanical models.
Doyle et al. evaluated wall stress via FEA using validated realistic nonlinear material properties
in an attempt to correlate wall stress to diameter and centerline asymmetry, finding that posterior
wall stress increases with anterior centerline asymmetry by 38% [25].

Raghavan et al.

incorporated a nonlinear biomechanical computational AAA tissue model in place of a
conventional model into the reconstructed geometries to evaluate wall stress distribution, finding
the PWS to be posterior in all cases and showing that AAA volume is a better predictor of PWS
than maximum diameter, and therefore possibly a better indicator of rupture risk [26]. Wang et
al. evaluated the effect of ILT on wall stress using a nonlinear large deformation algorithm
6

finding that PWS was reduced in the presence of ILT. The results showed an association
between increase in ILT presence and decreases in peak wall stress of 6% - 38% from a range of
2

2

2

30 – 44 N/cm to a range of 28 – 37 N/cm , compared to 12 N/cm for a healthy aorta [20].
Georgakarakos et al. evaluated PWS in the presence and absence of an ILT and its relation to
geometric parameters such as torsion, tortuosity, and mean curvature. They found a positive
correlation between PWS and ILT volume, and in the presence of an ILT found significant
correlation with the degree of centerline tortuosity and maximum diameter [27]. Speelman et al.
focused on preserving wall calcifications, or mineralization of sections of the arterial wall,
throughout 3-D model reconstruction and analyzing their effect on wall-stress [28].

Wall

calcifications are an interesting part of AAA pathology which develop from unique biochemical
conditions, but as they were rare in our dataset they were not modeled and considered. FEA
methods quantifying wall stress have proven more effective than the maximum diameter
measurement alone at predicting rupture risk. While it is clear AAA maximum diameter is not
the measurement best correlated to rupture risk, it is not completely clear what is the best
measure either.
1.1.2.3 Morphological Predictors of Arterial Wall Stress

There are a number of limitations when analyzing the clinical methods of rupture risk
estimation. Simple parameters are not perfect measures and do not necessarily account for
patient-specific information.

This has prompted researchers to investigate more sophisticated

parameters, but such investigations are also confronted with limitations. For example, it is
known that the assumed 2 cm normal infrarenal aorta diameter can actually vary from 1.5-cm to
2.5-cm between individuals, but that a difference of 1 cm diameter can indicate much different

7

rupture risks [4], [9], [13]. Advanced geometrical parameter definitions & attainment algorithms
will lead to more accurate risk predictions, increased sensitivity to patient-specific concerns, and
form a foundation for defining more reliable indicators of rupture risk. Much research has been
done to identify suitable geometric parameters to this end.

Martufi et al. evaluated geometric indices and validated an in-house automated software
package, VESSEG, suggesting use of a geometrical index incorporating multiple morphological
features [13]. Shum et al. evaluated geometric indices and regional variations in wall thickness
based on novel segmentation algorithms [29]. A number of fluid structure interaction (FSI)
simulations have also been performed in an attempt to link wall stress distribution to
hemodynamic factors. Xenos et al. analyzed FSI with a focus on ruptured aneurysm [30].
Bluestein et al. analyzed FSI in the presence of ILT [31]. Fillinger et. al determined the wall
stress distribution as a result of geometry and blood pressure using a hyperelastic nonlinear
model [32]. Wolters et al. developed a novel patient-specific mesh generation algorithm to
generate FSI meshes based on CT scan data [33]. Leung et al. compared a FSI model to a solid
stress model [34]. Papaharilaou et al. performed a decoupled FSI-FEA simulation for estimating
wall stress [35].

While not directly studying AAAs and CT scans, other works have contributed to the
general discussion surrounding AAAs. Borghi et al. reconstructed 3-D geometries of thoracic
aortic aneurysms from MRI by scanning the same region using different parameters optimized
for resolution of the aneurysm lumen, thrombus, and wall to reconstruct a 3-D model detailing
each component accurately [36]. Steinman et al. reconstructed a 3-D large intracranial aneurysm

8

model and performed a finite element computational fluid dynamics simulation of the 3-D
pulsatile velocity field [37].

Current work makes frequent use of isolated patient-specific 3-D AAA geometries
reconstructed from CT scans [13], [20], [23], [25–35]. Patient-specific 3-D aneurysm geometries
are also useful for the study of thoracic aneurysms and large intracranial aneurysms [36], [37].
Conventional FEA wall stress studies for idealized and patient-specific geometries are more
predictive of rupture than AAA maximum diameter, but many of these studies assume arterial
tissue to be isotropic and linearly elastic as well as disregard vertebral interaction of an AAA
during disease progression [23], [27], [28], [38–40]. The desire in the research community to
use patient-specific models is clear, and the inclusion of longitudinal data and surrounding tissue
would serve to directly enrich the findings of any of these existing works.
1.1.2.4 Intraluminal Thrombus

It has been shown that there is a correlation between an AAA diameter increase and ILT
growth [41], [42]. It has been more recently suggested that there is not only an association, but
that the ILT plays an active role in AAA pathology via mechanical stress shielding and pressure
reduction [20], [43]. Such findings are important given that the basic mechanical definition of
rupture of an AAA is when the wall stress at a surface exceeds the wall strength at that surface.
Hinnen et al. extracted actual ILT samples and implanted them in model aneurysms composed of
rubber and paraffin. The systolic circulatory pressure was measured as well as the systolic
pressure at various depths within the ILT from the wall and it was found that the ILT served to
reduce pressure as distance to the wall decreased. At 1 cm, 2 cm, and 3 cm from the wall
pressure ratios between the systolic circulatory pressure and systolic ILT pressure were 0.9 ±
9

0.09, 0.86 ± 0.10, and 0.81 ± 0.09 respectively [43]. The results carried two-fold importance:
they experimentally verified the theory of stress-shielding presented by Wang [20], and they
highlight the importance of considering ILT for FSI. The proper inclusion and modeling of ILT
is another source of enrichment for current works, but there are significant challenges in such
incorporation. Most technical studies of patient-specific AAAs reference CT data, which does
not offer clear boundaries between the vessel wall and ILT. Idealized geometries often assume a
uniform wall thickness ranging from 1.2 mm to 1.8 mm. High resolution CT scans can offer
sub-millimeter resolution on the order of at least 0.5 mm / pixel, which does not offer a clear
benefit over idealized thicknesses when other sources of error are taken into account. As ILTs
become increasingly larger their accurate modeling becomes more realistic and worthwhile,
however, as the error due to CT resolution is absolute and only a concern for low thicknesses.
1.1.2.5 Discussion

The weighing of rupture risk against surgical risk is the bottom line for treatment
decision, and therefore a proper understanding of long-term AAA pathology and biomechanics is
paramount to accurately estimating rupture risk. The last two decades have been a time of
significant advancement for such understanding, leading to FEA / computational biomechanics
studies utilizing patient-specific AAA geometries.

Computational biomechanics considers

vascular adaptations associated with AAA pathology through growth and remodeling (G&R)
mathematical models.

Conventional FEA is useful for studying present and short term

deformations of patient-specific AAA geometries, but its inherent exclusion of vascular adaption
limits its potential for long-term studies. The stress analyses of patient-specific AAA geometries
producible through computational biomechanics have proven more effective than standard
clinical measures at predicting rupture risk thereby affirming its potential usefulness, and
10

perhaps necessity, as a rupture risk assessment tool. One significant drawback hindering the
incorporation of these methods clinically is their inability to explain why rupture occurs most
often at the lateral / posterolateral regions of an AAA. There is no doubt that future expansion of
studies not already doing so would benefit from the inclusion of patient specific geometries. It is
clear from these studies that a generally accepted set of guidelines and methods exist for
segmentation and model generation. For this study we additionally include the generation of the
vertebrae model adjacent to the AAA for global registration.
1.1.3 Surrounding Tissues

The pathology of AAAs is such that enlargement of the aorta is associated with disease
progression.

As an AAA enlarges it becomes subject to increased reaction forces from

surrounding tissue, especially from the vertebrae of the spine. The maximum diameter criterion
may serve to distinguish stable and unstable aneurysms to an extent, but cannot necessarily
explain the degree of stability. Back pain can be associated with AAAs due to aneurysm
expansion exerting pressure on the lumbar nerve roots as a result of vertebral contact. It is also
well known that there is an association between the presence of an AAA and vertebral erosion as
well as evidence shown for ruptures tending towards the posterolateral surface of the aorta.
There are not many studies that have been done considering the growth and remodeling
interactions of an AAA, the lumbar region of the spine, the ILT, and surrounding tissues. Vorp
et al. suggested that the “limitation of posterior expansion caused by the vertebral column might
result in preferential anterior expansion of the aneurysmal wall and an asymmetric
configuration” [44]. Watton et al. developed the first mathematical model to predict aneurysm
growth including spinal contact based on this suggestion, but it was not paired with a

11

longitudinal study or additional tissue constraints [45]. In this study it is shown that different
patients display different degrees of interaction with surrounding tissues. For a subset of the
patients, however, a marked effect on AAA geometry is observed.

Clearly defined

morphological parameters describing these interactions would be of significant value.

The vena cava is the major vein for returning blood back to the heart. The superior vena
cava returns blood from the upper half of the body to the heart while the inferior vena cava
returns blood from the lower half of the body to the heart. The inferior vena cava, spine, renal
system, and aorta interact as shown (Figure 3).

Superior
Mesenteric
Artery

Left Renal Vein

Left Kidney

Right Kidney

AAA

Figure 3. Near the lumbar vertebrae (white) the inferior vena cava (blue) runs opposite to the
aorta (red) on its right lateral side. The inferior vena cava and the aorta are connected, most
notably for AAAs, to the renal system (purple). The left renal vein crosses anteriorly over the
aorta inferior to the superior mesenteric artery.

12

All renal arteries and renal veins segment and terminate at the respective kidneys. In
healthy individuals the two vessels do not have substantial physical interaction, but in the
presence of an AAA the inferior vena cava presents an asymmetric right-sided growth barrier to
AAA geometrical evolution. Mean arterial blood pressure, which has a pulsatile flow, in healthy
individuals usually ranges from 90 – 100 mm Hg compared to the non-pulsatile venous side of
the circulation which is approximately 84 mm Hg less; mean venous blood pressure ranges from
approximately 10 mm Hg at the venules to approximately 0 mm Hg at the heart [1]. Blood flow
in the inferior vena cava is strongly influenced by the intraabdominal pressure (IAP), which
typically ranges from 2 – 10 mm Hg [1]. When considering the influence of the inferior vena
cava as a growth barrier of an AAA the mean venous pressure is thus far less influential than the
biomechanical resistance the tissue itself provides.

Each kidney of the renal system is sheathed by a fibrous external capsule and surrounded
by fatty connective tissue which, in combination with attached vessels and fascia, holds each
kidney in place [1]. The kidneys have some freedom of movement to react to movement of the
diaphragm, but otherwise remain relatively stable. As an AAA progresses it is possible for renal
arteries or veins to become stressed or distended, though migration of the vessels is resisted by
the other tissues serving to anchor the kidneys in place, most notably the left renal vein crossing
anteriorly over the aorta.

The spine acts as a rigid boundary for AAA growth, and can therefore affect the
remodeling of the aorta [44–46]. The interaction between vertebral contact & ILT position is
also extremely important. It is believed that the circumferential position of an ILT can lead to

13

substantially different AAA growth & remodeling, especially with regards to accelerated tissue
degradation. An anteriorly growing ILT may grow relatively unconstrained, while a posteriorly
growing ILT presents additional considerations. Such an ILT could cause a general vessel
centerline migration of the associated AAA away from the spine, as well as expose the posterior
side of the associated AAA to compressive forces of the spine as a growth barrier. It is known
that the arterial wall becomes thinner and weaker at areas contacting ILT, so it should be
considered whether such an area is also contacting the spine and therefore under additional force.
Additionally such a situation favors curvature flattening of the vessel at the region of spinal
contact and such a geometrical profile changes the stress experienced by the vessel, especially at
circumferential regions of significant and sudden curvature change.
1.2 Objectives

This thesis will detail an investigation of AAA surface evolution for a longitudinal CT
scan dataset progressing from qualitative to quantitative methods.

The identification of

influential morphological features coupled with automated algorithms to obtain parameters
representing those features will lend insight into the link between AAA pathology and AAA
biomechanics and thereby provide a strengthening tool for computational biomechanics. In the
selection of the parameters it is important that they should be intrinsically understood by medical
professionals as integrated study of AAA pathology and biomechanics requires clear discourse
between medical and research professionals. It is known that some morphological features arise
from contact between an evolving AAA and vertebrae, but the impact of these effects on rupture
risk is not fully known. The geometrical parameters investigated in this work will reflect
interactions with surrounding tissues, most notably the lumbar vertebrae. A longitudinal study

14

serves to further justify the suitability of selected morphological parameters by investigating
their respective evolution throughout AAA progression. Therefore the objectives of the study
are as follows:

-

Process a collection of CT scans to obtain fully segmented and registered 3-D models of
the vertebral and arterial surfaces.

-

Identify morphological features from the reconstructed CT scans and perform analysis
akin to clinical measurements.

-

Perform a relatively defined parameter study based on arterial centerlines generated by a
developed program to assess morphological feature evolution during AAA progression.

Chapter 2 will introduce a working knowledge of medical image modalities such CT,
MRI, and US. The advantages and disadvantages of each modality will be discussed with
respect to AAA diagnosis, monitoring, and surgical planning. The processing methodology
employed for our dataset will be covered in which 3-D models are reconstructed from CT scans.
The sum of 3-D models for a patient form a longitudinal model set. Three dimensional models
of the AAA luminal surface, AAA tissue surface, and lumbar vertebrae will be developed
following established methods of image segmentation. Models belonging to longitudinal subsets
of the data will additionally undergo an image registration process to find the spatial transform
that maps their positions and orientations with respect to the lumbar vertebrae, which is assumed
relatively unchanging over time. Image registration allows for an accurate investigation of the
true spatial differences between scans at different times.

15

Chapter 3 will examine the CT scan slices concurrently with the generated 3-D models to
identify morphological features influenced by surrounding tissue. It is theorized that certain
morphological parameters may be more suitable for indicating or estimating rupture risk than the
currently used maximum diameter criterion as such features may be indicative of unique sources
of increased wall stress observable through FEA.

Each CT scan will be imported to and

visualized through biomedical imaging software for analysis and subsequent determination of
suitable morphological parameters. Measurements of the maximum tissue and luminal diameters
will be taken directly from the CT scans according to established clinical methodology.

Chapter 4 will cover parameterization of the surface.

The centerline generation,

smoothing, and parameterization of the surface will be described in detail.

After

parameterization of the surface is complete, the method for extracting morphological
measurements from the surface will be covered. The centerlines serve as a reference, or an axis
to parameterize the surface to, for each model which allows for parameterization of the surface
and the relative definition of geometrical parameters.

It is theorized that interaction with

surrounding tissue, particularly the lumbar vertebrae, vena cava, and renal system, plays an
important role in the evolution of the AAA surface geometry.
1.3 Contributions

It is believed that the evaluation of morphological parameters throughout AAA
progression that can be used to quantify AAA surface evolution with respect to surrounding
tissues is the pivotal next step in advancing the general understanding of AAA biomechanics.

16

Through collaboration with Dr. Whal Lee at Seoul National University Hospital (SNUH) in
South Korea we have obtained a high resolution CT scan dataset for this purpose.

The influence of surrounding tissues is understudied in current works. Immediately upon
visually reviewing the CT scans it could be seen that surrounding tissues, particularly the spine,
were significantly interacting with the AAAs during their progression. A FEA or computational
biomechanics study would typically involve the internal arterial pressure and vessel pre-stretch
as the sole mechanical driving forces of surface evolution, but the addition of a growth barrier
such as the lumbar vertebrae would introduce significant reaction forces. Considering the
lumbar vertebrae would change the outcome of such simulations and offer more insight as to the
complex biomechanics in effect under in vivo conditions. Consideration of the mechanical
influence of other tissues such as the inferior vena cava, adipose tissue, and connective tissue
would also lend insight into the link between AAA pathology and biomechanics, but the spine is
focused upon for its clear contribution.

The dataset contains a total of 51 CT scans from 11 patients with a mean surveillance
interval of 373 days between longitudinal scans allowing for the unique investigation of AAA
surface evolution. Such data allows not only for the investigation of morphological parameter
values, but their evolution as well.

Studies including longitudinal data sets are seriously

underrepresented as access to such data collections is rare, especially at high-resolutions.
Typically 3-D meshes developed from patient-specific data are based on single CT scans of
individual patients with no additional follow-up CT images for study, unlike access to a
longitudinal dataset provides.

The major contribution of this study is to showcase the

17

importance of considering surrounding tissue, particularly the spine, when performing stress
studies and geometrical parameter indexing of ideal and patient-specific AAA geometries. The
developed centerline algorithm is proposed to be superior to standard algorithms which calculate
the 2-D transverse polyline centroids to develop a centerline [13], [23], [25], [29], [47–50].

18

CHAPTER 2. Biomedical Image Processing: Segmentation and Registration
2.1 Introduction
2.1.1 Biomedical Image Processing Background

This chapter serves to highlight the processes involved in fully visualizing and modeling
the anatomical structure of an in-vivo patient-specific AAA. Medical professionals performing
diagnosis, continued monitoring, or surgical planning of patients require clear visualization of
the AAA geometry. Research groups performing morphological feature surveys, FEA, or growth
and remodeling simulations of patient-specific AAA geometries need to have quantitative
information about the geometry. In standard clinical practice the qualitative imaging resulting
from US is most commonly used for diagnosis and continued monitoring. Use of CT and/or
MRI scanning is typically reserved for surgical planning, though they offer the ability to record
quantitative geometry data for delayed decision making and analysis [51–53]. From a qualitative
image, whether from US, CT, or MRI, an accurate production of quantitative geometry
measurements can result, though the AAA maximum external tissue diameter is the standard
measure of interest. In clinical practice measurement is usually performed through the use of
calipers on real-time US imaging [51–53].

Technical studies trying to quantify various

parameters tend to focus on 3-D analyses, which involve additional processing of the data,
usually CT scans, before analysis (Figure 4).

19

Data
Obtain CT scans

Render data anonymous

Import to software

Segmentation
Vertebrae

Lumen

Tissue

3-D Modeling
Generate 3-D models

Remesh

Register

Technical Computing
Parameterization

Simulation

Analysis

Figure 4. A process diagram highlighting key elements of the biomedical imaging process from
data obtainment to generating technical computing models.
A brief review of the relevant biomedical imaging modalities is first presented. This
serves to justify the use of CT scans seen by most studies and introduce key terminology
associated with each modality. After this introduction the methods utilized for segmentation and
3-D modeling / image registration are covered.
2.1.2 Biomedical Imaging Modalities

To properly diagnose, perform continued monitoring, and perform surgical planning for
an AAA, noninvasive imaging modalities are employed clinically; Acharya et al. provides an
excellent resource to become familiar with noninvasive imaging modalities from which much of
this summary is derived [54]. Noninvasive imaging modalities can visualize internal anatomical
structures of the human body accurately and quickly, even providing real-time imagery in some

20

cases. With respect to diagnoses of AAAs US, MRI, and CT are the most widely favored
modalities. A summary of each will be presented and discussed.

US refers to sound waves with high temporal frequencies.

US imaging involves

producing US signals, directing them to a particular area in the body with handheld probes, and
receiving the signal at an appropriate detector array. Image reconstruction takes places in realtime on a monitor. As the signal passes through different layers of tissue and encounters tissue
layer boundaries partial reflection and partial transmission of the signal occurs. There are a wide
variety of electronic transducer arrays that can be employed to acquire different dimensionality
and quality of this qualitative imagery, but the same basic principles apply for each setup. Many
setups allow for the 2-D real-time visualization of structures within the body, which is
particularly useful for cardiology. 3-D US was patented in 1987 but still sees contention for
widespread clinical use mainly due to cost, training, and data storage [55]. US imaging is
capable of not only acquiring qualitative imagery, but quantitative diagnostic information.
Doppler Echocardiography, for example, can provide vector measurements of blood flow
direction and magnitude. The major advantage of US is real-time feedback and equipment that is
compact, lightweight, and of comparatively low-cost which also requires less specialized
personnel compared to other imaging modalities’ equipment.

MRI is based on Nuclear Magnetic Resonance (NMR), which is a phenomenon
associated with physical chemistry. NMR was discovered in 1946 by Bloch and Purcell, who
were later awarded the Nobel Prize in 1952 for the discovery, which was coincidentally the same
year that aortic replacement surgery was introduced by DuBost [11], [54]. For an MRI scan a

21

magnetic field is generated and maintained enclosing the scanning area. MRI thus requires
patients to be subject to extremely powerful magnetic fields within a gandry during the scanning
process, which makes factors such as implants, pacemakers, and claustrophobia important to
consider as contraindications. Radio Frequency (RF) energy is introduced into the scanning area
to excite nuclei within the field, which releases RF energy as the MR signal to be recorded. The
ideal element to perform MRI with is Hydrogen, which is also the most common element in
biological tissue. MRI is capable of collecting images at any plane through a volume and is
excellent for soft-tissue resolution. Automatic Multi Planar Reconstruction (MPR) of the image
takes place after the scan data is reconstructed allowing for 3-D viewing of internal anatomical
structures. MRI is extremely expensive and not widely available, rendering it an uncommon
choice for AAA diagnosis and monitoring protocols. Usually MRI is only utilized for AAAs if
an individual requires a CT scan but has contraindications for CT scanning.

CT scanning utilizes X-rays to acquire diagnostic images. A conventional radiograph
produces a single 2-D image by placing the scanning area of a patient between an X-ray emitter
and X-ray receiver. X-rays sent through the scanning area are linear and lose energy as they pass
through tissue. The resultant X-ray energy is recorded for the field of view. The result is a
radiograph, which is only a 2-D projection of a 3-D scanning area. CT, on the other hand, is
capable of returning 3-D information as a series of 2-D slices based on the same principles. An
X-ray source on one side of the scanning area propagates X-rays to an X-ray detector on the
opposite side of the scanning area parallel through the axial plane. The detector is capable of
electronically storing the information. Relative motion is applied to the system between the
scanning apparatus and the scanning area such that the X-ray source and detector remain in the

22

axial plane, but at some fixed angle from their previous positions. The emission and detection
process is repeated according to the parameters set by a given scanning algorithm. A total of
180° of rotation is required to fully image the plane in two dimensions. There exist many
reconstruction algorithms, but in essence they all manage to reconstruct 2-D tissue density
information from the saved collection of X-ray projection information. The 2-D imaged plane
represents a finite 2-D matrix of tissue density values known as an axial or transverse slice.
Slices have a thickness dimension, however, as well as a finite resolution as a result of the
electronic digitization. Therefore each element of a slice matrix actually represents a voxel, or 3D pixel, with dimensions fixed by the scanning algorithm and reconstruction algorithm. Some
algorithms are inherent to the scanner used in that they produce identical resolution of data every
time, but some scanners exist such as the Dynamic Spatial Reconstructer that are capable of
dynamic algorithms to meet spatial resolution demands. CT scans are the method of choice for
surgical planning of an AAA after US monitoring has identified the need for treatment [51–53].

There are advantages and disadvantages to US, MRI, and CT with respect to AAA
diagnosis. US has the advantage of mobile compact equipment such that patient comfort is
minimally affected whereas MRI and CT both require the patient to remain motionless during
scanning, in addition to often presenting a claustrophobic environment which dissuades some
patients from being comfortable with the process. US also presents both a low cost up-front for
the medical provider and for the patient relative to MRI and CT which present much higher
costs. US is additionally able to visualize motion of anatomical structures through real-time
imaging and extract quantitative information such as blood flow direction and magnitude. If a
medical professional wishes to quantify the maximum tissue diameter of an AAA for diagnosis

23

or continued monitoring, US provides a quick, low-cost, noninvasive method of accomplishing
that goal, and studies have shown it is only marginally less reliable than CT for maximal
diameter measurements [56]. Motion of the vessel from the cardiac cycle introduces intraoperator variability in maximal AAA tissue diameter measurements from US, however, and there
is no standardized method to reduce this variability [57]. US sees widespread use in the United
States as a diagnostic and monitoring tool as it provides clear information about AAA size.
Obesity and the presence of bowel gas can render US diagnosis of an AAA less effective,
however. AortaScan, a portable 3D ultrasound device which measures abdominal aortic tissue
diameter automatically, was investigated against conventional CT scanning for the diagnosis of
AAAs in 44 known AAA cases. Conventional CT diagnosis by trained screening personnel was
able to correctly diagnose 43 / 44 cases, while AortaScan missed the diagnosis in 8 / 44 cases
and reported 13 / 44 false positives [58]. This may not be a fault of US, but rather illustrates the
shortcomings of automated diagnosis vs. the judgment of trained screening personnel. US is an
ideal choice for visualizing small stable AAAs, but becomes inadequate when advanced
investigation and/or surgical treatment planning is required.

CT and MRI are also effective at diagnosis, but considered unnecessary for small stable
AAAs due to their higher cost, decreased patient throughput, lessened availability, and principle
concerns. CT involves X-rays which are inherently harmful to tissue, and though permissible
levels have been set, exposure should be limited as much as possible. In the case of CT for
AAAs an intravenous dye is injected into the aorta to render the luminal volume easily
differentiable from surrounding tissue. In some cases patients may have an allergic reaction or
experience irritation of the kidneys due to the dye, especially patients with chronic kidney

24

concerns. For MRI the necessity of strong magnetic fields defines absolute contraindications for
scanning. Such magnetic fields could disrupt the proper function of a pacemaker for example, as
well as affect other implants in undesirable or fatal ways. US is the preferred clinical tool for
initial diagnosis and continued monitoring due to low cost, wide availability, compact and
mobile equipment, and minimal patient discomfort.

If AAA disease progression warrants

additional investigation and/or surgical planning CT scans with intravenous dye then become the
preferred scanning method due to the clear and absolute depiction of volumetric information.
CT scans are otherwise considered unnecessarily time and cost consuming for small stable
AAAs even though they offer superior qualitative imagery [56].

If a patient has

contraindications for CT scanning, MRI is then an adequate alternative for surgical planning.
This review serves to illustrate the rarity of a longitudinal CT scan collection of small AAAs and
the quality of information that can be extracted.
2.2 Methods
2.2.1 Data

A biomedical software, Mimics® (Materialise, Leuven, Belgium), is used to process the
CT scans, which are axial/transverse slice stacks [59]. The software interface allows for viewing
of up to 3 voxel domain, or slice, windows simultaneously through MPR and a 3-D rendering
window. MPR is a general term that refers to the interpretation of additional planes based on
known planes. Standard MPR results in the additional view of the sagittal and coronal planes
from the transverse planes, but other methods exist that generate curvilinear planes based on
arterial centerlines for example. Other commercial biomedical imaging software packages exist,
as well as some in-house programs such as that used by Martufi et al. (VESSEG v.1.0.2,
25

Carnegie Mellon University, Pittsburgh, Pennsylvania) [13]. For the current study a collection
of high resolution CT scans of patients with AAAs was obtained through collaboration with Dr.
Lee’s clinical group at SNUH according to an approved institutional IRB protocol; any
identifiable information that could be traced back to the patient was removed prior to our
obtainment of the dataset.

Supplemental information was also provided: blood pressure

measurements, age, and gender for each scan. Typically clinical small AAA monitoring is done
using US, and a collection of CT scans of this caliber for small AAAs is rare, especially in the
United States. Each scan features a contrast agent in the form of an intraluminal dye which
renders the luminal volume easily differentiable from nearby tissue.

We obtained a total of 51 sets of CT images representing the longitudinal study of 11
patients termed patients A - K, with 2 to 9 sets of images from each patient and a mean
surveillance interval of 373 days. The combination of patient ID and scan number will be
adopted for this study such that the first scan from patient G would be known as scan G1, the
second scan from patient G would be G2, etc. Of these 51 scans were included 3 post-EVAR
scans, 2 thoracic aortic aneurysm scans, and 1 scan which failed to capture the AAA region
completely.

Discounting the aforementioned scans reduces the dataset to a total of 45 sets of CT scans
representing the longitudinal study of 10 patients, with 2 to 7 sets of images from each patient
and a mean surveillance interval of 401 days. Each CT scan has dimensions of 512 px × 512 px
in the transverse plane and a 0° gandry tilt. The voxel resolution varies between scans as a result
of the scanner brand/model, scanning algorithm, and reconstruction algorithm employed. The

26

total number of slices varies as well but usually includes the entire torso of the patient,
presumably allowing for the additional diagnosis of thoracic aneurysms. It is possible for the
geometry data to be scaled, but scaling of the entire dataset was at 100% upon obtainment and
never altered or further considered henceforth. Patient information and a statistical summary of
the scans are shown by Tables 1 and 2. More detailed information for every scan is located in
Appendix A (Tables A.1 and A.2).

Table 1. Patient age, number of scans, and sex after discounting irrelevant scans.
Patient
A
B
C
D
F
G
H
I
J
K

# Scans
2
3
2
3
6
6
7
6
5
5

Age
68
71
69
63
78
65
68
66
54
62

Sex
M
M
M
F
F
M
M
M
M
M

Table 2. Statistical summary of the CT dataset discounting irrelevant scans. The surveillance
interval refers to the time between individual scans, such as between J1 and J2 or J4 and J5. The
total study length refers to the total time between the first and last scan of a patient, such as
between J1 and J5.
# of Patients
Mean age
Minimum # scans
Maximum # scans
Mean # scans
Minimum surveillance interval
Maximum surveillance interval
Mean surveillance interval
Minimum total study length
Maximum total study length
Mean total study length

27

10
66.4
2
7
4.5
84
1713
401
182
3325
1502

days
days
days
days
days
days

Some scans were discounted from this study for different reasons, which are important to
justify before moving on. Scan F1 is not focused on the abdominal region, and the scan is
therefore fully discounted (Figure 5). Scans H8 – H9 are post-EVAR, and the presence of an
ELG is clearly present (Figure 6). The focus of this study is AAAs pre-EVAR, and these scans
are therefore fully discounted. Scan K6 appears to be post-EVAR but it is difficult to discern as
no ELG is present, though it is possible a full aortic vessel graft is present (Figure 7). While the
cause is unclear, it is clear that scan K6 does not represent natural progression of an AAA and
the scan is therefore fully discounted. Scans E1 – E2 clearly feature a thoracic aneurysm only
and the focus of this study is AAAs, so the scans discounted (Figure 8). Figures 5 - 8 clearly
describe the reason for discounting the aforementioned scans and are detailed below. These
figures feature 3-D results from the methodology covered in Chapter 2 for the sake of visual aid,
but the same conclusions are clear from visual analysis of the CT scans without employing the
methodology covered in Chapter 2.

Figure 5. A 3-D left side view of F1 after initial processing. The lower scan boundary is
equivalent to the bottom of the image which clearly does not contain the abdominal cavity nor
the full presence of an AAA. The scan F1 was therefore discounted from this study on the basis
of insufficient data.

28

Figure 6. A 3-D transparent frontal view of the tissue surface (green, left) and the ELG (blue,
left), and the luminal surface (red, right) and the ELG (blue, right) of scan H8. The presence of
an ELG indicates this to be data from a post-EVAR scan, which is not the focus of this study.
The scans H8 – H9 were therefore discounted from this study on the basis of irrelevant data.

Figure 7. a 3-D view overlaid on a respective CT scan coronal plane for scan K4 (left) and for
scan K6 (right). The shape of the AAA from scan K6 does not represent natural AAA disease
progression; it instead more closely represents a post-EVAR geometry. The focus of this study
is pre-EVAR AAAs and K6 is therefore fully discounted.
29

Figure 8. A 3-D view overlaid on a respective CT scan coronal plane for scan E1. A thoracic
aneurysm is present, but no AAA is present. The focus of this study is AAAs, so the scans are
discounted.
2.2.2 Segmentation

The software interface allows for viewing of up to 3 voxel domain windows and a 3-D
rendering window. The voxel domains, or slices, can be viewed from the sagittal, coronal, and
transverse planes expectedly. Each voxel represents a finite volume with a single intensity value
based on the HU (Hounsfield Unit) scale, a measure of tissue density. A brief summary of some
core tools of the biomedical imaging software package is covered first. The application of these
tools is then described in some detail to produce specific voxel domains, more commonly
referred to as masks, representative of desired anatomical structures through a process known as
segmentation.

Thresholding refers to setting an upper and lower bound based on voxel intensity values,
in this case representative of tissue density, and separating voxels within the bounds from voxels
30

outside the bounds to binarize a subset of the data. Essentially a range is identified on the tissue
density histogram associated with the sum of all the voxels, of which many preset ranges exist, to
develop a mask. Cropping is a standard operation for any type of image processing in which
new borders are defined within an image and segments of the image outside these bounds are
removed. In this biomedical imaging application it works no differently, save for the operation
is performed in 3-D to voxels instead of 2-D to pixels. Region Growing is a unique operation
which can greatly increase overall processing speed by examining the 2-D or 3-D connectivity of
voxels. A single voxel is selected as a source and only voxels connected by the specified
dimension of connectivity which are already within the mask remain within the mask while all
others are removed. This can serve to eliminate disconnected voxels from a mask acting as
background noise as well as more advanced applications. Morphologic Operations define a set
of advanced tools that utilize generated masks created through other means. An existing mask
can be dilated, eroded, opened, or closed in three dimensions automatically by a set voxel value.
For example, a voxel cube with dimensions 3×4×5 dilated by 1 voxel would become a voxel
cube with dimensions 5×6×7. Additionally this tool can limit operations to only occur within a
specified mask, effectively applying user-specified bounding constraints.

The first step to 2-D mask generation of the arterial lumen & vertebrae is the proper
selection of a HU intensity value range, better known as thresholding. When selecting the
threshold limits it is important to consider the fact that any range could be tightened following
any future operations if deemed necessary later, but there is no way to widen the range after
other operations. The arterial lumen contains a contrast agent which renders it easily identifiable
for a select threshold for each scan, usually > 226 HU. Applying a threshold selects every pixel

31

matching the criteria and assigned them to a mask. The threshold range for the vertebrae, on the
other hand, is not always straightforward to assign.

It is well known that bone density

differences exist between individuals for a variety of reasons such as age, sex, nutritional intake,
and lifestyle. The threshold for the vertebral region of interest therefore has to be adjusted from
a base range of > 226 HU on a scan-by-scan basis. There is additionally a fill holes option for
thresholding that is used to fill empty voxels within the region of interest.

The goal of

processing the data within Mimics was to export surface information, which would have been
erroneous with internal surfaces caused by holes within regions of interest.

Upon completion of thresholding a number of tools are used. First and foremost the
region of interest is specified by a virtual bounding box and cropped. Connective pixels are then
identified between regions of interest and virtually severed; for example, the intercostal arteries
of the aorta are virtually severed from the spine. This process is done with manual editing tools
in both 2-D and 3-D, which are time-intensive user-based subjective operations.

Upon

successfully severing connective pixels the region growing tool can be used on the region of
interest to remove background artifacts. Employing this process for the vertebral and luminal
regions results in two masks which represent the internal anatomical structures. An example of a
segmented mask of lumen in a coronal and transverse slice is shown by Figure 9.

32

Figure 9. A segmented mask for the arterial lumen (red) and arterial tissue (green) of scan B2
shown for a coronal slice (left) and transverse slice (right). The transverse slice (right) highlights
the rectangular nature of voxels which occupy a segmented mask. The arterial tissue is seen to
contain a wide region on the upper-left, or left anterior, side indicating the presence of an ILT.

The mask generation for the tissue domain is much different than for the luminal and
vertebral domains, and also somewhat depends on the prior generation of each. The vertebral
and luminal domains each occupy a distinct region of the tissue density spectrum of the entire
image. The arterial tissue, on the other hand, shares a widely populated region of the tissue
density spectrum with other soft and/or muscular tissues. Intuitively, however, it is known that
the tissue domain exists immediately adjacent to the luminal domain. Thus a working domain
can be defined based on the tissue density and the tissue thickness at its thickest point through
morphology operations.

These operations are based on the existing luminal model and

discussion with Mimics® representatives reinforced the advantages of this method as it is
repeatable and reliable.

The preset thresholding range for the arterial tissue was > -25 HU. The region of interest
was then further defined by a cropping operation.
33

The current tissue mask thus contains

vertebral and luminal voxels. The vertebral mask was subtracted from the tissue mask via a
Boolean operation to quickly and efficiently remove voxels pre-defined as vertebral voxels. The
tissue mask now contains luminal voxels as well as other voxels > -25 HU. A new mask
containing the resulting voxels from dilating the luminal mask by a user-defined number of
voxels is created. This user-defined number represents the maximum thickness between the
luminal domain and the tissue domain bounds. This new mask will serve to automatically limit
the resulting voxels to those already existent within the arterial tissue mask. The arterial tissue
mask was thus discarded as the controlled dilation mask replaced it. This process was attempted
repeatedly to find an ideal dilation value which minimized the inclusion of excess voxels.
Connective pixels were then identified between regions of interest and virtually severed. This
process was done with manual editing tools in both 2-D and 3-D, which are time-intensive userbased subjective operations. For the arterial tissue this process was much more prone to operator
bias. CT scans do not offer excellent soft tissue differentiation, and often boundaries were
approximated as best as possible. A final region growing operation successfully removed any
remaining disconnected voxels.
2.2.3 3-D Modeling

Completed 2-D masks represent voxel domains, where every voxel is a 3-D rectangle.
The direct 3-D interpretation of a mask would therefore appear as the sum of many connected
rectangular volumes, which would feature a very rough and noisy surface. The masks are
modeled in 3-D through pre-defined and/or custom methods to produce surfaces in Standard
Tessellation Language (STL) format instead, which describes an unstructured surface by the unit
normal and vertices of triangles in 3-D Cartesian coordinates. The STL quality is still dependent

34

upon the voxel resolution, however, which differs between every scan. A generated STL model
highlighting the surface structure is shown by Figure 10.

Figure 10. A 3-D AAA lumen STL reconstructed from scan B2. STL format describes a
triangularly meshed surface, as can be clearly seen. The data exported from the 3-D models for
Chapters 3 & 4 uses point cloud format. This format will only contain the vertex information
from STL models.

A wide variety of options exist for generating and refining an STL model from a mask,
thus STL generation from masks was standardized according to the a set of parameters (Table 3).
Given that the voxel resolution and therefore the STL resolution is not standardized no
standardized smoothing or wrapping criteria could be clearly established past the initial STL
model generation. Laplacian smoothing was applied via a subjective number of iterations from 0
- 10 to each STL at a constant smoothing factor of 0.7 (Table 4). The goal of smoothing is to
reduce surface irregularities, spikes, and noise with qualitative visual inspection as the means to
gage the reduction.

35

Table 3. Optimal STL model generation parameters used for every STL model generated in the
study. Contour interpolation results in a smoother model with less gaps and is recommended for
medical CT applications [59]. Matrix Reduction is set to 1, equivalent to no matrix reduction,
which inherently nullifies the need to select a matrix reduction algorithm as well (not shown).
Interpolation
Smoothing

Triangle Reduction

Contour instead of Gray Value
Iterations: 2
Smooth factor: 0.3
Reducing mode: Advanced Edge
Tolerance: Automatic (mm)
Edge angle: 10 degrees
Iterations: 3

Table 4. The number of employed smoothing iterations at a constant smoothing factor of 0.7 on
each STL model after optimal model generation. The smoothing was limited to a maximum of
10 iterations to avoid significant loss of detail.

Scan
A1
A2
B1
B2
B3
C1
C2
D1
D2
D3
F2
F3
F4
F5
F6
F7

Smooth Iterations
Tissue
Lumen
10
6
8
5
10
8
10
6
10
5
8
4
9
5
3
7
10
4
7
7
9
10
10

Scan
G1
G2
G3
G4
G5
G6
H1
H2
H3
H4
H5
H6
H7

Smooth Iterations
Tissue
Lumen
3
3
3
3
8
3
6
4
4
8
6
6
8
6

36

Scan
I1
I2
I3
I4
I5
I6
J1
J2
J3
J4
J5
K1
K2
K3
K4
K5

Smooth Iterations
Tissue
Lumen
4
4
4
5
5
5
8
3
8
4
10
4
10
4
10
4
10
4
4
5
5
4

For each scan surface models for the arterial lumen, arterial tissue, and adjacent vertebrae
are generated, then smoothed as judged necessary. The 3-D STL models can then be registered.
Registration refers to translating, rotating, and scaling a data set to a common position with a
reference data set.

Within the software three relevant modes of registration exist.

Image

Registration requires landmark point pairs to be identified in a reference image and secondary
image. The landmark points can be identified in the sagittal plane, coronal plane, and transverse
plane. The 4×4 transformation matrix that best relates the sets of landmark points is then
automatically calculated and applied to the secondary image.

Point Registration requires

landmark point pairs to be identified in a reference position and desired position. The landmark
points can be identified in the sagittal plane, coronal plane, transverse plane, and in 3-D. The
4×4 transformation matrix that best fits the set of landmark points is then automatically
calculated and applied to the selected STL.

Global Registration simply requires the

identification of a fixed reference STL and a moveable STL. The distance field between the
STLs is then automatically minimized over the course of a user-defined number of iterations.
Each iteration of global registration calculates a 4×4 transformation matrix that further reduces
the distance field and automatically applies it.

Global registration is chosen for its high

accuracy, ease of use, repeatability, and objectivity.

Global registration is an iterative process based on minimizing the distance field between
STLs. The sum of all of a patient’s scans are first imported to a common Mimics workspace, for
example scans A1 & A2 of patient A. The existing STLs for a given scan are then linked such
that transformation of one model would affect the linked models identically, but minimization of
the distance field is unaffected by the linked models. In essence this means the distance field is

37

minimized between two models, and the resulting transformation matrix is applied to any
additional desired models. A fixed reference vertebral STL and a second moveable vertebral
STL are chosen then the moveable vertebral STL is manually translated and rotated to be in as
identical a position as possible to the reference vertebral STL by eye. The two vertebral models
are selected and global registration is applied for 1000 iterations to register the moveable set of
STLS according to the transformation matrix describing the minimum distance field. This
process was very quick (<1 minute) and was repeated for all STLs from a patient on a common
reference vertebral STL.

The resulting STLs of a patient include the vertebral STL from each scan aligned with
vertebral STL of the first scan, and all associated lumen and tissue models positioned
automatically based on this alignment.

This registration method allows for the unique

visualization of the true spatial differences of the AAA geometry and offers insight into the
surface evolution of an AAA.
2.3 Results

Lumen models were generated for 45 scans, and tissue models were generated for 14
scans. Figures 11 and 12 represent lumen STL models from each patient from the front coronal
plane and right sagittal plane views, respectively. These models include the renal arteries, iliac
arteries, sometimes the mesenteric arteries, and occasionally the celiac trunk. The region of
interest for Chapter 3 onwards is the infrarenal aorta, and only the infrarenal region is generated
for all 45 scans as shown by Figure 13.

38

Figure 11. Front coronal plane view of 3-D luminal STL models generated for all 11 patients.
From top-left to bottom-right in reading order the STLs originate from scans: A2, B3, C2, D3,
E2, F7, G6, H7, I6, J6, K5.

39

Figure 12. Right sagittal plane view of 3-D luminal STL models generated for all 11 patients.
From top-left to bottom-right in reading order the STLs originate from scans: A2, B3, C2, D3,
E2, F7, G6, H7, I6, J6, K5.

40

Figure 13. A 3-D STL model of a reconstructed infrarenal AAA luminal surface as seen from the
left sagittal plane (left), front coronal plane (center), and right sagittal plane (right). A bounding
box indicates the infrarenal region specifically to be exported to a technical computing software
package.

An example of two registered vertebral models from scans I1 & I3 is shown (Figure 14).
It is clear that the models positions match extremely closely given that they are independently
generated from different CT scans.

41

Figure 14. An example of two registered vertebral models from scans I1 & I3 as viewed from the
left (left), front (center), and right (right). Each model is the same color, and has been registered
to match the position of its pair. Near the top and bottom of the models clipping can sometimes
be observed, which highlights the boundaries of a model (red). In general, however, it is clear
that the models’ positions match extremely closely given that they are independently generated
from different CT scans. The associated lumen and tissue models are also automatically
registered (not shown) according to the transformation matrix generated by minimizing the
distance field between vertebral models.

42

Tissue models are very time consuming to produce, and as such they are only produced in
the infrarenal aorta with the absence of any branches. Tissue models were generated for 14 of 45
scans based first on the order the data was received and second on the number of scans available
per patient. Eventually all of the data will be processed. Tissue models underwent the same
registration process concurrently as the lumen as shown by Figure 15. Additionally the arterial
wall thickness can be visualized through the transparent overlay of a respective lumen and tissue
model as shown by Figure 16.

Figure 15. Transparent and registered front coronal plane view (left) and right sagittal plane view
(second from left) of 3-D luminal STL models for B1 and B3. Transparent and registered front
coronal plane view (second from right) and right sagittal plane view (right) of 3-D tissue STL
models for B1 and B3. The lumen and tissue models presented are not to scale, as the lumen
models are inherently always smaller than the tissue models.

43

Figure 16. Transparent and registered front coronal plane view (left) and right sagittal plane view
(right) of 3-D luminal and tissue STL models for B3.

Regions of significant vertebral interaction with the vessels were identified for some
cases. The region with significant vertebral interaction displays flattening of vessel curvature in
the cross-sectional view (Figures 17 and 18). The vertebral surface can feature ossifications
originating and projecting from the spine known as osteophytes. The presence of osteophytes on
the spine increases the complexity of the relationship between an AAA and the spine, but not the
nature of the relationship; osteophytes serve as a growth barrier and have an influence on vessel
flattening (Figures 17 and 18). The prevalence and variability of osteophyte growth on the spine
is patient-specific and must be considered on a case-by-case basis.

44

Figure 17. A transverse slice un-zoomed (upper-left) and zoomed (upper-right) from scan A2. A
3-D reconstruction of the AAA tissue geometry for scan A2 from the right (bottom-left) and an
opened version of the AAA tissue geometry for scan A2 from the right (bottom-right). Red
circles indicate linked areas between images. Within the red circles it is clear that the AAA
tissue is interacting with the spine and that the local geometry of the tissue is affected.

45

Figure 18. A transverse slice un-zoomed (upper-left) and zoomed (upper-right) from scan A2. A
3-D reconstruction of the AAA tissue geometry for scan A2 from the left (bottom-left) and an
opened version of the AAA tissue geometry for scan A2 from the left (bottom-right). Red circles
indicate linked areas between images. Within the red circles it is clear that the AAA tissue is
interacting with the spine and that the local geometry of the tissue is affected.

The inferior vena cava was modeled for scan B2 to visualize its interactions with an AAA
(Figure 19). The shape of the aorta deviated from standard anatomical positions below the left
renal vein and we speculate that the left renal vein can serve as an anchor to the infrarenal during
46

the aneurysm progression. It may also be possible that the physical constraint of the left renal
vein over the vessel and the tethering of the renal arteries provide a strong confinement, or an
anchor at the region of the aorta, below which the aorta gradually bends as the lesion progresses.

Figure 19. Reconstructed 3-D models for the vertebrae (white), AAA lumen (red), and inferior
vena cava tissue/lumen (blue) as seen from the left (left), front (center), and right (right) for scan
B2.
2.4 Discussion

The 3D models generated from the biomedical software are relatively good, but there is
some degree of error inherent to segmentation and modeling. The resolution is approximately 0.7
mm per pixel, which is high for CT images, but some detail can still be lost in areas of thin
tissue, such as arterial tissue compressed against the spine (Table A.1). CT scans are also not
instantaneously obtained for the whole region. For this collection the entire torso was scanned,
usually within 10 seconds. The resulting CT scan is therefore taken to be an average of the
geometry in vivo, which is important to consider when observing the aorta given it is constantly
undergoing cyclic deformation by means of the cardiac cycle.
47

It is shown from the results of Figures 11 and 12 that the reconstructed models represent
AAA geometries accurately. Although smoothing was applied subjectively, the surface finish of
each model appears to have similar quality. Showing all results for 45 lumen models, 45
vertebral models, and 14 tissue models is not practical and, hence, a subset of images was chosen
per patient to be representative. As the infrarenal aorta begins at the renal branches it is clear for
each case that the vessel is at some angle from the z-axis. Patients C and J do not support this
trend as strongly, but upon examination it is clear that all vessels present some degree of
tortuosity. This tortuosity is three-dimensional, as the same trend can be observed in the coronal
and sagittal planes. A significant presence of an ILT is also likely to be a factor affecting the
tortuosity.

The majority of AAAs contain an ILT layer, though the thickness and circumferential
position of the ILT varies significantly between patients. For example patient J presents a
relatively small predominantly anterior ILT while patient D presents a relatively large
predominantly left lateral ILT. There is no uncontested cause and effect relationship between
patient specific AAA geometries / hemodynamics and ILT formation sites. Regardless of the
factors influencing ILT growth sites, the presence or absence of an ILT affects the biomechanical
relationship between the AAA and surrounding tissues; the spine serves as a rigid growth barrier
to AAA expansion, and if ILT is present in the circumferential region in contact with the spine
then disease progression is affected [44–46]. Arterial tissue and ILT carry different mechanical
properties, which directly leads to different interface mechanics when considering vertebral /
vessel interaction. It is observed from Figures 17 and 18 that extended areas of contact between
tissue and vertebrae present marked vessel flattening within the cross section. This vessel

48

flattening was observed both in the case of vessel contact with normal vertebral tissue and vessel
contact with vertebral ossifications. The inferior vena cava is adjacent to the abdominal aorta
and may interact with an AAA.

For the lesion shown in Figure 19, however, no clear

morphological features representing interaction between an AAA and the inferior vena cava are
observed, although the idea of the left renal vein serving to anchor an AAA in place is clear
when visualized.

From these observations it is hypothesized that the observed change of the centerline,
eccentricity and asymmetry of the vessel, and the curvature of the vessel cross section polylines
are morphological features representing the effects of surrounding tissue.

49

CHAPTER 3. Morphological Features Influenced by Surrounding Tissue
3.1 Introduction

The pathology of AAA is such that enlargement of the aorta is associated with disease
progression. The simplest measure of aortic enlargement is maximum diameter. Aneurysms less
than 5 cm to 5.5 cm in diameter are considered small and aneurysms over 5 cm to 5.5 cm are
considered large. Surgical treatment is usually only recommended for large aneurysms unless
patient-specific concerns warrant otherwise. As covered in Chapter 2, US is the usual imaging
modality used for screening and monitoring of small AAAs, while CT is used when intervention
is deemed necessary. Regardless of imaging modality, values of the maximum diameter can be
extracted according to four criteria: plane of acquisition, axis of measurement, caliper placement,
and diameter selection [51]. Critical reviews of the clinical obtainment of AAA maximum
diameters are offered by Long et al., Ferket et al., and Powell et al. from which much of this
review is derived [51–53]. The plane of acquisition, axis of measurement, and caliper placement
variability between studies is described.

The plane of acquisition refers to whether the maximum diameter is obtained from the
anatomical (absolute) reference planes or the aortic (relative) planes. The anatomical reference
planes are the sagittal, coronal, and transverse planes. The aortic reference planes are normal to
the arterial centerline, such that the (x, y, z) coordinate basis differs as the centerline is traversed,
especially for tortuous vessels; a plane normal to the centerline is known as a longitudinal plane.
US can be aligned with either reference by a skilled operator, while CT scanning is automatically

50

aligned to the transverse plane as the native plane. CT MPR usually generates the additional
anatomical planes, the sagittal and coronal planes, but other MPR algorithms exist to generate
planes for the aortic reference. Non-standard MPR algorithms for CT scans inherently suffer
from a loss of resolution however, as the scanning algorithm cannot be adjusted from the
anatomical standard planes to match unique MPR algorithms. MRI scanning has the unique
advantage of being able to take measurements in any plane, and additional MPR algorithms exist
to take advantage of this.

The axis of measurement refers to the chosen axis within to the cross-sectional plane of
measurement for maximum AAA diameter. The anteroposterior and transverse axes refer to the
sagittal and transverse axes respectively.

The transverse cross-sectional view allows for

simultaneous multi-axis measurement while the sagittal and coronal views only allow for a
single axis of measurement. For US measurements are performed during scanning, while CT
MPR offers delayed analysis where any axes can be used for measurement.

After identifying the acquisition plane and axis of measurement the caliper positioning
must be defined. The caliper positioning refers to which surface of the AAA wall to bound the
maximum diameter measurement to whether with calipers for an US or digital measurements for
a CT reconstruction. The positioning can be internal, external, leading edge to leading edge, and
outer anterior to inner posterior. A total of 56 studies were critically reviewed by Long et al.,
Ferket et al., and Powell et al. to identify the choices of plane of acquisition, axis of
measurement, and positioning of the calipers defined by different clinical groups [51–53]. The
results of their studies are summarized in Table 5.

51

Table 5. Summary of the definition of techniques employed over the survey of 56 studies by
approximate percentage [51–53]. Some groups reported multiple methods of measurement to be
in use rather than a single method, resulting in total percentages not equal to 100%.
Plane of Acquisition
Not Specified 51%
Longitudinal 34%
Transverse 30%
Axis of Measurement
Anteroposterior 63%
Not Specified 37%
Transverse 30%
Positioning of the Calipers
Not Specified 61%
External 27%
Internal 4%
Leading edge to leading edge 2%
Outer anterior to inner posterior 4%
These results clearly show that no undisputed standard methodology can be extracted
from the literature to determine maximum AAA diameter.

An axis of measurement

(anteroposterior) was the only variable that was more often specified than not.
3.2 Method

From reviewing the literature it is clear that no standardized method of diameter
measurement exists. In general 4 steps are identified and followed, however. The plane of
acquisition is specified, the measurement axis is specified, the caliper positioning is specified,
and the final maximum diameter is obtained as either a maximum obtained value or a multi-trial
average value. For obtaining maximum diameters from the CT standard MPR of anatomical
planes the transverse plane was selected as the plane of acquisition, as 30% of reviewed studies
suggest. The longitudinal plane which 34% of studies suggest is subjectively defined or requires

52

centerline generation to perform relative MPR and attain true longitudinal planes to the
centerline. An example transverse measurement plane is shown by Figure 20.

Figure 20. A sagittal slice (left) of scan J5 and a transverse slice (right) of scan J5 with the
segmented lumen/tissue (red/green). The blue line indicates an anteroposterior maximum AAA
tissue diameter measurement axis in the transverse plane (left, right), the solid red line indicates
an anteroposterior maximum AAA tissue diameter measurement in the longitudinal plane (left),
the dashed red line indicates a crude longitudinal axis approximation (left), and the orange line is
simply the AAA tissue diameter measurement at 90° in the transverse plane to the
anteroposterior measurement axis. The caliper positioning for this case was external/tissue.

The measurement axis specified is anteroposterior unless it is clear that the vessel is
noticeably asymmetric in which case multiple axes were tested and the axis corresponding to the
maximal measurement was chosen. The caliper positioning for obtaining tissue values was
external while the caliper positioning for obtaining lumen values was internal. Specifying the
axis of measurement for the lumen additionally considers the presence of ILT and its asymmetry.
Measuring both the maximum diameter of the lumen/tissue and the diameter on the axis 90° to it
within the plane of acquisition allows for a simplistic calculation of eccentricity based on the
53

diameter ratio. Asymmetric ILT can lead to very different axes of measurement for the
maximum AAA lumen diameter as shown by Figure 21.

Figure 21. A transverse slice of scan D2 with the segmented lumen (red), maximum
internal/lumen diameter axis (blue) and axis 90° to the maximum internal/lumen diameter (red).
It is clear that anteroposterior measurements will not yield maximal results in the presence of an
eccentric ILT. This is the reason that the maximum diameter axis was used and anteroposterior
was not strictly used in all cases.
3.3 Results

Maximum diameter measurements for the lumen and tissue were recorded from
anatomical plane CT MPRs for each patient. For patient C two aneurysmal sacs are recorded,
the primary is referred to as Patient C, while the secondary is referred to as Patient C (sac)
(Figure 22).

54

Patient C (sac)

 Patient C

Figure 22. 3D CT reconstruction of patient C as viewed from the right. Two aneurysmal sacs are
present and the primary is referred to as Patient C, while the secondary is referred to as Patient C
(sac).

The values obtained from measuring the AAA maximum tissue diameter in the transverse
plane preferring an anteroposterior measurement axis unless asymmetry is present, with external
caliper positioning, and no value averaging are plotted by Figure 23 (Table A.3). The linear
approximations of annual expansion rate are calculated by fitting the data linearly and the results
are shown by Table 6.

55

Maximum Transverse Diameter Growth of Arterial Tissue
75
A

65

B

60

C

55

C (sac)

50

D

45

F

40

G

35

H

30

Max Diameter (mm)

70

I
J

25
0

500

1000

1500
2000
Days

2500

3000

3500

K

Figure 23. The maximum transverse diameter of the arterial tissue vs. scanning date for each
patient.

Table 6. Linear approximation of annual transverse maximum AAA tissue diameter growth and
2
R fit value.
Patient
A
B
C
C (sac)
D
F
G
H
I
J
K

Dmax linear increase (cm / year)
0.3010
0.3931
0.2462
0.1208
0.5410
0.4796
0.4670
0.2427
0.1888
0.2542
0.4064

56

2

R
1.0000
0.9569
1.0000
1.0000
1.0000
0.9804
0.9542
0.9536
0.9253
0.9381
0.9406

The annual expansion rates of the tissue were additionally calculated for successive scans
and plotted as shown by Figure 24. A measure of eccentricity of the tissue was taken as shown
by Equation 1 and plotted for each scan as shown by Figure 25.
( )

Annual Expansion Rate of Maximum AAA Tissue Diameter Growth
Calculated on a Successive Scan Basis

Annual Growth (cm / year)

1.6
1.4

A

1.2

B

1.0

C

0.8

C (sac)
D

0.6

F

0.4

G

0.2

H

0.0
Scans 6-7

Scans 5-6

Scans 4-5

Scans 3-4

Scans 2-3

I
Scans 1-2

-0.2

J
K

Figure 24. The annual expansion rate of maximum AAA tissue diameter growth calculated on a
successive scan basis for each patient.

57

Eccentricity of AAA Tissue on a Successive Scan Basis
1.20
A

1.16

B

1.14

Eccentricity

1.18

C

1.12

C (sac)

1.10

D

1.08

F
G

1.06

H

1.04

I

1.02
Scan 7

Scan 6

Scan 5

Scan 4

Scan 3

Scan 2

Scan 1

1.00

J
K

Figure 25. The tissue eccentricity based on eccentricity = (Dmax) / (D @ 90° to Dmax) for each
patient scan.

The values obtained from measuring the AAA maximum lumen diameter in the
transverse plane preferring an anteroposterior measurement axis unless asymmetry is present,
with internal caliper positioning, and no value averaging are plotted by Figure 26 (Table A.4).
The linear approximations of annual expansion rate are calculated by fitting the data linearly and
the results are shown by Table 7.

58

Maximum Transverse Diameter Growth of Arterial Lumen
65
A

Max Diameter (mm)

60

B

55

C
50

C (sac)
D

45

F

40

G
35

H

30

I
J

25
0

500

1000

1500
2000
Days

2500

3000

3500

K

Figure 26. The maximum transverse diameter of the arterial lumen vs. scanning date for each
patient.

Table 7. Linear approximation of annual transverse maximum AAA lumen diameter growth and
2
R fit value.
Patient
A
B
C
C (sac)
D
F
G
H
I
J
K

Dmax linear increase (cm / year)
0.0151
0.0183
0.0059
0.0065
0.0141
0.0027
0.0140
0.0006
0.0041
0.0038
0.0030

59

R2
1.0000
0.9705
1.0000
1.0000
0.9879
0.3166
0.9694
0.0542
0.9954
0.7776
0.6189

The annual expansion rates of the lumen were additionally calculated for successive
scans and plotted as shown by Figure 27. A measure of eccentricity of the lumen was taken as
shown by Equation 1 and plotted for each scan as shown by Figure 28.

Annual Expansion Rate of Maximum AAA Lumen Diameter Growth
Calculated on a Successive Scan Basis
1.2

Annual Growth (cm / year)

1.0

A

0.8

B

0.6

C

0.4

C (sac)

0.2

D

0.0

F

-0.2

G

-0.4

H

-0.6
Scans 6-7

Scans 5-6

Scans 4-5

Scans 3-4

Scans 2-3

I
Scans 1-2

-0.8

J
K

Figure 27. The annual expansion rate of maximum AAA lumen diameter growth calculated on a
successive scan basis for each patient.

60

Eccentricity of AAA Lumen on a Successive Scan Basis
1.6
A

Eccentricity

1.5

B
C

1.4

C (sac)
D

1.3

F
G

1.2

H
1.1

I
K

Scan 7

Scan 6

Scan 5

Scan 4

Scan 3

Scan 1

1.0

Scan 2

J

Figure 28. The lumen eccentricity based on eccentricity = (Dmax) / (D @ 90° to Dmax) for each
patient scan.
3.4 Discussion

The results of the maximum tissue diameter study from Figure 23 are straightforward in
that they all appear relatively linearly increasing. Other methods of fit were tested and found to
2

have lower R coefficients when compared to linear fits. Some positive trend was expected as
AAA diameters continually increase during disease progression. A decrease in AAA maximum
transverse diameter of the arterial tissue is only observed in the case from I1 – I2. The value
appears to be approximately 3mm less than expected from the trend. Most likely the plane of
2

acquisition was improperly judged and this value is not representative. The R fitting values
from Table 6 obtained by fitting linear trend lines to the tissue diameters are excellent at an
average of 0.9681 ± 0.0287. The actual rates of increase mostly all fall below 0.5 cm/year which
61

indicate relative stability of the AAAs. The estimated rate of increase for patient D, however,
was over 0.54 cm / year. Some studies suggest less than 1 cm / year is stable, others suggest less
than 0.5 cm / year is stable. Regardless, it is unclear why this patient’s AAA, with a higher
growth rate and the fact that monitoring began at a maximum tissue diameter over 5 cm, was
allowed to progress as far as it did without treatment. On that same note, inspection of the
maximum tissue diameter values clearly shows at least 21 / 49 measurements over 5 cm in
maximum tissue diameter. Such measurements would suggest the dataset equally represents
small and large AAAs, even though all are pre-EVAR. The successive scan calculated growth
rate did not yield any interesting trends, aside from the possibility of some measurement error
between scans F2 – F3 yielding an uncharacteristically high number. The scanning algorithms
and resultant voxel resolutions are quite different between scans, leading to different standard
deviations of measurement (Tables A.1 and A.2). For smaller vessel diameters the resolution of
the voxels influences the results more significantly, such as between scans F2 – F3.

By

measuring the basic eccentricity for each scan it was hoped a clear trend would emerge
reinforcing the idea of asymmetric growth, but no clear trends emerged across patients.

The results of the maximum lumen diameter from Figure 26 showed increasing diameters
over time, as expected. For patients G & F it is interesting to note the large dip in diameter
between 2000 – 2500 days. This dip is most likely due to growth of an ILT reducing the lumen
volume at a rate faster than the entire vessel can expand. Fitting linear trend lines to the
2

maximum diameter growth resulted in R fitting parameters averaging 0.7901 ± 0.3276 (Table
7). Therefore the maximum lumen diameters did not fit linear trends as well as the maximum
tissue diameters did. The successive scan calculated expansion rates offered both positive and
62

negative expansion rates. No doubt this is a combination of the chosen acquisition plane and
relative growth of an ILT.

Measuring maximum lumen diameter changes may be more

beneficial if ILT measurements are made in the same plane concurrently. This would offer
insight to the influence of ILT on maximum lumen diameter measurements. The eccentricity
showed some favorable results and the range of eccentricity values recorded suggest significant
variation amongst patients. Lumen eccentricity outweighed tissue eccentricity in that at least 13
lumen eccentricities were greater than the maximum tissue eccentricity. If the maximum tissue
diameter and maximum lumen diameter were always measured in the same plane the ILT would
represent the difference between the values after accounting for wall thickness. This would then
no longer represent the maximum lumen diameter, however, as the maximum lumen diameter
and maximum tissue diameter often occur in different planes and axes.

The methods followed suffer from subjectivity. Experienced screening personnel have an
in-depth intuitive understanding of how to make these measurements according to their
institution’s guidelines to achieve reliable and repeatable results. No such personnel aided in the
attainment of these measurements, and the overall accuracy was most likely decreased. No
information could be found on established protocol for the analysis of longitudinal data. Perhaps
there exist some unknown conventions to ensure a logical progression of recorded values that
trained personnel are aware of.

As longitudinal studies become more widespread proper

protocols need to be established to standardize the interpretation of the results.

63

CHAPTER 4. Parameterization of Luminal and AAA Surfaces
4.1 Introduction

For an AAA in vivo there are multiple tissues contacting its boundary, none of which
have been previously, fully considered for their effect throughout disease progression.
Morphological features such as maximum diameter, asymmetry, cross-sectional polyline
curvature flattening, tortuosity, and eccentricity are significantly influenced by both surrounding
tissue and hemodynamic factors. In order to quantify either the combined or separate influence
of such factors during disease progression a precise characterization of aneurysm’s geometrical
evolution is needed. Multiple methods for geometrical parameterization of abdominal aortic
aneurysms (AAAs) have been previously developed using individual patient CT scan data but
the focus has been mainly on the association of such geometrical parameters with the rupture risk
and the efficacy of the parameterization is not fully investigated for a longitudinal study yet [60].
For this study a series of 3D models for AAAs in longitudinal studies, an arterial centerline
generation algorithm (Appendices B.1 – B.6), and a geometric parameterization procedure for
the arterial surfaces are developed.

An iterative algorithm was developed to generate a centerline for an arbitrary arterial
surface model by collecting the center points of maximally inscribed spheres that fit within the
surface boundaries at fixed travel intervals. The number of inscribed spheres is dependent on the
choice of travel interval. Using the sphere center points and 4th order polynomial base functions,
a smooth line approximation is made. The surface data is then parameterized to the centerline by
defining an (s, θ) coordinate system where s is the travel length in mm along the centerline and 0

64

≤ θ ≤ 2π radians. We then calculate r for each surface point where r is the minimum distance
from the centerline to the point on the surface at a given (s, θ). The parameterized surface r(s, θ)
is then approximated by normal shape functions with fitting coefficients using a nonlinear least
squares method. The change of the surface over time is calculated by the difference between
these surface functions. The approximated surface can be regenerated in Cartesian coordinates.
Alongside the parameterization technique a set of code was developed to characterize
morphological features of the surface. The maximum diameter at each plane, the simplified
eccentricity, and polyline curvature were then determined.
4.2 Methods
4.2.1 Importing the data

A 3-D point-cloud model is imported to Matlab. The point cloud model is essentially a
subset of an STL model developed in Chapter 2 in that point cloud models are comprised of
solely the vertices of STL models. The imported data is a matrix containing 3 columns for the x,
y, and z coordinates across a number of rows representing the total number of coordinates. The
model represents a closed surface with capped vessel ends. When defining the coordinate
system inherited from the biomedical image processing methodology care must be taken.
Depending on the scan, significantly different origin points were observed with respect to the
STL positions. Registration serves to identify relative positions of one data set relative to
another data set, but not across all sets. There is no single coordinate system to represent all
scans, but a set of planes parallel to the x-y, y-z, and x-z planes for every scan can be identified
(Figure 29). The unit of measure is millimeter.

65

x-y : Transverse
y-z : Sagittal
x-z : Coronal

Figure 29. The coordinate planes inherited from the anatomical planes where the coordinate axis
is such that the x-y plane is parallel to the transverse plane (green), the y-z is parallel to the
sagittal plane (red), and the x-z plane is parallel to coronal plane (blue).

The coordinates of M unordered points in the model can be defined by Equation 2.
(

)

,

-

( )

One requirement for the algorithm is that the vessel ends should be an open surface to
accurately approximate centerline points at the model bounds. Therefore the volume of points is
first truncated transversely immediately near the top and bottom to open the vessel ends with
minimal removal of other points. Two truncation planes are identified from the top: ( max(zN)2.5) mm, and from the bottom: (min(zN)+ 2.5) mm, respectively. Only coordinates of qN within
these bounds remain, and are processed through the rest of the calculations. Thus it is most

66

efficient to now consider that M be redefined as the total remaining points after truncation. The
only variables carried forward are qN for N = [1, …, M].

4.2.2 Centerline Generation Algorithm

The centerline generation algorithm serves to generate a centerline in the threedimensional space. The previous method used most for centerline generation is the one that
simply collects the centroids of transverse polylines [13], [23], [25], [29], [47–50].

The

algorithm proposed in this work generates a centerline by fitting maximally inscribed spheres
within the surface model constrained to planes approximately normal to the centerline at
relatively evenly spaced intervals (Figure 30).

Figure 30. The arterial surface (blue) with two maximally inscribed spheres (red) and their
respective constraint planes (green) shown. The centerline (black) is shown within the lumen
model for scan B2 as seen from an isometric point of view (left), parallel to the x-z or coronal
plane (center), and the y-z or sagittal plane (right).

67

The result is a set of n points pi = (xi, yi, zi) for i = [1, …, n] that is constructing the
centerline. The algorithm begins by selecting a set of points (q1, …, qK) whose z values are
bounded by min(qN(z)) and (min(qN(Z))+g), where g is a constant variable defining a small gap.
The first point p1 is then located to the middle of the set of points (q1, …, qK) in the x-y plane
and at the bottom in the z axis as shown by Equation 3.
( )

( )

( )

( )

( )

( )

.

[

( )/

( )
]

The distances dN for N = [1, …, M] to all points on the surface model qN for to the center point
p1 are calculated by Equation 4.
√(

( )

( ))

(

( )

( ))

(

( )

( ))

( )

The resulting vector A corresponding to ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ at N for min(dN) is then recorded. The center
point p1 is then translated in the opposite direction of the projection of A onto the x-y plane C,
simply the x and y components of A in this case, a fixed distance t where t is a user-defined
variable, usually 1 – 2 mm. If the absolute position of the center point has not changed by 1 mm
over 15 iterations then the algorithm halves the translation amount and continues. This process
of reducing the translation amount is repeated until a maximum number of iterations is met.
When this process is completed the point p1 represents the approximate center point of a
maximally inscribed sphere at the bottommost x-y plane of the volume. Essentially the z-

68

coordinate of a sphere is locked, and the sphere is allowed to grow and translate in the x-y plane
until it reaches a maximum size and the center point p1 is recorded. For the second center point
p2 the final (x,y) coordinates of p1 are used, but vertical translation kv is applied where v is a
user-input. Otherwise the estimation of p2 continues the same as p1.

For each successive point referring to Figure 31 is extremely useful to visualize the
terminology used. For each successive point for n = 3 onward the normal of the vector between
the previous two center points B = pn-1 – pn-2 is calculated and serves as the translation direction
and v serves as the translation magnitude to determine the initial distance between the previous
center point pn-1 and the current center point pn. This serves to better place the initial guess to
more quickly and accurately converge on the best value of pn. Using B to place the initial guess
of pn is especially important for tortuous vessels. The point on the surface of qN at a minimum
distance to the center point pn is calculated as well as the vector from pn to qN as A such that the
minimum distance is equivalent to ||A||. The center point pn is then translated in the opposite
direction of C, the projection of A onto the plane normal to B, by 1 mm where C is defined by
Equation 5.
(

‖ ‖
‖ ‖

69

)

( )

If the absolute position of the center point has not changed by 1 mm over 15 iterations
then the algorithm halves the translation amount and continues. This process of reducing the
translation amount is repeated until a maximum number of iterations is met.

Figure 31. The lumen point cloud model for scan B2 colored by z-coordinate values with a red
sphere at approximately (155, -180, -163) (left), and the sphere itself colored by z-coordinate
values with vectors indicating the orientations of vectors A (blue), B (black), and C (red) for a
sample center point (right).

This process is repeated until the z-value of an initial guess for pn exceeds the maximum
z-value of the model bounds, at which point the z-coordinate of pn is reset to the maximum zvalue and the final center point is recalculated using x-y plane projections in place of projections
on the plane normal to B. The result is a collection of n points defined by pi = (xi, yi, zi) for i =
[1, n] in Cartesian coordinates which approximate points on the centerline of the vessel, all other
variables are discarded moving forward. The points attained from the centerline generation

70

algorithm represent approximations subject to surface roughness, approximation functions are
now applied to attain a smooth centerline function and derive points on that line. The only
variables carried forward to the next step are qN for N = [1, …, M] and pi for i = [1, …, n]. An
example of centerline approximation points is shown by Figure 32.

Figure 32. The final collection of centerline points plotted within the lumen point cloud model
for scan B2 colored by z-coordinate value as seen from an isometric point of view (left), parallel
to the x-z or coronal plane (center), and the y-z or sagittal plane (right).

71

4.2.3 Smooth Centerline Approximation

Based on the existing set of points from the centerline generation algorithm some new
terms can be introduced. The successive distance between each of n points defined by pi = (xi,
yi, zi) for i = [1, …, n] is defined by resultant vectors as shown by Equation 6. The total distance
traversed along the centerline from p1 to pi is defined as li by Equation 7.

√.

/

.

/

.

,

/

,

∑

-

-

( )

( )

These quantities are shown for clarification in a diagram for n = 7 by Figure 33.

p7
p6
p5

l7
l6

d5

l5

d4

p4
p3

d6

l4

d3

l3

d2
p2

l2
d1

l1

p1

Figure 33. The meaning of p, d, and l can be hard to ascertain from equations alone. The
variable d represents the point to point distances between p values. The variable l represents the
sum of d values thus far at a given l.

72

The sum of the point to point distances ln represents the length of the centerline. We
define an integer set of m integers existing from [0, ln]. Therefore if ln = 149.63 for example,
then m = 149 and 150 points will be defined on the smooth line approximation for sk where s ϵ
+

ℝ0 and k = [0, m].

In this case a discrete set of sk values is defined to produce and analyze a

discrete set of longitudinal planes.

From this definition of an s parameter a smooth

approximation function can be generated based on the set of points pi. This method is adapted
from a previous method employed in our lab by Zeinali-Davarani et al. [60]. Depending on a set
of conditions outlined in Table 8 an appropriate interpolation function element
constructed and its partial derivative

( )

calculated for each sk.

73

( ) is

Table 8. Conditional table for shape function assignment of smooth line approximation.
( )

Condition
(

) (

)

(

) (

)

Else

0
(

) (

)

(

) (

)

Else

0
(

) (

)

(

) (

)

Else

0
(

) (
(

)

) (

Else

)

0
(

) (
(

)

) (

Else

)

0

( ) are constructed and their partial

After the interpolation function elements
derivatives

( )

are calculated for each s the approximation function elements

derivatives

( )

are then normalized as

( )

( )

.(

( )

)(∑

( ) and

( )⁄

by Equations 8 and 9.

( )
∑
( ))/
(∑i
74

( )

( )
(.

( ) and

( )/ (∑

( ))

( )

))

( )

( ) is multiplied by pi and the sum of these multiplications at each

The interpolation function

s results in a point Xk = (xk, yk, zk) such that the total collection of points represent points that lie
on the smooth line approximation of the centerline as shown by Figure 34 by Equation 10.
( )

The derivative approximation function

∑

( )

( )

(

)

is also multiplied at each s by pi. The

normalized sum of these multiplications for a given s results in a vector N1 by Equation 11.

∑
‖

( )

(

)

‖

A unique vector N1 exists for each s and represents a centerline tangent vector at Xk =
(xk, yk, zk). From N1 at each s additional vectors N2 and N3 can be calculated to form unique
longitudinal coordinate bases at each point Xk = (xk, yk, zk). The only variables carried forward
to the next step are qN, Xk, s, and N1.

75

Figure 34. The final collection of centerline points plotted within the point-cloud lumen model of
scan B2 colored by z-coordinate value and the fitted centerline (black) as seen from an isometric
point of view (left), parallel to the x-z or coronal plane (center), and the y-z or sagittal plane
(right).
4.2.4 Longitudinal Acquisition Plane Generation

As described in the previous section, tangent vectors N1 exist for each s originating from
the point Xk for k = [1, m]. In order to fully define a longitudinal acquisition plane two
additional coordinate basis vectors N2 and N3 must be calculated (Figure 35). The transverse
plane, on the other hand, uses the known Cartesian standard basis {i, j, k}. The longitudinal
plane at each s is defined by projection of the x-axis onto the plane perpendicular to N2, and
therefore N3 = N1 × N2 by Equation 12.
(
(

‖

76

)
)‖

(

)

Figure 35. The vectors N1 (red), N2 (blue) and N3 (green) are plotted for each s originating from
the points Xk for k = [1, m] for the lumen model of scan B2.

Points lying on the longitudinal planes defined by (N1, N2, N3) can then be determined
through use of the dot product. For each s the vectors from Xk to all points on the surface model
qN, or ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , are calculated for all combinations for k = [1, m] and N = [1, M]. A dot product
inequality is evaluated and the vector results are stored for valid pairs of k & N as rh by
Equations 13 and 14.
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗

(

)

⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ( ) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ( ) ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ( )

(

)

77

h is the number of points Xk which lie in longitudinal planes. As rh is additionally populated the
associated surface points qN are stored in a new point-cloud volume vh representing only points
that lie on the longitudinal planes defined by (N1, N2, N3).

An example clarifying these

variables is shown by Figure 36.

Figure 36. The imported point-cloud lumen model composed of points qN for scan B2 colored by
z-coordinate value with points on the smooth approximated centerline shown (black, left), and
the subset of points vh colored by z-coordinate value from the point-cloud lumen model for scan
B2 which lie on longitudinal planes defined by (N1, N2, N3) (right). For example, the vectors in
rh from Xk to the associated points on the surface vh at the plane defined by (N1, N2, N3) for k =
40 are shown in green (right).

78

The determination of points lying on the longitudinal planes at s only using dot product
thresholding reports extra erroneous values when planes intersect the surface at more than one
boundary. To compensate for this phenomenon the magnitude of each vector from rh for a
longitudinal plane at s is calculated and analyzed. If a vector’s magnitude is greater than the
minimum vector magnitude in that plane by some threshold, say 5 mm, the vector is removed
from rh. A more statistically justified approach was originally implemented, but was found
unreliable across geometries. The method described requires a small degree of user-guided
iteration to ensure accurate results in place of full statistical automation, but was deemed an
appropriate trade-off with regards to development time (Appendices B.1 – B.6). The magnitudes
of vectors within rh represent the radius from the centerline points Xk to the surface points vh.
The mean and maximum radius at each longitudinal plane can be found as a useful measure of
growth (Figure 37).

79

Figure 37. The radii values are calculated for all h by ||rh||. The mean is then calculated at each
longitudinal plane for the lumen of scan B2 (left), and the circles of mean longitudinal radii at
each longitudinal plane can be visualized for the lumen of scan B2 (right).
4.2.5 Parameterization of Cross Sections into r(s, θ)

At each longitudinal plane defined by (N1, N2, N3) a collection of vectors rh exists in
Cartesian coordinates that describe the distance from the centerline point Xk to surface points qN.
To more efficiently analyze the data on a longitudinal plane basis a transformation to polar
coordinates takes place. In polar coordinates the magnitudes of rh vectors represent r. By
considering a suite of dot products between each rh vector, N1, and N2 within each longitudinal
plane at a given s the angular values θ within that plane associated with each rh can be obtained.
First a set of 4 calculations determines if the value of θ is 0, π/2, π, or 3π/2. If the dot product

80

does not correspond to one of those exact values then the region θ lies in can be determined and
θ subsequently calculated. This process is clarified by Table 9 and Figure 38.

Table 9. Conditional table for calculation of θ from a given rh vector. First it is determined
whether θ lies on the N2, N3, -N2, or -N3 axis. If it does not lie on an axis than the region it lies
in is determined and θ is calculated appropriately.

1

⁄

st

⁄
(

)
(

or
2

nd

(
(

)

)
)
(

or
(

)

81

)

Figure 38. Schematic representation for the calculation of θ from a given rh vector. First it is
determined whether θ lies on the N2, N3, -N2, or –N3 axis (left), and if it does not lie on an axis
then the region it lies in is determined and θ is calculated appropriately (right).
4.2.6 Coordinate transformation from r(s, θ) to (x, y, z)

Successfully parameterizing the data from (x, y, z) space to r(s, θ) space allows for unique
surface visualization of the cross-sectional radii along the vessel circumference relative to the
vessel centerline parameter s (Figure 39).

82

Figure 39. Parameterized surface data (red) from scan B2 shown as r(s, θ) surface.

From the parameterized coordinate system any point can be reconstructed within the
Cartesian coordinate system through coordinate transformations. Each longitudinal plane is first
assumed to share an origin with the Cartesian origin to determine the relative rotation of the
plane. Given that <N2, N3, N1> represents the unit vector basis at each s, and that <ei, ej, ek>
represents the unit vector basis in Cartesian coordinates a rotation transformation tensor between
each basis is attainable.

In reality there exist translation transformation vectors at each

longitudinal plane defining translation from the origin. These vectors are derived from the
centerline point Xk for each s. This two-step transformation allows for the reconstruction of the
surface from r(s, θ) coordinates to Cartesian coordinates at each s first by rotation and then by
translation. The transformation of any given r(s, θ) point at a given s to any given (x, y, z) point
is thus given by Equation 15.

83

[ ]

[

][

]

( )
[ ( )]
( )

(

)

It is apparent that reconstruction of the points into Cartesian coordinates is
straightforward and directly applicable to any set of data in longitudinal planes. Therefore any
fit or approximation of the parameterized surface is directly transformable back to Cartesian
coordinates (Figure 40).

Figure 40. Subset of the original point cloud data which lie on longitudinal planes (left) and the
data reconstructed from the r(s, θ) data (right).
4.2.7 Surface Approximation

After using the parameterization presented in the previous section, the distance of each
point on the surface from the centerline, r, is plotted with respect to two variables {s, θ} as seen
in Figure 39. In this section, we present a smooth numerical approximation of the surface with a
function r(s, θ).

The numerical code for the approximation in MatLab can be found at

84

Appendices B.7 and B.8. Figure 41 shows the schematic diagram for the parameterization and
the functional approximation.

1. Parameterize from
Cartesian to r(s, θ)

3. Fit surface with
shape functions

2. Reconstruct in
Cartesian coordinates

4. Reconstruct in
Cartesian coordinates

Figure 41. A sample arterial lumen surface from scan C2 (left) is shown to undergo the steps of
parameterization (center) followed by surface approximation (right) with the optional
reconstruction to Cartesian coordinates shown for clarity in both cases.

The surface is approximated by a linear combination of fitting parameters and shape
functions, i.e., (

)

∑

(

), where the shape functions are defined at Ns x (Nθ - 1)

nodal points. Nθ and Ns are the number of grids in the domain along the θ and s axes. The
parameterized data is special in the sense that θ spans from 0 to 2π radians, but the continuity
between θ = 0 radians and θ = 2π radians should be ensured. For this reason, the shape functions
at θ = 2π become not independent and the shape functions are defined from 1 to (Nθ - 1) in the
circumferential axis. Given the variable θ is connected at 0 and 2π, the shape function must

85

reflect this. For this, let us define new variables dθi and dsj by Equation 16. The shape function
th

at k node is then defined by Equation 17.
*|

|

|

|+

.

|

|

/

(

)

(

)

β1 and β2 are two fitting parameters, which requires optimization for the given number of shape
functions. Finally the resulting approximations

can be obtained by using the least squares

method. To optimize β1 and β2 an iterative solver, fminsearch [61], in MatLab is employed to
minimize the squared sum of the residual error to a convergence tolerance of 0.0001 based on
initial guesses of β1 = 0.009 and β2 = 1.5.

This allows for the determination of points

approximated by shape functions from each r(s, θ) value based on the residual error between the
data and the approximation value (Figures 42 and 43). From this optimization any point of the
parameterized surface can then be defined by the resulting shape functions (Figure 44).

86

Figure 42. Parameterized surface data (red) and the respective points from the best fit surface
(blue) from scan B2 shown as r(s, θ) surfaces.

Figure 43. Plot of the residual error between the parameterized surface data and the respective
points from the best fit surface from scan B2 shown as an r(s, θ) surface.

87

Figure 44. Parameterized surface data (red) and the respective cross-sectional lines from the best
fit surface (blue) from scan B2 shown as r(s, θ) surfaces.
4.2.8 Eccentricity and Polyline Curvature

Eccentricity and polyline curvature are two morphological features which can be
examined on the cross-sectional level. Mathematically the eccentricity of a conic-section is
based on the ratio between the semimajor axis and the distance along the semimajor axis at
which the focus lies. Some alternatives exist which are simplified and carry slightly different
formulations, though all represent the tendency of a conic-section towards or away from
circularity. A clear measure of eccentricity intrinsically understood by medical professionals is
simply the ratio between the diameter of a cross-section’s major and minor axes. This can be
automated by finding the largest distance between two points of a cross-section iteratively as the
major axis, and then identifying the axis joining the point pair closest to 90° to the major axis as
the minor axis (Figure 45).

88

(x2, y2)
(x3, y3)

(x1, y1)

(xc, yc)

Figure 45. Major axis (left, red) and minor axis (left, green) of a cross-section (blue) used for
eccentricity determination. Terminology of the 3-point polyline curvature method is detailed
(right).

To define polyline curvature, a two dimensionally derived value, many methods were
available. The chosen method was by defining the curvature at a given point (x2, y2) to be the
inverse of the radius of a circle fit to that point and two neighboring points (x1, y1) and (x3, y3)
by first finding the center of the circle (xc, yc) by Equations 18 and 19.
(

(

)

)

(

)

(

)

89

(

)

(

)

(

)

After identifying the center of the circle based on 3 points the radius, and therefore the
inverse of the radius, can be calculated by considering the center of the circle (xc, yc) and any of
the points (x1, y1), (x2, y2), or (x3, y3) by Equations 20 and 21.
√(

)

(

)

(

)

(

)

The processes for eccentricity and polyline curvature were both automated (Appendices
B.9 and B.10)
4.3 Results
4.3.1 Centerline generation

The centerline generation algorithm was ran for 59 cases and it was able to successfully
generate centerline points for 55 / 59 cases, or ≈ 93% of cases with the default input parameters.
The failed models were lumen models from scans D3, F5, H3, and K3 with their error detailed in
Appendix C (Figures C.1 – C.4). By tweaking the input parameters to generate fewer points at a
higher number of iterations per point models were then successfully generated for scans F5, H3,
and K3. Scan D3 featured a highly irregular anomaly that was outside the scope of the algorithm
completely. The centerlines of patient H are highlighted below (Figures 46 and 47).

90

Figure 46. Isometric view of the centerlines (left) and a transverse view (right) for scans H1 – H5
where black, blue, green, red, cyan, magenta represent H1, H2, H3, H4, H5, H6 respectively. In
this case the general region of the spine is the x-z plane at y = 100 such that the arrow indicates
migration away from the spine.

91

Figure 47. x-z/frontal/coronal view of the centerlines (left), and y-z/sagittal view of the
centerlines (right) for scans H1 – H5 where black, blue, green, red, cyan, magenta represent H1,
H2, H3, H4, H5, H6 respectively. In the x-z plane (left) migration of the vessel centerline is
taking place towards the left, and in the y-z plane (right) migration is also taking place towards
the left. In this case the general region of the spine is the x-z plane at y = 100 such that the arrow
indicates migration away from the spine.
4.3.2 Eccentricity and Curvature

The results of the eccentricity ratio and 3-point polyline curvature parameter
measurements for the arterial lumen surface of B2 are shown (Figures 48 and 49). It is clear that
eccentricity of cross-sections differ along s in an easily understood manner such that more
elliptical cross-sections have higher eccentricities. The polyline surface curvature is a bit less
92

intuitive. Lower values represent flatter surfaces, and were expected at regions of contact with
the spine. This appears to be present at two points in the proximal posterior and left lateral sides.
There are 4 visible stripes of extreme curvature representing θ = 0, π/2, π, and 3π/2. Due to a
lack of foresight the method used did not account for sinusoidal behavior at these values, though
it very easily could.

1.4
1.35
1.3

1.5

1.25

1.4

1.2

1.3

1.15
1.2

1.1
1.05

1.1

1.0
Figure 48. Eccentricity parameter values displayed for B2 lumen before surface approximation
with a cutoff of 1.4 (left) and after surface approximation (right).

93

0.1
0.08
0.06
0.04
0.02
0
Figure 49. Curvature parameter values displayed for B2 lumen after surface approximation with
a cutoff of mean + 1 standard deviation.
4.4 Discussion
4.4.1 Centerline Generation Algorithm

The centerline generation method yielded results comparable to those found using a trial
of the MedCAD module within Mimics®, a commercial grade centerline generation code. Those
results cannot be shown because access to them was lost when module access expired. The
developed algorithm can handle complex paths, and the number of center points is easily
adjustable. Estimation of an individual maximally inscribed sphere usually showed excellent
convergence within 200 iterations. The centerline generation technique is automated, repeatable,
and relatively insensitive to outward saccular growth. If a different file format than STL was
used, preferably a meshed surface, the top and bottom vessel ends could be more easily excluded
from the algorithm providing a small increase in data near each vessel end as well as rendering
the method insensitive to vessel end surface orientations.

It is important to smoothly

approximate the centerline for the longitudinal plane acquisition method presented. Longitudinal
94

planes are defined normal the centerline, and even a slightly noisy centerline would yield
inaccurate results.

It should be noted that other centerline generation algorithms have been developed for
use with the carotid bifurcation which are based on maximally inscribed spheres, but the
proposed method was developed independently [62–64].

The definition of a centerline is

straightforward and unarguable for ideal geometries such as tubes, but as geometries become
more complex the definition and method of acquiring the centerline becomes more obscure. For
example, a maximally inscribed sphere method is insensitive to extreme saccular growth,
whereas a transverse polyline centroid method would reflect every change in such a growth;
protrusions of the surface inward cause a direct decrease of the size of a maximally inscribed
sphere near that boundary, whereas protrusions outward have no effect. There is no answer as to
which of these methods is correct, and instead the focus should be on which one yields more
meaningful and justifiable results for the purpose of a given study. For example, in the proposed
method the radius is more conservatively defined as the distance from the centerline; the distance
from the centerline can be affected and interpreted in a number of unique ways. If the anterior
side of the vessel expands and the posterior side remains stable, the distance from the centerline
on each side would also increase, which is an important concept to be aware of. As a result the
distance to the centerline can be a misleading parameter unless properly understood.

The initial guess for iteratively fitting a maximally inscribed sphere for a given crosssection also limits the final result. By enforcing a dependency on the previously generated points
of the centerline for the generation of the next point, the overall path of the centerline can be

95

predicted and maintained in an intuitive way. However, this can lead to scenarios where slightly
different inputs yield vastly different results. In order to fully automate centerline generation and
render reliable and repeatable results more advanced measures must be developed.

Such

measures might involve multiple quick passes of centerline generations to accurately predict the
most likely general course of the final centerline to generate using the detailed algorithm. An
additional important consideration is how to accurately compare intra-patient results; each
surface from a longitudinal study has a unique centerline, and therefore each parameterized
surface is parameterized to a different centerline, so the difference in the surfaces are not directly
related with the local change. Nevertheless, this geometric characterization may allow, to some
degree, the separation of local expansion from global shifting (i.e., movement of the centerline)
of the aneurysm shape. A shortcoming of the method presented here then stems from the fact that
individual surface points cannot be individually tracked from one CT surface to the next except a
few identifiable markers.

If two centerlines are generated from two scans for the same patient at different times
they may have different degrees of tortuosity, not to mention aneurysm expansion. Both of these
changes can lead to centerlines with longer total arc lengths, which makes translation of
geometric changes through parameterization difficult. For instance, one surface may have 125
longitudinal planes, while the next may have 140 for the same region. The changes observed are
also regionally based; the proximal or distal end of an aneurysm may see large changes while
another may not, resulting in an uneven change of the spatial distribution of the interpreted
longitudinal planes. Understanding these changes in such a way that their interpretation can be
automated is vital to advancing the utility of longitudinal parameter estimation.

96

4.4.2 Parameterization and Surface Approximation

Another utility of surface parameterization is that various quantities at the aneurismal
wall or on the surface, such as wall thickness, surface curvature, or wall stresses from
biomechanical analyses, can be given as functions of two spatial parameters (s, θ) and time t.
FEA, FSI, or computational fluid simulations could be done with a parameterized arterial surface
mesh and the resulting data would represent both geometrical and mechanical values across the
surface. This will allow future work to utilize techniques such as linear regression or Gaussian
process learning to analyze the effect of various factors on aneurysm expansion and stress and to
predict the time evolution. Many computational biosolid studies use idealized geometries, but
when coupled with this automated surface parameterization, new possibilities emerge to analyze
stress distribution or rupture risk over different geometries.

A major drawback of the proposed method is the reduction of the number of surface
points. For example, the lumen surface of scan B2 before parameterization contains 17,247
points, while after parameterization it contains only 5,106 points. This reduction is due to the
condition that a point must pass a particular inequality to be considered as a point lying on a
given longitudinal plane. Alternatively the s parameter of a given surface point could be
approximated by the minimum distance from the point to the centerline as done by ZeinaliDavarani et al. [60]. This method results in an s value for every point on the surface, but as the
surface strays from an ideal cylinder complications arise to render this method imperfect as well.

97

There are significant drawbacks associated with the degree of reduction seen by the
proposed method. If an FEA, FSI, or fluid simulation was performed and mechanical values are
recorded for each point of the utilized surface mesh, then the reduction of points on the surface
results in a direct loss of data. Less data renders statistical methods less powerful, especially
with respect to surface approximation and correlation analyses. An advanced parameterization
method which preserves every point on the surface and has the speed of the minimum-distance
method as well as the robustness of the longitudinal plane method would be ideal. In order to
accomplish this surface meshes would have to replace surface STLs as the input of choice. The
information a connectivity matrix provides is invaluable for ensuring a reliable algorithm can be
developed.

The surface approximation algorithm itself is quite robust, although some optimization is
required for it to reach its full potential. A test case was run for every combination of the
number of shape functions in each dimension to find an optimum pair. The quality of the result
was based on minimizing the sum of the squared error of the residuals (SSE). The resulting
values are first sorted by Nη and then Nθ for 100 cases to show the general trend of the resulting
SSE values (Figure 50, Table C.1). The resulting values are then sorted by the SSE values
themselves (Figure 51, Table C.2). The important conclusion is that the number shape functions
in each dimension do not evenly contribute to improved accuracy, and that unique optimum
combinations exist on a case by case basis between accuracy and computational time. One goal
of surface approximation is to translate the surface data into a form suitable for rigorous
statistical computations absent from biomechanical rationale to push the limits of what statistics
alone can accomplish before reinforcing it with biomechanical rationale.
98

Shape Function Parameter Selection vs. SSE
60000
50000

SSE

40000

30000
20000
10000
0
1

11

21

31
41
51
61
71
Shape Function Combination #

81

91

Figure 50. SSE values curve sorted by Nη and Nθ for 100 cases of shape function surface fitting.

Shape Function Parameter Selection vs. SSE
60000
50000

SSE

40000

30000
20000
10000
0
1

11

21

31
41
51
61
71
Shape Function Combination #

81

91

Figure 51. SSE values curved sorted by SSE value for 100 case of shape function surface fitting.
4.4.3 Morphological Parameters

The proposed methods for determining parameter values of eccentricity and polyline
curvature are inadequate. The measure of eccentricity suffers from major drawbacks: there is no
constraint that the major and minor axes intersect at the given cross-sections origin, the axes are
99

only approximately 90° from one another, the direction of the major axis is subject to surface
irregularities, and there is no continuity of axes from one cross-section to the next. All of these
problems stem from the fact that point-cloud data was used instead of surface data; by
constraining the eccentricity/asymmetry axes to pre-defined unstructured points along the
polyline of a cross-section there is no way to remedy the drawbacks listed. In order to properly
define the eccentricity axes a more structured distribution of data for a given cross-section is
required. If the pre-fit and post-fit values are compared numerically it is clear how the results are
drastically improved with clear continuity of eccentricity trends across cross-sections for the
ideally structured post-fit surface (Figure 52). The following conditions are able to be enforced
for the post-fit surface that cannot be enforced for the pre-fit surface: the major and minor axes
intersect at a given cross-section’s origin, the axes are exactly 90° from one another, the surface
is smooth rendering the direction of the major axis less influenced by surface irregularities, and
there is a natural continuity and progression of eccentricity across cross-sections as a virtue of
utilizing normal shape functions. The proposed method of eccentricity calculation is therefore
only really meaningful when paired with the post-fit surface of a given AAA.

100

Figure 52. Eccentricity of the arterial lumen surface of B2 pre-fit (blue) and post-fit (red).

The second analyzed morphological parameter, the polyline curvature, is also limited in
its interpretation. First and foremost a small revision to the method would eliminate the highvalue stripes seen at θ = 0, π/2, π, and 3π/2. After such a fix it must also be critically analyzed as
to the best set of 3 points to use to determine the curvature of a given point. For the results
shown a set of a given point and the points immediately neighboring it on either side are used for
the curvature calculation. The selection of the additional points for the calculation affects results
and should be considered very carefully. The results shown are for a post-fit approximated
surface, which is mathematically smooth. The polyline curvature method yields nonsensical
results for the pre-fit surface due to surface spikes and irregularities and no discernable trends
101

could be identified. The selection of curvature calculation points for the pre-fit data is extremely
difficult to determine to yield meaningful results.

It is clear that both of the morphological parameters described are not easily obtained for
pre-fit surfaces. Whether this highlights the importance of a proper approximation algorithm or
the shortcomings of the methods themselves is hard to distinguish, but it is clear the results show
more intuitively understood trends for post-fit surfaces.

The post-fit surfaces are still

approximations, however, so the results have to be presented with that caveat. It is believed that
the surface approximation presented is worthwhile and actually enhances the readability and
reliability of results for the methods proposed.
4.4.4 Future Work

The algorithms presented are not finalized and this study has many possibilities for
expansion to serve as a useful tool for clinical management of AAAs or computational
biomechanical studies. There exist many places to further automate and optimize the existing
methods, but some methods presented may benefit from a total redesign. For example, it must
be clarified that eccentricity and asymmetry are intertwined but distinctly different. Eccentricity
is simply a measure of the skewedness of a given cross-section, but the asymmetry can be best
analyzed by considering the directions of the major and minor axes utilized to calculate
eccentricity.

The currently proposed method is highly influenced by surface irregularities.

Alternatives exist such as the principal axes method and calculating the second moment of inertia
of a cross-section. These methods are not as heavily influenced by surface irregularities and
carry a more two-dimensionally robust justification than simply identifying the maximum

102

diameter pair of points constrained to the origin of a cross-section. Likewise the polyline
curvature method presented is simplistic and an increase in dimensionality would be ideal.
Considering the 3-D curvature of the surface instead of polyline curvature constrained to
individual planes would offer a more intuitive global understanding of the curvature across the
aneurysm and provide a more robust solution for unique AAA surfaces.

Aside from the various improvements that could be done to the existing parameterization
method proposed, a lot still must be done to fully bridge the gap between AAA pathology and
biomechanics. The consideration of the interactions with surrounding tissues of AAAs will serve
to enrich the biomechanical understanding of AAA pathology. Other biomechanical factors still
exist of course; for example, it is not understood why different demographics such as women vs.
men exhibit much different risk for AAA development. Ongoing studies have been examining
the unique blood flow due to anatomical differences between men and women in the abdomen
and suggest that hemodynamics may play a role in AAA development. Other studies have
suggested that hormones in the body and bloodstream may contribute to AAA risk as well. At
this point it is clear that not only do the engineering and medical communities need to cooperate
but biochemistry, statistics, and biology experts should become involved to further the
understanding of AAAs. The mechanical understanding and simulation of AAAs has reached a
critical point that must be overcome. The logical next step for the research community is the
widespread adoption of longitudinal datasets to understand the mechanical and geometrical
evolution of AAAs over time considering mechanical, statistical, biological, and biochemical
factors between AAAs and surrounding tissues.

103

APPENDICES

104

APPENDIX A

Detailed information and measurements of CT scan data

105

Table A.1. Detailed scan resolution information for the CT data set studied.
Scan
A1
A2
B1
B2
B3
C1
C2
D1
D2
D3
E1
E2
F1
F2
F3
F4
F5
F6
F7
G1
G2
G3
G4
G5
G6
H1
H2
H3
H4
H5
H6
H7
H8
H9
I1
I2
I3
I4
I5

Pixel Size (mm)
0.656
0.645
0.680
0.678
0.777
0.688
0.684
0.605
0.635
0.578
0.617
0.572
0.477
0.375
0.510
0.531
0.545
0.539
0.656
0.783
0.670
0.730
0.639
0.684
0.633
0.494
0.547
0.723
0.678
0.688
0.652
0.625
0.703
0.684
0.547
0.625
0.645
0.676
0.549

# Slices
1003
618
931
951
1073
961
684
375
596
621
863
899
124
133
441
381
411
730
588
2096
1020
805
510
642
463
249
409
693
666
653
640
504
505
453
468
653
287
659
443

Slice increment (mm)
var.
var.
var.
0.699
var.
var.
1.000
1.499
1.000
0.699
0.699
0.699
2.000
2.000
1.249
var.
1.500
0.699
1.000
var.
0.699
0.699
1.000
1.000
var.
var.
var.
var.
0.699
0.699
0.699
1.000
1.000
var.
1.250
1.000
var.
0.699
1.000

106

Slice thickness (mm)
1.000
1.000
1.000
1.000
1.000
1.000
1.000
2.000
1.250
1.000
1.000
1.000
3.000
3.000
1.250
3.200
2.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
3.200
3.200
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.250
1.250
3.200
1.000
1.000

Table A.1 (cont’d)
I6
J1
J2
J3
J4
J5
K1
K2
K3
K4
K5
K6

0.582
0.684
0.621
0.658
0.650
0.652
0.742
0.793
0.637
0.660
0.758
0.750

608
732
673
650
456
654
926
701
717
995
1086
787

0.699
var.
0.699
var.
1.000
0.699
0.699
1.000
var.
0.699
0.699
0.699

1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.250
1.000
1.000
1.000
1.000

Table A.2. Detailed scanner manufacturers and models, date of scans, and total time elapsed
from initial scanning.
Scan
A1
A2
B1
B2
B3
C1
C2
D1
D2
D3
E1
E2
F1
F2
F3
F4
F5
F6
F7
G1
G2
G3
G4
G5
G6
H1
H2

Manufacturer
SIEMENS
Philips
SIEMENS
SIEMENS
SIEMENS
SIEMENS
Philips
SIEMENS
GE_MEDICAL_SYSTEMS
SIEMENS
SIEMENS
SIEMENS
SIEMENS
SIEMENS
GE_MEDICAL_SYSTEMS
Philips
SIEMENS
SIEMENS
Philips
SIEMENS
SIEMENS
SIEMENS
Philips
Philips
Philips
Marcon
Philips

Model
Definition
Brilliance 64
Definition
Definition
SOMATOM Definition
Definition
Brilliance 64
Sensation 16
LightSpeed Ultra
Definition
Definition
Definition
SOMATOM PLUS 4
SOMATOM PLUS 4
LightSpeed Ultra
Mx8000
Sensation 16
Sensation 16
Brilliance 64
Sensation 16
Sensation 16
Sensation 16
Brilliance 64
Brilliance 64
Brilliance 64
Mx8000
Mx8000
107

Date
9/11/2009
3/12/2010
10/30/2008
5/21/2009
11/2/2010
3/9/2009
12/1/2010
7/31/2004
12/7/2005
4/18/2008
4/4/2008
10/13/2008
12/15/1999
3/9/2001
12/2/2003
2/24/2004
8/17/2004
12/2/2005
3/6/2007
4/20/2005
10/13/2005
4/13/2006
10/9/2006
4/4/2007
10/11/2007
10/25/2001
11/18/2002

Elapsed
0
0 y, 6 m, 1 d
0
0 y, 6 m, 21 d
2 y, 0 m, 3 d
0
1 y, 6 m, 22 d
0
1 y, 4 m, 7 d
3 y, 8 m, 18 d
0
0 y, 6 m, 9 d
0
1 y, 2 m, 22 d
3 y, 11 m, 17 d
4 y, 2 m, 9 d
4 y, 8 m, 2 d
5 y, 11 m, 17 d
7 y, 2 m, 19 d
0
0 y, 5 m, 23 d
0 y, 11 m, 24 d
1 y, 5 m, 19 d
1 y, 11 m, 15 d
2 y, 5 m, 21 d
0
1 y, 0 m, 24 d

Table A.2 (cont’d)
H3
H4
H5
H6
H7
H8
H9
I1
I2
I3
I4
I5
I6
J1
J2
J3
J4
J5
K1
K2
K3
K4
K5
K6

SIEMENS
SIEMENS
SIEMENS
SIEMENS
Philips
Philips
Philips
GE_MEDICAL_SYSTEMS
GE_MEDICAL_SYSTEMS
Philips
SIEMENS
Philips
SIEMENS
SIEMENS
SIEMENS
SIEMENS
Philips
SIEMENS
SIEMENS
GE_MEDICAL_SYSTEMS
N/A
N/A
N/A
N/A

Sensation 16
Definition
Definition
Sensation 16
Brilliance 64
Brilliance 64
Brilliance 64
LightSpeed Ultra
LightSpeed Ultra
Mx8000
Sensation 16
Brilliance 64
Definition
Sensation 16
Sensation 16
Definition
Brilliance 64
SOMATOM Definition
Sensation 16
LightSpeed Ultra
N/A
N/A
N/A
N/A

108

7/28/2007
7/9/2008
6/29/2009
6/16/2010
12/2/2010
3/9/2011
3/19/2011
8/26/2003
9/3/2004
9/9/2005
8/4/2006
8/3/2007
7/1/2009
11/2/2006
11/23/2007
11/28/2008
11/27/2009
5/15/2010
10/5/2005
5/20/2006
8/10/2007
8/19/2008
8/8/2009
12/1/2009

5 y, 9 m, 3 d
6 y, 8 m, 14 d
7 y, 8 m, 4 d
8 y, 7 m, 22 d
9 y, 1 m, 7 d
9 y, 4 m, 12 d
9 y, 4 m, 22 d
0
1 y, 0 m, 8 d
2 y, 0 m, 14 d
2 y, 11 m, 9 d
3 y, 11 m, 8 d
5 y, 10 m, 5 d
0
1 y, 0 m, 21 d
2 y, 0 m, 26 d
3 y, 0 m, 25 d
3 y, 6 m, 13 d
0
0 y, 7 m, 15 d
1 y, 10 m, 5 d
2 y, 10 m, 14 d
3 y, 10 m, 3 d
4 y, 1 m, 26 d

Table A.3. Tissue Measurements from individual patient CT scans.

Time t
Dmax

Scan 4

Scan 5

Scan 6

Scan 7

H
0

I
0

J
0

K
0

52.02

45.03

44.85

31.68

40.20

42.22

29.82
632

49.59
494

450

40.42
176

40.81
389

29.81
374

39.29
386

42.18
227

45.37

33.91

59.21

44.69

49.01

48.80

30.32

42.69

43.46

41.68

32.25

55.68
1357

42.74
1448

41.47
358

44.78
2102

28.82
745

41.50
757

42.42
674

54.53

72.10

58.79

49.90

56.14

35.02

44.09

45.68

54.58

70.03

56.46
1532

43.78
537

50.30
2449

32.28
1074

43.79
1121

42.22
1049

62.41

52.18

59.23

36.16

46.74

53.00

90° to Dmax
Time t

56.96
1707

49.84
714

52.87
2804

34.28
1438

45.67
1290

51.94
1403

62.73

53.38

61.74

37.78

50.07

57.53

90° to Dmax
Time t
Dmax

61.92
2179

49.73
904

54.49
3156

36.61
2136

49.30

55.24

69.11

58.08

67.28

41.63

90° to Dmax
Time t

67.92
2638

52.08

59.53
3325

37.98

Dmax

Scan 3

G
0

Dmax

Scan 2

B
0

C
0

C (sac)
0

D
0

57.31

46.21

41.11

31.82

90° to Dmax
Time t

50.43
182

44.21
203

39.53
632

Dmax

Scan 1

A
0

58.81

50.06

90° to Dmax
Time t
Dmax

55.42

44.31
733

73.10

69.20

90° to Dmax

72.23

60.22

90° to Dmax
Time t
Dmax

109

F
0

Table A.4. Lumen Measurements from individual patient CT scans.

Time t
Dmax

Scan 4

Scan 5

Scan 6

Scan 7

H
0

I
0

J
0

K
0

45.88

39.87

38.48

28.86

35.62

35.39

24.53
632

40.12
494

450

33.94
176

36.26
389

27.08
374

33.23
386

27.29
227

39.80

29.51

50.85

40.81

40.80

40.81

30.32

35.27

38.98

38.40

26.05

32.29
1357

39.33
1448

37.03
358

39.92
2102

28.82
745

33.78
757

28.74
674

52.37

64.69

48.48

43.72

40.77

32.14

37.65

38.32

49.34

45.40

36.43
1532

38.52
537

31.90
2449

28.80
1074

33.04
1121

30.27
1049

49.09

47.54

33.97

33.67

37.92

41.24

90° to Dmax
Time t

37.52
1707

43.80
714

28.05
2804

32.17
1438

36.05
1290

28.91
1403

49.58

50.22

38.58

35.12

41.01

39.94

90° to Dmax
Time t
Dmax

39.56
2179

45.25
904

31.64
3156

33.47
2136

39.54

29.32

44.17

51.20

41.18

37.55

90° to Dmax
Time t

39.96
2638

43.95

34.22
3325

34.54

Dmax

Scan 3

G
0

Dmax

Scan 2

B
0

C
0

C (sac)
0

D
0

54.24

38.36

36.05

25.38

90° to Dmax
Time t

47.91
182

37.73
203

34.05
632

Dmax

Scan 1

A
0

56.99

44.40

90° to Dmax
Time t
Dmax

52.46

39.65
733

48.89

45.30

90° to Dmax

40.76

36.35

90° to Dmax
Time t
Dmax

110

F
0

APPENDIX B

MatLab code developed for automatic parameterization of AAAs

111

Appendix B.1. step_0.m, part of the Matlab code written to generate the centerlines
clear all; close all; clc;
a = dir('*.txt');
global A B i Xs Xstan s N1 N2 N3 perp2 meanR radius circ theta C D num_data
for i = 1:length(a)
A{i,1} = load(a(i).name); A{i,2} = a(i).name(1:2); A{i,3} = a(i).name(4);
B{i,2} = A{i,2}; B{i,3} = A{i,3};
end
for i = 1:length(a)

%

step_1;
step_2;
% step_3;
% step_4;
% step_5;
% step_6;

end

112

Appendix B.2. step_1, the core functionality of the centerline generation algorithm
%% script for the centerline generation
min_z = min(A{i,1}(:,3));
max_z = max(A{i,1}(:,3));
k = 0;
for j = 1:length(A{i,1})
if A{i,1}(j,3) > (min_z + 2.5) && A{i,1}(j,3) < (max_z - 2.5)
k = k + 1;
A_temp(k,1) = A{i,1}(j,1); A_temp(k,2) = A{i,1}(j,2); A_temp(k,3) = A{i,1}(j,3);
end
end
A{i,1} = A_temp;
itenum = 250;
S_n = 2;
min_z = min(A{i,1}(:,3)); max_z = max(A{i,1}(:,3));
S = 0;
S_chk = max_z - 1;
while S_chk < max_z
S = S + 1;
htext = sprintf(['Calculating Center Point for s = ',num2str(S),' ...']);
clear A2 d CoT k; k = 0;
if S == 1
for l = 1:length(A{i,1});
if A{i,1}(l,3) >= min_z && A{i,1}(l,3) <= min_z + 1
k = k+1;
A2(k,:) = A{i,1}(l,:);
end
end

113

Appendix B.2 (cont’d)
B{i,1}(S,1) = sum(A2(:,1))/length(A2); B{i,1}(S,2) = sum(A2(:,2))/length(A2); B{i,1}(S,3)
= min_z;
elseif S == 2
B{i,1}(S,1) = B{i,1}(S-1,1); B{i,1}(S,2) = B{i,1}(S-1,2); B{i,1}(S,3) = B{i,1}(S-1,3) +
S_n;
else
PQ_S1 = [B{i,1}(S-1,1)-B{i,1}(S-2,1), B{i,1}(S-1,2)-B{i,1}(S-2,2), B{i,1}(S-1,3)B{i,1}(S-2,3)];
PQ_S2 = sqrt((B{i,1}(S-1,1)-B{i,1}(S-2,1))^2 + (B{i,1}(S-1,2)-B{i,1}(S-2,2))^2 +
(B{i,1}(S-1,3)-B{i,1}(S-2,3))^2);
PQ_S = PQ_S1 / PQ_S2;
B{i,1}(S,1) = B{i,1}(S-1,1) + S_n*PQ_S(1); B{i,1}(S,2) = B{i,1}(S-1,2) + S_n*PQ_S(2);
B{i,1}(S,3) = B{i,1}(S-1,3) + S_n*PQ_S(3);
end
itesize(S) = 1;
h = waitbar(0,htext);
for j = 1:itenum
for k = 1:length(A{i,1})
d(k) = sqrt((A{i,1}(k,1)-B{i,1}(S,1))^2 + (A{i,1}(k,2)-B{i,1}(S,2))^2 + (A{i,1}(k,3)B{i,1}(S,3))^2);
end
[d_X,d_Y] = min(d);
[X,Y,Z] = sphere(20); %Makes a dummy unit sphere
X = d_X*X; Y = d_X*Y; Z = d_X*Z; %Makes the unit sphere have radius = minmum
distance above
X = X + B{i,1}(S,1); %Moves the sphere's center to the centroid
Y = Y + B{i,1}(S,2);
Z = Z + B{i,1}(S,3);
PQ_1 = [A{i,1}(d_Y,1)-B{i,1}(S,1), A{i,1}(d_Y,2)-B{i,1}(S,2), A{i,1}(d_Y,3)-B{i,1}(S,3)];
114

Appendix B.2 (cont’d)
PQ_2 = sqrt((A{i,1}(d_Y,1)-B{i,1}(S,1))^2 + (A{i,1}(d_Y,2)-B{i,1}(S,2))^2 +
(A{i,1}(d_Y,3)-B{i,1}(S,3))^2);
PQ = PQ_1 / PQ_2;
if S==1
B{i,1}(S,1) = B{i,1}(S,1) - PQ(1)*itesize(S); B{i,1}(S,2) = B{i,1}(S,2) PQ(2)*itesize(S);
elseif S==2
B{i,1}(S,1) = B{i,1}(S,1) - PQ(1)*itesize(S); B{i,1}(S,2) = B{i,1}(S,2) PQ(2)*itesize(S);
else
proj_PQ = cross(PQ_S,(cross(PQ,PQ_S)/norm(PQ_S)))/norm(PQ_S);
proj_PQ = proj_PQ /norm(proj_PQ);
B{i,1}(S,:) = B{i,1}(S,:) - proj_PQ*itesize(S);
end
x_c2(j,S) = B{i,1}(S,1); y_c2(j,S) = B{i,1}(S,2); z_c2(j,S) = B{i,1}(S,3);
d2(j,S) = d_X;
if j >= 15
CoT(j-14,S) = sqrt((x_c2(j,S)-x_c2(j-14,S))^2 + ((y_c2(j,S)-y_c2(j-14,S))^2));
if CoT(j-14,S) < itesize(S)
itesize(S) = itesize(S) / 2;
end
end
%
if j==1
%
figure;scatter3(A{i,1}(:,1),A{i,1}(:,2),A{i,1}(:,3),1,'CData',A{i,1}(:,3),'marker','.'); hold on
%Plots the AAA colorfully
%
plot3(B{i,1}(S,1),B{i,1}(S,2),B{i,1}(S,3),'g.','Markersize',15); %Plots the
centroid as a green dot
%
title(A(i,2));xlabel('x direction (mm)');ylabel('y direction (mm)');zlabel('z
direction (mm)');axis image %labels the graph
%
hold
on;plot3(A{i,1}(d_Y,1),A{i,1}(d_Y,2),A{i,1}(d_Y,3),'b.','Markersize',20) %Marks the
minimum distance point with a blue dot
%
plot3(X,Y,Z,'ro');axis image %Shows the sphere before any iterations
if S~=1 && S ~= 2 && j==itenum
figure;scatter3(A{i,1}(:,1),A{i,1}(:,2),A{i,1}(:,3),1,'CData',A{i,1}(:,3),'marker','.');
hold on %Plots the AAA colorfully
115

Appendix B.2 (cont’d)
plot3(B{i,1}(S,1),B{i,1}(S,2),B{i,1}(S,3),'g.','Markersize',15); %Plots the centroid
as a green dot
title(A(i,2));xlabel('x direction (mm)');ylabel('y direction (mm)');zlabel('z direction
(mm)');axis image %labels the graph
hold on;plot3(A{i,1}(d_Y,1),A{i,1}(d_Y,2),A{i,1}(d_Y,3),'b.','Markersize',20)
%Marks the minimum distance point with a blue dot
plot3(X,Y,Z,'ro');axis image %Shows the sphere after all iterations
figure;
plot3(B{i,1}(S,1),B{i,1}(S,2),B{i,1}(S,3),'g.','Markersize',15); %Plots the centroid
as a green dot
title(A(i,2));xlabel('x direction (mm)');ylabel('y direction (mm)');zlabel('z direction
(mm)');axis image %labels the graph
hold on;plot3(A{i,1}(d_Y,1),A{i,1}(d_Y,2),A{i,1}(d_Y,3),'b.','Markersize',20)
%Marks the minimum distance point with a blue dot
scatter3(X(:),Y(:),Z(:),'Cdata',Z(:),'marker','.');axis image %Shows the sphere after
all iterations
plot3([B{i,1}(S,1),B{i,1}(S,1)+proj_PQ(1)*d_X],[B{i,1}(S,2),B{i,1}(S,2)+proj_PQ(2)*d_X],[B
{i,1}(S,3),B{i,1}(S,3)+proj_PQ(3)*d_X],'r','Markersize',15); %graphs projection vector
plot3([B{i,1}(S,1),B{i,1}(S,1)+PQ(1)*d_X],[B{i,1}(S,2),B{i,1}(S,2)+PQ(2)*d_X],[B{i,1}(S,3),
B{i,1}(S,3)+PQ(3)*d_X],'b','Markersize',15); %graphs min_dist vector
plot3([B{i,1}(S,1),B{i,1}(S,1)+PQ_S(1)*d_X],[B{i,1}(S,2),B{i,1}(S,2)+PQ_S(2)*d_X],[B{i,1}
(S,3),B{i,1}(S,3)+PQ_S(3)*d_X],'k','Markersize',15); %graphs normal vector
end
waitbar(j / itenum,h,htext)
S_chk = B{i,1}(S,3);
end
delete(h)
end
disp(['Refining Final S value: ',num2str(S),' ...'])
clear d d_X d_Y CoT h htext
B{i,1}(S,1) = B{i,1}(S-1,1) + S_n*PQ_S(1);
B{i,1}(S,2) = B{i,1}(S-1,2) + S_n*PQ_S(2);
B{i,1}(S,3) = max_z;

116

Appendix B.2 (cont’d)
itesize(S) = 1;
for j = 1:itenum
Appendix 6 (cont’d)
for k = 1:length(A{i,1})
d(k) = sqrt((A{i,1}(k,1)-B{i,1}(S,1))^2 + (A{i,1}(k,2)-B{i,1}(S,2))^2 + (A{i,1}(k,3)B{i,1}(S,3))^2);
end
[d_X,d_Y] = min(d);
PQ_1 = [A{i,1}(d_Y,1)-B{i,1}(S,1), A{i,1}(d_Y,2)-B{i,1}(S,2), A{i,1}(d_Y,3)-B{i,1}(S,3)];
PQ_2 = sqrt((A{i,1}(d_Y,1)-B{i,1}(S,1))^2 + (A{i,1}(d_Y,2)-B{i,1}(S,2))^2 +
(A{i,1}(d_Y,3)-B{i,1}(S,3))^2);
PQ = PQ_1 / PQ_2;
B{i,1}(S,1) = B{i,1}(S,1) - PQ(1)*itesize(S); B{i,1}(S,2) = B{i,1}(S,2) - PQ(2)*itesize(S);
x_c2(j,S) = B{i,1}(S,1); y_c2(j,S) = B{i,1}(S,2); z_c2(j,S) = B{i,1}(S,3);
d2(j,S) = d_X;
if j >= 15
CoT(j-14,S) = sqrt((x_c2(j,S)-x_c2(j-14,S))^2 + ((y_c2(j,S)-y_c2(j-14,S))^2));
if CoT(j-14,S) < itesize(S)
itesize(S) = itesize(S) / 2;
end
end
disp(['S value: ',num2str(S),' ... '])
disp(['Iteration #: ',num2str(j),' / ',num2str(itenum)])
end
figure;scatter3(A{i,1}(:,1),A{i,1}(:,2),A{i,1}(:,3),1,'CData',A{i,1}(:,3),'marker','.');
hold on; axis image;
scatter3(B{i,1}(:,1),B{i,1}(:,2),B{i,1}(:,3));
plot3(B{i,1}(:,1),B{i,1}(:,2),B{i,1}(:,3));
clear A_temp CoT PQ proj_PQ PQ_1 PQ_2 PQ_S PQ_S1 PQ_S2 S S_chk _n X Y Z d d2 d_X
d_Y itenum itesize j k l max_z min_z x_c2 y_c2 z_c2

117

Appendix B.3. step_2.m of the centerline generation algorithm This step generates smooth
points along the centerline by an approximation function
X = B{i,1}';
Lbp=zeros(1,length(X)-1); %Length between points: calculates distance between each point.
Lbp(1)=distance between X(1) and X(2), etc.
for j=1:(length(X)-1)
Lbp(j)=(((X(1,(j+1))-X(1,j))^2)+((X(2,(j+1))-X(2,j))^2)+((X(3,(j+1))-X(3,j))^2))^(1/2);
end
B{i,4} = sum(Lbp); % Length of the centerline X or B{i,1}
l=zeros(1,length(X)); %the distance along the centerline at each point. l(1) is the distance from
X(:,1) to X(:,1). l(2) is the distance from X(:,1) to X(:,2), etc.
for j=1:length(Lbp)
l(j+1)=sum(Lbp(1:j));
end
phi=zeros(1,(length(l)));%a group of base functions that depend on s. Use phis functions to fit a
centerline.
dphi_ds=zeros(1,length(l)); %derivative of phis with respect to s. Used to find tangent vectors at
each s for use as a reference vector later.
Xsarray=zeros(1,length(X),3); %contains the x,y,z elements that need to be summed to create
Xs. It is the result of X.*nphis (nphis is a normalization of the phis variable).
Xstanarray=zeros(1,length(X),3); %contains the x,y,z elements of X.*dnphidss.
s=l(1):1:l(end); %longitudinal distance along the centerline
Xs=zeros(1,length(s),3); %function of s that defines the coordiantes of the centerline at any s. It
is a smooth 4th order approximation function of the centerline.
Xstan=zeros(1,length(s),3); %contains the x,y,z coordinates that, at each s, will create a vector
from the origin to the aforementioned coordinates parrallel to the tangent vector at that s.
j = length(l);
for k=1:length(s)
for m=1:j
if m == 1
if ( s(k) >= l(1) ) && ( s(k) <= l(3) )
phi(m) = (s(k)-l(3))^2*(s(k)+l(3))^2/((l(1)-l(3))^2*(l(1)+l(3))^2);
dphi_ds(m) = (2*(s(k)-l(3))*(s(k)+l(3))^2+2*(s(k)-l(3))^2*(s(k)+l(3)))/((l(1)l(3))^2*(l(1)+l(3))^2);
else
phi(m) = 0.0;
dphi_ds(m)=0.0;
end
118

Appendix B.3 (cont’d)
elseif m == 2
if ( s(k) >= l(1) ) && ( s(k) <= l(4) )
phi(m) = (s(k)+l(2))^2*(s(k)-l(4))^2/((l(2)+l(2))^2*(l(2)-l(4))^2);
dphi_ds(m) = (2*(s(k)+l(2))*(s(k)-l(4))^2+2*(s(k)+l(2))^2*(s(k)l(4)))/((l(2)+l(2))^2*(l(2)-l(4))^2);
else
phi(m) = 0.0;
dphi_ds(m)=0.0;
end
elseif m == j-1
if ( s(k) >= l(j-3) ) && ( s(k) <= l(j) )
phi(m) = (s(k)-l(j-3))^2*(s(k)-2*l(j)+l(j-1))^2/((l(j-1)-l(j-3))^2*(l(j-1)-l(j))^2);
dphi_ds(m) = (2*(s(k)-l(j-3))*(s(k)-2*l(j)+l(j-1))^2+2*(s(k)-l(j-3))^2*(s(k)-2*l(j)+l(j1)))/((l(j-1)-l(j-3))^2*(l(j-1)-l(j))^2);
else
phi(m) = 0.0;
dphi_ds(m)=0.0;
end
elseif m == j
if ( s(k) >= l(j-2) ) && ( s(k) <= l(j) )
phi(m) = (s(k)-l(j-2))^2*(s(k)-2*l(j)+l(j-2))^2/((l(j)-l(j-2))^4);
dphi_ds(m) = (2*(s(k)-l(j-2))*(s(k)-2*l(j)+l(j-2))^2+2*(s(k)-l(j-2))^2*(s(k)-2*l(j)+l(j2)))/((l(j)-l(j-2))^4);
else
phi(m) = 0.0;
dphi_ds(m)=0.0;
end
else
if ( s(k) >= l(m-2) ) && ( s(k) <= l(m+2) )
phi(m) = (s(k)-l(m-2))^2*(s(k)-l(m+2))^2/((l(m)-l(m-2))^2*(l(m)-l(m+2))^2);
dphi_ds(m) = (2*(s(k)-l(m-2))*(s(k)-l(m+2))^2+2*(s(k)-l(m-2))^2*(s(k)l(m+2)))/((l(m)-l(m-2))^2*(l(m)-l(m+2))^2);
else
phi(m) = 0.0;
dphi_ds(m)=0.0;
end
end
end
n_phi=phi/(sum(phi)); %Normalize
n_dphi_ds=((dphi_ds*sum(phi))-(phi*sum(dphi_ds)))/(sum(phi)^2); %normalize
for m=1:j
119

Appendix B.3 (cont’d)
for n=1:3
Xsarray(1,m,n)=X(n,m)*n_phi(m);
Xstanarray(1,m,n)=X(n,m)*n_dphi_ds(m);
end
end
for n=1:3
Xs(1,k,n)=sum(Xsarray(1,:,n));
Xstan(1,k,n)=sum(Xstanarray(1,:,n));
end
end
for j=1:(length(Xs)-1)
Lbp2(j)=(((Xs(1,(j+1),1)-Xs(1,j,1))^2)+((Xs(1,(j+1),2)-Xs(1,j,2))^2)+((Xs(1,(j+1),3)Xs(1,j,3))^2))^(1/2);
end
Lbp3 = (((Xs(1,length(Xs),1)-Xs(1,1,1))^2)+((Xs(1,length(Xs),2)Xs(1,1,2))^2)+((Xs(1,length(Xs),3)-Xs(1,1,3))^2))^(1/2);
tort(i) = sum(Lbp2(:)) / Lbp3;
figure('Color',[1 1
1]);scatter3(A{i,1}(:,1),A{i,1}(:,2),A{i,1}(:,3),1,'CData',A{i,1}(:,3),'marker','.');
hold on;
% scatter3(B{i,1}(:,1),B{i,1}(:,2),B{i,1}(:,3));
% plot3(B{i,1}(:,1),B{i,1}(:,2),B{i,1}(:,3));
plot3(Xs(:,:,1),Xs(:,:,2),Xs(:,:,3),'k','LineWidth',2);
tort_string = sprintf('Tortuosity = %f',tort(i));
title(tort_string);
axis image;
B{i,5}(:,:) = Xs(1,:,:);
clear j k l m n Lbp Lbp2 Lbp3 dphi_ds n_dphi_ds n_phi phi

120

Appendix B.4. step_3.m, the section of the Matlab code which calculates the normal vectors
%% Find N1, N2, N3 and perpendicular vector planes at each point s along the centerline.
%% Normal to centerline basis
%N2: is a unit vector that is a projection of the x axis unit vector perpendicular to vector N1.
%N3: is a unit vector that is the result of the cross product of N1 X N2.
N1=zeros(1,length(s),3); %Xstan vectors as unit vectors: the unit tangent vector at each s along
the centerline approximation.
for j=1:length(Xstan)
for k=1:3
N1(1,j,k)= Xstan(1,j,k)/(sqrt((Xstan(1,j,1)^2)+(Xstan(1,j,2)^2)+(Xstan(1,j,3)^2)));
end
end
graphN1=zeros(1,2,3); %A vector from the nth point along s, to the coordinates of the N1 vector
added to the nth point along s.
figure;
for j=1:length(Xs)
for k=1:3
graphN1(1,:,k)=[Xs(1,j,k),(Xs(1,j,k)+N1(1,j,k))];
end
plot3(graphN1(1,:,1),graphN1(1,:,2),graphN1(1,:,3),'r');
hold on
end
%% N2 Calculations: the projection of the x axis onto the perpendicular plane at point s along
the centerline
k=zeros(1,2,3); %serves as a unit vector in the x direction.
k(1,2,1)=1;
N2_num = zeros(1,length(N1),3); %the numerator of the equation to compute N2
N2_den = zeros(1,length(N1)); %the denominator of the equation needed to compute N2
N2=zeros(1,length(N1),3);
for j=1:length(N1) %This for loop calculates N2 for each s
for m=1:3
N2_num(1,j,m)=(k(1,2,m)(((k(1,2,1)*N1(1,j,1))+(k(1,2,2)*N1(1,j,2))+(k(1,2,3)*N1(1,j,3)))*N1(1,j,m)));
121

Appendix B.4 (cont’d)
end
N2_den(j)=sqrt((N2_num(1,j,1)^2)+(N2_num(1,j,2)^2)+(N2_num(1,j,3)^2));
for m=1:3
N2(1,j,m)=N2_num(1,j,m)/N2_den(j);
end
end
for j=length(N2) %makes N2 vectors into unit vectors
N2_den(j)=sqrt((N2(1,j,1)^2)+(N2(1,j,2)^2)+(N2(1,j,3)^2));
for m=1:3
N2(1,j,m)=N2(1,j,m)/N2_den(j);
end
end
for j=1:length(Xs) %Graphs each N2 at it's respective s
for m=1:3
graphN2(1,:,m)=[Xs(1,j,m),(Xs(1,j,m)+N2(1,j,m))];
end
plot3(graphN2(1,:,1),graphN2(1,:,2),graphN2(1,:,3),'b')
end
%% N3 Calculations
N3=zeros(1,length(N2),3); %the result of the cross product between N1 and N2
for j=1:length(N2) %calculates N3
N3(1,j,1)=(N1(1,j,2)*N2(1,j,3))-(N1(1,j,3)*N2(1,j,2));
N3(1,j,2)=-((N1(1,j,1)*N2(1,j,3))-(N1(1,j,3)*N2(1,j,1)));
N3(1,j,3)=(N1(1,j,1)*N2(1,j,2))-(N1(1,j,2)*N2(1,j,1));
end
for j=1:length(Xs) %graphs N3 at it's respective s
for m=1:3
graphN3(1,:,m)=[Xs(1,j,m),(Xs(1,j,m)+N3(1,j,m))];
end
plot3(graphN3(1,:,1),graphN3(1,:,2),graphN3(1,:,3),'g')
end
%Manipulates and labels the first figure properly.
grid on
axis image
title('Approximated Centerline with N1, N2, and N3 Vectors at each point s','fontsize',12)
122

Appendix B.4 (cont’d)
xlabel('X direction (mm)')
ylabel('Y direction (mm)')
zlabel('Z direction (mm)')
%% Perpindicular vector determination
perp=zeros(length(A{i,1}),length(s),3); %a structured array that contains x,y,z coordinates of
perpendicular vectors going from a designated s to a known point on the surface
for j=1:length(N1) %Performs the dot product between the tangent vector at each s with the
vector going from that point s to each known surface point to determine where the dot product is
close to zero. This indicates the vectors are perpendicular and the information for that vector is
stored in perp.
for m=1:length(A{i,1})
k=[A{i,1}(m,1)-Xs(1,j,1), A{i,1}(m,2)-Xs(1,j,2),A{i,1}(m,3)-Xs(1,j,3)];
if norm(k(3))<15
dotprod=(N1(1,j,1)*k(1))+(N1(1,j,2)*k(2))+(N1(1,j,3)*k(3));
h=norm(k);
if dotprod/h<.01 && dotprod/h>-.01
for n=1:3
perp(m,j,n)=A{i,1}(m,n)-Xs(1,j,n);
end
end
end
end
end
for j=1:length(s) %Eliminates all of the dot products that did not result in perpendicular vectors
to vastly decrease the size of the perp vector and store it as perp2
k=0;
for m=1:length(A{i,1})
if perp(m,j,:)~=0
k=k+1;
for n=1:3
perp2(k,j,n)=perp(m,j,n);
end
else
end
end
end
k=size(perp2);
figure
123

Appendix B.4 (cont’d)
for m=1:k(1)%Graphs the perpendicular perp2 vectors at their respective s
for j=1:k(2)
if perp2(m,j,1)==0 && perp2(m,j,2)==0 && perp2(m,j,3)==0
else
for n=1:3
graphperp(:,:,n)=[Xs(1,j,n),(Xs(1,j,n)+perp2(m,j,n))];
end
plot3(graphperp(:,:,1),graphperp(:,:,2),graphperp(:,:,3),'g')
hold on
end
end
end
axis image
%% Radii Checks
radius=zeros(k(1),k(2)); %the magnitude or radius of a given perp2 vector at a given s
for m=1:k(2) %This for loop finds the magnitude, or radius, of each particular perpendicular
vector
for j=1:k(1)
radius(j,m)=sqrt((perp2(j,m,1)^2)+(perp2(j,m,2)^2)+(perp2(j,m,3)^2));
end
end
% radius_script
for m=1:k(2)
for j=1:k(1)
if radius(j,m)==0
radius(j,m)=NaN;
end
end
end
for j=1:k(2) %check of radius irregularities
for m=1:k(1)
if radius(m,j)>(min(radius(:,j))+5)
perp2(m,j,:)=0;
end
end

124

Appendix B.4 (cont’d)
end
clear radius
for m=1:k(2)
for j=1:k(1)
radius(j,m)=sqrt((perp2(j,m,1)^2)+(perp2(j,m,2)^2)+(perp2(j,m,3)^2));
end
end
meanR=zeros(1,length(s)); %the mean radius at a given s
k=size(perp2);
for m=1:length(s)%Finds the mean radius at each s
h=0;
for j=1:k(1)
if radius(j,m)~=0
h=h+1;
end
end
meanR(m)=sum(radius(:,m))/h;
end
maxR_all(i) = max(meanR(:));
figure;
for m=1:k(1)%Graphs the perpendicular perp2 vectors at their respective s
for j=1:k(2)
if perp2(m,j,1)==0 && perp2(m,j,2)==0 && perp2(m,j,3)==0
else
graphperp(:,:,1)=[Xs(1,j,1),(Xs(1,j,1)+perp2(m,j,1))];
graphperp(:,:,2)=[Xs(1,j,2),(Xs(1,j,2)+perp2(m,j,2))];
graphperp(:,:,3)=[Xs(1,j,3),(Xs(1,j,3)+perp2(m,j,3))];
plot3(graphperp(:,:,1),graphperp(:,:,2),graphperp(:,:,3),'g')
hold on
end
end
end
l=0;
for m=1:k(1)%Graphs the perpendicular perp2 vector endpoints at their respective s
for j=1:k(2)
125

Appendix B.4 (cont’d)
if perp2(m,j,1)==0 && perp2(m,j,2)==0 && perp2(m,j,3)==0
else
l=l+1;
raw_str(l,1,i)=[(Xs(1,j,1)+perp2(m,j,1))];
raw_str(l,2,i)=[(Xs(1,j,2)+perp2(m,j,2))];
raw_str(l,3,i)=[(Xs(1,j,3)+perp2(m,j,3))];
hold on
end
end
end
figure; subplot(1,2,1)
plot3(Xs(1,:,1),Xs(1,:,2),Xs(1,:,3),'k');
hold on;
scatter3(Xs(1,:,1),Xs(1,:,2),Xs(1,:,3),5,'k');
scatter3(A{i,1}(:,1),A{i,1}(:,2),A{i,1}(:,3),1,'CData',A{i,1}(:,3),'marker','.');
title('Raw (x,y,z) data');axis image;hold off;
subplot(1,2,2)
plot3(Xs(1,:,1),Xs(1,:,2),Xs(1,:,3),'k');
hold on;
scatter3(Xs(1,:,1),Xs(1,:,2),Xs(1,:,3),5,'k');
scatter3(raw_str(:,1,i),raw_str(:,2,i),raw_str(:,3,i),1,'CData',raw_str(:,3,i),'marker','.')
title('s Normal Plane (x,y,z) data');axis image; hold off;
clear N2_den N2_num dotprod graphN1 graphN2 graphN3 graphperp h j k l m n perp raw_str

126

Appendix B.5. step_4.m of the developed centerline generation algorithm
%% Find the magnitude of each individual perpendicular vector at a given s, and plot the mean
radius with respect to s.
theta=0:pi/20:(2*pi); %contains the vector with angle values for the parametric computation of
the circles
circ=zeros(1,length(theta),3); %contains information on the mean radius and orientation of the
plane at a particular s so that a circle can be plotted to approximate the surface.
figure;
for j=1:length(meanR) %Plots circles of uniform radii corresponding to the radius at the
particular perpendicular point.
for m=1:length(theta)
circ(1,m,1)=(meanR(j)*cos(theta(m))*N2(1,j,1))+(meanR(j)*sin(theta(m))*((N1(1,j,2)*N2(1,j,3)
)-(N1(1,j,3)*N2(1,j,2))))+Xs(1,j,1);
circ(1,m,2)=(meanR(j)*cos(theta(m))*N2(1,j,2))(meanR(j)*sin(theta(m))*((N1(1,j,1)*N2(1,j,3))-(N1(1,j,3)*N2(1,j,1))))+Xs(1,j,2);
circ(1,m,3)=(meanR(j)*cos(theta(m))*N2(1,j,3))+(meanR(j)*sin(theta(m))*((N1(1,j,1)*N2(1,j,2)
)-(N1(1,j,2)*N2(1,j,1))))+Xs(1,j,3);
end
plot3(circ(1,:,1)',circ(1,:,2)',circ(1,:,3)')
hold on
end
%Manipulates and labels the second figure properly.
axis image
grid on
title('Representation of Surface with Perpendicular Cross-sections A1','fontsize',12)
xlabel('X direction (mm)')
ylabel('Y direction (mm)')
zlabel('Z direction (mm)')
%Graphs the mean radius with respect to distance s along the centerline.
figure
plot(meanR,s)
%Manipulates and labels the third figure properly.
axis image
grid on
title('Mean radius with respect to s','fontsize',12)
xlabel('Mean radius (mm)')
ylabel('Longitudinal direction s (mm)')
clear j m
127

Appendix B.6. step_5.m of the developed centerline algorithm
perp2dotN=zeros(length(perp2),length(s),3); %contains dot product of perp2 with N2 and N3
depending on the 2nd or 3rd structured array in. The first structured array is empty
theta=zeros(length(perp2),length(s)); %contains the angle theta in radians that the perpendicular
vector makes with the N2 vector.
k=size(perp2);
for m=1:k(2)%a series of conditions and inverse tangent function formulas to determine which
angle theta should be assisgned to which perpendicular vector. The perpendicular plane was
divided up into four different regions that had their own conditions and formulas.
for j=1:k(1)
perp2norm =
[perp2(j,m,1)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]),perp2(j,m,2)/norm([perp2(j,m,1),p
erp2(j,m,2),perp2(j,m,3)]),perp2(j,m,3)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)])];
N2norm = [N2(1,m,1),N2(1,m,2),N2(1,m,3)];
N3norm = [N3(1,m,1),N3(1,m,2),N3(1,m,3)];
%
perp2dotN(j,m,2)=(N2(1,m,1))*(perp2(j,m,1)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]))+(
N2(1,m,2))*(perp2(j,m,2)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]))+(N2(1,m,3))*(perp2(j
,m,3)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]));
%
perp2dotN(j,m,3)=(N3(1,m,1))*(perp2(j,m,1)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]))+(
N3(1,m,2))*(perp2(j,m,2)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]))+(N3(1,m,3))*(perp2(j
,m,3)/norm([perp2(j,m,1),perp2(j,m,2),perp2(j,m,3)]));
perp2dotN(j,m,2)=dot(N2norm,perp2norm);
perp2dotN(j,m,3)=dot(N3norm,perp2norm);
perp2dotN(j,m,4)=dot(-N2norm,perp2norm);
perp2dotN(j,m,5)=dot(-N3norm,perp2norm);
if radius(j,m)~=0
if perp2dotN(j,m,2) == 1
theta(j,m) = 0;
elseif perp2dotN(j,m,3) == 1
theta(j,m) = pi/2;
elseif perp2dotN(j,m,4) == -1
theta(j,m) = pi
elseif perp2dotN(j,m,5) == -1
128

Appendix B.6 (cont’d)
theta(j,m) = 3*(pi/2);
elseif perp2dotN(j,m,2) > 0 && perp2dotN(j,m,3) > 0
theta(j,m) = acos(perp2dotN(j,m,2));
elseif perp2dotN(j,m,2) < 0 && perp2dotN(j,m,3) > 0
theta(j,m) = acos(perp2dotN(j,m,2));
elseif perp2dotN(j,m,2) > 0 && perp2dotN(j,m,3) < 0
theta(j,m) = -acos(perp2dotN(j,m,2)) + (2*pi);
elseif perp2dotN(j,m,2) < 0 && perp2dotN(j,m,3) < 0
theta(j,m) = -acos(perp2dotN(j,m,2)) + (2*pi);
else
end
else
end
end
end
%% Graph r(theta) along an s axis in 3D, fit an approximated surface over it
j=size(perp2);
snew=zeros(j(1),j(2)); %new array s that has an s value in each place holder that a radius of a
perp vector would appear.
for m=1:j(2)
snew(:,m)=s(m);
end
num_data=0;
% figure
for j=1:length(s) %Graph the radius wrt theta along an s axis to represent the 3D surface with a
straightened centerline and split about the N2 vector.
for m=1:length(perp2)
if theta(m,j)~=0
num_data=num_data+1;
%
plot3((theta(m,j)),snew(m,j),radius(m,j),'.r')
A{i,4}(num_data,1)=snew(m,j);
A{i,4}(num_data,2)=theta(m,j);
A{i,4}(num_data,3)=radius(m,j);
%
hold on
End
129

Appendix B.6 (cont’d)
end
end
% %Manipulates and labels the fourth figure properly.
% title('Raw data','fontsize',12)
% xlabel('Angle theta (rad)')
% ylabel('Longitudinal distance s (mm)')
% zlabel('Radius (mm)')
% set(gca,'XTick',0:pi/4:2*pi)
% grid on
%% Clearing variables
clear N2norm N3norm j k m perp2dotN perp2norm

130

Appendix B.7. step_6a of the developed approximation code to determine fitting coefficients
function [LSE] = step6a(beta)
global s num_data A i n_eta n_theta
n_node = (n_eta+1)*(n_theta);
node = zeros(n_node,2);
dy = s(end)/n_eta;
dx = (2*pi)/n_theta;
count = 0;
for m=1: n_eta+1
if m== (n_eta+1)
y_node = s(end);
else
y_node = (m-1)*dy;
end
for j=1:n_theta
x_node = (j-1)*dx;
count = count+1;
node(count, 1) = x_node;
node(count, 2) = y_node;
end
end
K = zeros(n_node, n_node);
F = zeros(n_node, 1);
for p=1:num_data
zeta = A{i,4}(p,1);
theta = A{i,4}(p,2);
for k = 1: n_node
x0_k=[node(k,2), node(k,1)];
ds = abs(zeta-x0_k(1));
dtheta = min( abs(theta-x0_k(2)), 2*pi-abs(theta-x0_k(2)));
phi_k = exp( -beta(1)*ds^2-beta(2)*dtheta^2);
F(k,1)=F(k,1)+phi_k*A{i,4}(p,3);
for n=1:n_node
x0_n=[node(n,2), node(n,1)];
ds = abs(zeta-x0_n(1));
dtheta = min( abs(theta-x0_n(2)), 2*pi-abs(theta-x0_n(2)));
131

Appendix B.7 (cont’d)
phi_n = exp( -beta(1)*ds^2-beta(2)*dtheta^2);
K(k,n)=K(k,n)+phi_k*phi_n;
end
end
end
alpha = K\F;
n_plot_theta =length(s);
n_plot_s = length(s);
s2=0;
theta=0;
data_predict = zeros(size(A{1,4},1), size(A{1,4},2));
count=0;
for m=1: length(A{i,4})
s2=A{i,4}(m,1);
theta=A{i,4}(m,2);
approx=0.0;
for k=1:n_node
x0=[node(k,2), node(k,1)];
ds2 = abs(s2-x0(1));
dtheta2 = min( abs(theta-x0(2)), 2*pi-abs(theta-x0(2)));
phi_k = exp( -beta(1)*ds2^2-beta(2)*dtheta2^2);
approx = approx+ alpha(k)*phi_k;
end
if approx < 0
approx = 0;
end
count=count+1;
data_predict(count, 1) = theta;
data_predict(count, 2) = s2;
data_predict(count, 3) = approx;
end
C{i,1}(:,1)=data_predict(:,2);
C{i,1}(:,2)=data_predict(:,1);
C{i,1}(:,3)=data_predict(:,3);
C{i,1}(1,4)=n_eta;
C{i,1}(1,5)=n_theta;
D{i,1} = zeros(size(A{1,4},1), size(A{1,4},2));
132

Appendix B.7 (cont’d)
D{i,1}(:,1) = A{i,4}(:,1);
D{i,1}(:,2) = A{i,4}(:,2);
D{i,1}(1,4)=n_eta;
D{i,1}(1,5)=n_theta;
for j = 1 : length(A{i,4})
D{i,1}(j,3) = C{i,1}(j,3) - A{i,4}(j,3);
end
abs_err = abs(mean(D{i,1}(:,3)));
D{i,1}(1,6)=abs_err;
LSE = sum(D{i,1}(:,3).^2);

133

Appendix B.8. step_6 of the developed approximation code
%%
global s num_data A i n_eta n_theta
n_eta = 10;
n_theta = 10;
options = optimset('Display','iter');
beta0 = [.9/100 1.5];
[beta, fval] = fminsearch(@step6a,beta0,options)
% eta_shapes = [1,2,3,4,5,6,7,8,9,10];
% theta_shapes = [1,2,3,4,5,6,7,8,9,10];
%
% for iii = 1:length(eta_shapes)
% for jjj = 1:length(theta_shapes)
%
global s num_data A i n_eta n_theta
%
n_eta = eta_shapes(iii);
%
n_theta = theta_shapes(jjj);
%
options = optimset('Display','iter');
%
beta0 = [.9/100 1.5];
%
[beta, fval] = fminsearch(@step6a,beta0,options)
%
kkk(iii,jjj)= fval
% end
% end

%% Surface Approximation
%Number of base functions.
C_count = 0;
for n_eta = n_eta
for n_theta = n_theta
C_count = C_count + 1;
n_node = (n_eta+1)*(n_theta);
node = zeros(n_node,2);
dy = s(end)/n_eta;
dx = (2*pi)/n_theta;
count = 0;
for m=1: n_eta+1
if m== (n_eta+1)
134

Appendix B.8 (cont’d)
y_node = s(end);
else
y_node = (m-1)*dy;
end
for j=1:n_theta
x_node = (j-1)*dx;
count = count+1;
node(count, 1) = x_node;
node(count, 2) = y_node;
end
end
K = zeros(n_node, n_node);
F = zeros(n_node, 1);
for p=1:num_data
zeta = A{i,4}(p,1);
theta = A{i,4}(p,2);
for k = 1: n_node
x0_k=[node(k,2), node(k,1)];
ds = abs(zeta-x0_k(1));
dtheta = min( abs(theta-x0_k(2)), 2*pi-abs(theta-x0_k(2)));
phi_k = exp( -beta(1)*ds^2-beta(2)*dtheta^2);
F(k,1)=F(k,1)+phi_k*A{i,4}(p,3);
for n=1:n_node
x0_n=[node(n,2), node(n,1)];
ds = abs(zeta-x0_n(1));
dtheta = min( abs(theta-x0_n(2)), 2*pi-abs(theta-x0_n(2)));
phi_n = exp( -beta(1)*ds^2-beta(2)*dtheta^2);
K(k,n)=K(k,n)+phi_k*phi_n;
end
end
end
alpha = K\F;
% Working method for approximated surface generation
% n_plot_theta =length(s);
135

Appendix B.8 (cont’d)
% n_plot_s = length(s);
% ds = s(end)/n_plot_s;
% dtheta = 2*pi/n_plot_theta;
% s2=0;
% theta=0;
% count = 0
% data_predict = zeros((n_plot_theta+1)*(n_plot_s+1), 3);
% for m=1: 1
% for j=1: n_plot_theta+1
%
s2= (m-1)*ds;
%
theta=(j-1)*dtheta;
%
approx=0.0;
%
for k=1:n_node
%
x0=[node(k,2), node(k,1)];
%
ds2 = abs(s2-x0(1));
%
dtheta2 = min( abs(theta-x0(2)), 2*pi-abs(theta-x0(2)));
%
%
phi_k = basis_cyl(s2, theta, x0, [beta1, beta2]);
%
phi_k = exp( -beta1*ds2^2-beta2*dtheta2^2);
%
approx = approx+ alpha(k)*phi_k;
%
end
%
if approx < 0
%
approx = 0;
%
end
%
count=count+1;
%
data_predict(count, 1) = theta;
%
data_predict(count, 2) = s2;
%
data_predict(count, 3) = approx;
% end
% end
% My method
n_plot_theta =length(s);
n_plot_s = length(s);
s2=0;
theta=0;
data_predict = zeros(size(A{1,4},1), size(A{1,4},2));
count=0;
for m=1: length(A{i,4})
s2=A{i,4}(m,1);
theta=A{i,4}(m,2);
approx=0.0;
136

Appendix B.8 (cont’d)
for k=1:n_node
x0=[node(k,2), node(k,1)];
ds2 = abs(s2-x0(1));
dtheta2 = min( abs(theta-x0(2)), 2*pi-abs(theta-x0(2)));
phi_k = exp( -beta(1)*ds2^2-beta(2)*dtheta2^2);
approx = approx+ alpha(k)*phi_k;
end
if approx < 0
approx = 0;
end
count=count+1;
data_predict(count, 1) = theta;
data_predict(count, 2) = s2;
data_predict(count, 3) = approx;
end
C{i,C_count}(:,1)=data_predict(:,2);
C{i,C_count}(:,2)=data_predict(:,1);
C{i,C_count}(:,3)=data_predict(:,3);
C{i,C_count}(1,4)=n_eta;
C{i,C_count}(1,5)=n_theta;
%Plots the known data and the approximated surface to see how good the fit is.
figure
plot3(A{i,4}(:,1), A{i,4}(:,2), A{i,4}(:,3),'r.','MarkerSize',4)
grid on
hold on;
plot3(C{i,C_count}(:,1), C{i,C_count}(:,2), C{i,C_count}(:,3),'b.','MarkerSize',4)
str_title = sprintf('Approximated (blue) vs. raw (red) for n eta = %d & n theta =
%d',n_eta,n_theta);
title(str_title,'fontsize',12,'fontname','Times')
xlabel('s','fontsize',12,'fontname','Times')
ylabel('\theta (deg)','fontsize',12,'fontname','Times')
zlabel('r (mm)','fontsize',12,'fontname','Times')
% Statistics of fit
D{i,C_count} = zeros(size(A{i,4},1), size(A{i,4},2));
D{i,C_count}(:,1) = A{i,4}(:,1);
D{i,C_count}(:,2) = A{i,4}(:,2);
D{i,C_count}(1,4)=n_eta;
D{i,C_count}(1,5)=n_theta;
137

Appendix B.8 (cont’d)
for j = 1 : length(A{i,4})
D{i,C_count}(j,3) = C{i,C_count}(j,3) - A{i,4}(j,3);
end
figure
plot3(D{i,C_count}(:,1), D{i,C_count}(:,2), D{i,C_count}(:,3),'g.','MarkerSize',4)
grid on
title('Residuals plot of approxmation fit','fontsize',12)
xlabel('s','fontsize',12,'fontname','Times')
ylabel('\theta (deg)','fontsize',12,'fontname','Times')
zlabel('r (mm)','fontsize',12,'fontname','Times')
abs_err = abs(mean(D{i,C_count}(:,3)))
D{i,C_count}(1,6)=abs_err;
clear n_node node dy dx count y_node y_node x_node K F beta1 beta2 hhh zeta theta x0_k
ds dtheta phi_k x0_n phi_n str_hhh alpha n_plot_theta n_plot_s s2 theta data_predict approx x0
ds2 dtheta2 phi_k
end
end
%% Surface Approximation take 2
% %Number of base functions.
%
% C_count = 0;
%
% for n_eta = n_eta
% for n_theta = n_theta
%
%
C_count = C_count + 1;
%
%
n_node = (n_eta+1)*(n_theta);
%
node = zeros(n_node,2);
%
dy = s(end)/n_eta;
%
dx = (2*pi)/n_theta;
%
%
count = 0;
%
for m=1: n_eta+1
%
if m== (n_eta+1)
%
y_node = s(end);
%
else
%
y_node = (m-1)*dy;
138

Appendix B.8 (cont’d)
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

end
for j=1:n_theta
x_node = (j-1)*dx;
count = count+1;
node(count, 1) = x_node;
node(count, 2) = y_node;
end
end
K = zeros(n_node, n_node);
F = zeros(n_node, 1);
for p=1:num_data
zeta = A{i,4}(p,1);
theta = A{i,4}(p,2);
for k = 1: n_node
x0_k=[node(k,2), node(k,1)];
ds = abs(zeta-x0_k(1));
dtheta = min( abs(theta-x0_k(2)), 2*pi-abs(theta-x0_k(2)));
phi_k = exp( -beta(1)*ds^2-beta(2)*dtheta^2);
F(k,1)=F(k,1)+phi_k*A{i,4}(p,3);
for n=1:n_node
x0_n=[node(n,2), node(n,1)];
ds = abs(zeta-x0_n(1));
dtheta = min( abs(theta-x0_n(2)), 2*pi-abs(theta-x0_n(2)));
phi_n = exp( -beta(1)*ds^2-beta(2)*dtheta^2);
K(k,n)=K(k,n)+phi_k*phi_n;
end
end
end
alpha = K\F;
% Working method for approximated surface generation
% n_plot_theta =length(s);
% n_plot_s = length(s);
% ds = s(end)/n_plot_s;
% dtheta = 2*pi/n_plot_theta;
139

Appendix B.8 (cont’d)
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

% s2=0;
% theta=0;
% count = 0
% data_predict = zeros((n_plot_theta+1)*(n_plot_s+1), 3);
% for m=1: 1
% for j=1: n_plot_theta+1
%
s2= (m-1)*ds;
%
theta=(j-1)*dtheta;
%
approx=0.0;
%
for k=1:n_node
%
x0=[node(k,2), node(k,1)];
%
ds2 = abs(s2-x0(1));
%
dtheta2 = min( abs(theta-x0(2)), 2*pi-abs(theta-x0(2)));
%
phi_k = exp( -beta1*ds2^2-beta2*dtheta2^2);
%
approx = approx+ alpha(k)*phi_k;
%
end
%
if approx < 0
%
approx = 0;
%
end
%
count=count+1;
%
data_predict(count, 1) = theta;
%
data_predict(count, 2) = s2;
%
data_predict(count, 3) = approx;
% end
% end
% My method
n_plot_theta =length(s);
n_plot_s = length(s);
ds = 1;
theta_div = 200; %Some even number, div per x-sec from 0 - 2*pi
dtheta = (2*pi) / theta_div;
data_predict = zeros(length(s)*(theta_div), length(s)*(theta_div));
count=0;
for j = 1 : length(s)
s2 = j-1;
for m = 1 : theta_div
theta = dtheta*m;
approx=0.0;
for k=1:n_node
x0=[node(k,2), node(k,1)];
140

Appendix B.8 (cont’d)
%
ds2 = abs(s2-x0(1));
%
dtheta2 = min( abs(theta-x0(2)), 2*pi-abs(theta-x0(2)));
%
phi_k = exp( -beta(1)*ds2^2-beta(2)*dtheta2^2);
%
approx = approx+ alpha(k)*phi_k;
%
end
%
if approx < 0
%
approx = 0;
%
end
%
count=count+1;
%
data_predict(count, 1) = theta;
%
data_predict(count, 2) = s2;
%
data_predict(count, 3) = approx;
%
end
%
end
%
%
C{i,C_count}(:,1)=data_predict(:,2);
%
C{i,C_count}(:,2)=data_predict(:,1);
%
C{i,C_count}(:,3)=data_predict(:,3);
%
C{i,C_count}(1,4)=n_eta;
%
C{i,C_count}(1,5)=n_theta;
%
%
%Plots the known data and the approximated surface to see how good the fit is.
%
%
figure
%
plot3(A{i,4}(:,1), A{i,4}(:,2), A{i,4}(:,3),'r.','MarkerSize',4)
%
grid on
%
hold on;
%
plot3(C{i,C_count}(:,1), C{i,C_count}(:,2), C{i,C_count}(:,3),'b.','MarkerSize',4)
%
str_title = sprintf('Approximated (blue) vs. raw (red) for n eta = %d & n theta =
%d',n_eta,n_theta);
%
title(str_title,'fontsize',12,'fontname','Times')
%
xlabel('s (mm)','fontsize',12,'fontname','Times')
%
ylabel('\theta (deg)','fontsize',12,'fontname','Times')
%
zlabel('r (mm)','fontsize',12,'fontname','Times')
%
%
% Statistics of fit
%
%
clear n_node node dy dx count y_node y_node x_node K F beta1 beta2 hhh zeta theta
x0_k ds dtheta phi_k x0_n phi_n str_hhh alpha n_plot_theta n_plot_s s2 theta data_predict
approx x0 ds2 dtheta2 phi_k
%
% end
% end
clear abs_err j k m n n_eta n_theta p
141

Appendix B.9. ecc_testing.m code developed for automatic eccentricity testing
sc = 0;
k = 0;
for j = 1:length(A{i,4})
if j == length(A{i,4})
k = k+1;
x_img(k,1) = A{i,4}(j,2);
x_img(k,2) = A{i,4}(j,3);
x_img(:,:) = sortrows(x_img(:,:),1);
x_img(k+1,:) = x_img(1,:);
kk = 0;
for l = 1:length(x_img)
for m = 1:length(x_img)
kk = kk + 1;
x_img_1(kk) = x_img(l,1);
x_img_2(kk) = x_img(l,2);
x_img_3(kk) = x_img(m,1);
x_img_4(kk) = x_img(m,2);
R_d(kk) = sqrt((x_img_2(kk)^2) + (x_img_4(kk)^2) (2*x_img_2(kk)*x_img_4(kk)*cos(x_img_1(kk)-x_img_3(kk))));
end
end
[R_max(sc+1),I_R] = max(R_d(:));
x_img_5 = x_img_1(I_R) + (pi/2);
x_img_6 = x_img_1(I_R) - (pi/2);
if x_img_5 > (2*pi)
x_img_5 = x_img_5 - (2*pi);
elseif x_img_5 < 0
x_img_5 = x_img_5 + (2*pi);
end
if x_img_6 > (2*pi)
x_img_6 = x_img_6 - (2*pi);
elseif x_img_6 < 0
142

Appendix B.9 (cont’d)
x_img_6 = x_img_6 + (2*pi);
end
for l = 1:length(x_img)
t_1(l) = (abs(x_img_5 - x_img(l))) / x_img(l);
t_2(l) = (abs(x_img_6 - x_img(l))) / x_img(l);
end
[t_1_val,t_1_ind] = min(t_1(:));
[t_2_val,t_2_ind] = min(t_2(:));
R_min(sc+1) = sqrt((x_img(t_1_ind,2)^2) + (x_img(t_2_ind,2)^2) (2*x_img(t_1_ind,2)*x_img(t_2_ind,2)*cos(x_img(t_1_ind,1)-x_img(t_2_ind,1))));
figure('Color',[1 1 1]);
polar(x_img(:,1),x_img(:,2),'b');
axis image;
hold on;
polar([x_img_1(I_R) x_img_3(I_R)],[x_img_2(I_R) x_img_4(I_R)],'r')
polar([x_img(t_1_ind,1) x_img(t_2_ind,1)],[x_img(t_1_ind,2) x_img(t_2_ind,2)],'g')
clear x_img x_img_1 x_img_2 x_img_3 x_img_4 x_img_5 x_img_6 I_R R_d t_1 t_2
t_1_val t_2_val t_1_ind t_2_ind
elseif A{i,4}(j,1) == sc
k = k+1;
x_img(k,1) = A{i,4}(j,2);
x_img(k,2) = A{i,4}(j,3);
elseif A{i,4}(j,1) ~= sc
x_img(:,:) = sortrows(x_img(:,:),1);
x_img(k+1,:) = x_img(1,:);
kk = 0;
for l = 1:length(x_img)
for m = 1:length(x_img)
kk = kk + 1;
x_img_1(kk) = x_img(l,1);
x_img_2(kk) = x_img(l,2);
x_img_3(kk) = x_img(m,1);
143

Appendix B.9 (cont’d)
x_img_4(kk) = x_img(m,2);
R_d(kk) = sqrt((x_img_2(kk)^2) + (x_img_4(kk)^2) (2*x_img_2(kk)*x_img_4(kk)*cos(x_img_1(kk)-x_img_3(kk))));
end
end
[R_max(sc+1),I_R] = max(R_d(:));
x_img_5 = x_img_1(I_R) + (pi/2);
x_img_6 = x_img_1(I_R) - (pi/2);
if x_img_5 > (2*pi)
x_img_5 = x_img_5 - (2*pi);
elseif x_img_5 < 0
x_img_5 = x_img_5 + (2*pi);
end
if x_img_6 > (2*pi)
x_img_6 = x_img_6 - (2*pi);
elseif x_img_6 < 0
x_img_6 = x_img_6 + (2*pi);
end
for l = 1:length(x_img)
t_1(l) = (abs(x_img_5 - x_img(l))) / x_img(l);
t_2(l) = (abs(x_img_6 - x_img(l))) / x_img(l);
end
[t_1_val,t_1_ind] = min(t_1(:));
[t_2_val,t_2_ind] = min(t_2(:));
R_min(sc+1) = sqrt((x_img(t_1_ind,2)^2) + (x_img(t_2_ind,2)^2) (2*x_img(t_1_ind,2)*x_img(t_2_ind,2)*cos(x_img(t_1_ind,1)-x_img(t_2_ind,1))));
%
figure('Color',[1 1 1]);
%
polar(x_img(:,1),x_img(:,2),'b');
%
axis image;
%
hold on;
%
polar([x_img_1(I_R) x_img_3(I_R)],[x_img_2(I_R) x_img_4(I_R)],'r')
%
polar([x_img(t_1_ind,1) x_img(t_2_ind,1)],[x_img(t_1_ind,2)
x_img(t_2_ind,2)],'g')

144

Appendix B.9 (cont’d)
clear x_img x_img_1 x_img_2 x_img_3 x_img_4 x_img_5 x_img_6 I_R R_d t_1 t_2
t_1_val t_2_val t_1_ind t_2_ind
k=1;
sc = sc+1;
x_img(k,1) = A{i,4}(j,2);
x_img(k,2) = A{i,4}(j,3);
end
end
ecc = R_max ./ R_min;
sc = 0;
for j = 1:length(A{i,4})
if A{i,4}(j,1) == sc
A{i,4}(j,6) = ecc(sc+1);
elseif A{i,4}(j,1) ~= sc
sc = sc+1;
A{i,4}(j,6) = ecc(sc+1);
end
end
clear k kk sc l m R_d

145

Appendix B.10. surf_curv.m code developed for automatic curvature determination
%% post fitting curvature
sc = -1;
count = 0;
for j = 1 : length(s)
sc = sc+1;
for k = 1 : theta_div
count = count+1;
if k == 1
R1 = C{i,C_count}((sc*theta_div)+theta_div,3);
R2 = C{i,C_count}((sc*theta_div)+k,3);
R3 = C{i,C_count}((sc*theta_div)+k+1,3);
t1 = C{i,C_count}((sc*theta_div)+theta_div,2);
t2 = C{i,C_count}((sc*theta_div)+k,2);
t3 = C{i,C_count}((sc*theta_div)+k+1,2);
elseif k == theta_div
R1 = C{i,C_count}((sc*theta_div)+k-1,3);
R2 = C{i,C_count}((sc*theta_div)+k,3);
R3 = C{i,C_count}((sc*theta_div)+1,3);
t1 = C{i,C_count}((sc*theta_div)+k-1,2);
t2 = C{i,C_count}((sc*theta_div)+k,2);
t3 = C{i,C_count}((sc*theta_div)+1,2);
else
R1 = C{i,C_count}((sc*theta_div)+k-1,3);
R2 = C{i,C_count}((sc*theta_div)+k,3);
R3 = C{i,C_count}((sc*theta_div)+k+1,3);
t1 = C{i,C_count}((sc*theta_div)+k-1,2);
t2 = C{i,C_count}((sc*theta_div)+k,2);
t3 = C{i,C_count}((sc*theta_div)+k+1,2);
end
x1 = R1*cos(t1);
x2 = R2*cos(t2);
146

Appendix B.10 (cont’d)
x3 = R3*cos(t3);
y1 = R1*sin(t1);
y2 = R2*sin(t2);
y3 = R3*sin(t3);
m1 = (y2-y1)/(x2-x1);
m2 = (y3-y2)/(x3-x2);
xc = (((m1*m2)*(y1-y3)) + (m2*(x1+x2)) - (m1*(x2+x3))) / (2*(m2-m1));
yc = ((-1/m1)*(xc-((x1+x2)/2))) + ((y1+y2)/2);
R = sqrt(((x2-xc)^2) + ((y2-yc)^2));
C{i,C_count}(count,7) = 1/R;
clear R1 R2 R3 t1 t2 t3 x1 x2 x3 y1 y2 y3 m1 m2 xc yc R
end
end

147

APPENDIX C

Results of centerline and surface approximation algorithm performance

148

Figure C.1. A 3-D model of the lumen for D3 as seen in the x-z plane (upper-left), a coronal, or
x-z, slice from scan D3 (upper-right), an x-z (bottom-left), and a y-z (bottom-right) view of the
centerline algorithm results for scan D3. It is clear that scan D3 features a geometric anomaly.
This anomaly is highly unusual and tends to force the centerline generation to follow the
observed path regardless of input parameters.
149

Figure C.2. An x-z (left) and y-z (right) view of the centerline algorithm results from scan F5. If
the final approximation of a centerline point differs greatly from the initial guess, the vector
between the previous centerline point and current centerline point becomes more extreme. In the
case when the vector points directly at a vessel boundary the algorithm will tend to fail. This
failure occurs because the minimum distance vector is nearly parallel to the axis normal to the
projection plane used for translations, resulting in near-zero translational adjustment vectors.
The algorithm assumes that at any given two points the vector between them is relatively parallel
to the centerline of the vessel.

150

Figure C.3. An x-z (left) and y-z (right) view of the centerline algorithm results from scan H5. If
the final approximation of a centerline point differs greatly from the initial guess, the vector
between the previous centerline point and current centerline point becomes more extreme. In the
case when the vector points directly at a vessel boundary the algorithm will tend to fail. This
failure occurs because the minimum distance vector is nearly parallel to the axis normal to the
projection plane used for translations, resulting in near-zero translational adjustment vectors.
The algorithm assumes that at any given two points the vector between them is relatively parallel
to the centerline of the vessel. For this case the algorithm was able to recover, but it resulted in a
large erroneous region before recovery.

Figure C.4. An x-z (left) and y-z (right) view of the centerline algorithm results from scan K3. It
is unclear why the centerline fluctuated from its expected center, but it is clearly erroneous and
requires adjustment of operating parameters.

151

Table C.1. SSE values from surface approximations sorted first by Nη and then by Nθ
combo #
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

Nη

Nθ
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
3
3
4
4
4
4
4
4
4
4
4

1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9

152

SSE
56956.98
56108.11
56312.28
55796.18
55296.88
55137.28
54892.25
54792.57
54611
54499.47
17815
17125.56
16736.05
16336.73
15430.52
15452.52
15157.9
15176.29
15049
15038.77
11652.39
10549.91
10432.27
9682.665
8776.14
8765.256
8547.275
8489.009
8269.801
8283.601
6487.158
6042.208
5676.069
5211.81
4139.079
4108.915
3904.917
3854.168
3747.15

Table C.1 (cont’d)
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

4
5
5
5
5
5
5
5
5
5
5
6
6
6
6
6
6
6
6
6
6
7
7
7
7
7
7
7
7
7
7
8
8
8
8
8
8
8
8
8
8

10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10

153

3703.919
5955.562
5603.217
5077.259
4528.495
3557.47
3509.808
3253.571
3216.315
3089.596
3044.364
5932.707
5311.066
5100.357
4378.236
3290.206
3269.365
2946.388
2868.631
2752.81
2713.908
6082.858
5466.559
5097.665
4436.855
3181.129
3249.602
2767.998
2699.992
2554.182
2480.305
5563.773
4903.801
4671.616
3896.41
2684.641
2635.93
2306.582
2204.442
2035.634
1974.096

Table C.1 (cont’d)
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

9
9
9
9
9
9
9
9
9
9
10
10
10
10
10
10
10
10
10
10

1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10

154

5612.163
5180.029
4673.011
3833.785
2583.552
2499.712
2132.18
1986.855
1846.959
1729.718
5579.099
5131.247
4651.369
3779.547
2608.741
2532.217
2078.223
1951.817
1823.693
1673.866

Table C.2. SSE approximations sorted based on SSE values
combo #
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

Nη
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
3
3
4
7
4
5
6
4
9
5
10
8

Nθ
1
3
2
4
5
6
7
8
9
10
1
2
3
4
6
5
8
7
9
10
1
2
3
4
5
6
7
8
10
9
1
1
2
1
1
3
1
2
1
1

155

SSE
56956.98
56312.28
56108.11
55796.18
55296.88
55137.28
54892.25
54792.57
54611
54499.47
17815
17125.56
16736.05
16336.73
15452.52
15430.52
15176.29
15157.9
15049
15038.77
11652.39
10549.91
10432.27
9682.665
8776.14
8765.256
8547.275
8489.009
8283.601
8269.801
6487.158
6082.858
6042.208
5955.562
5932.707
5676.069
5612.163
5603.217
5579.099
5563.773

Table C.2 (cont’d)
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

7
6
4
9
10
6
7
5
8
9
8
10
5
7
6
4
4
4
8
4
9
10
4
4
5
5
6
6
5
7
5
7
5
5
6
6
7
6
6
7
8

156

2
2
4
2
2
3
3
3
2
3
3
3
4
4
4
5
6
7
4
8
4
4
9
10
5
6
5
6
7
6
8
5
9
10
7
8
7
9
10
8
5

5466.559
5311.066
5211.81
5180.029
5131.247
5100.357
5097.665
5077.259
4903.801
4673.011
4671.616
4651.369
4528.495
4436.855
4378.236
4139.079
4108.915
3904.917
3896.41
3854.168
3833.785
3779.547
3747.15
3703.919
3557.47
3509.808
3290.206
3269.365
3253.571
3249.602
3216.315
3181.129
3089.596
3044.364
2946.388
2868.631
2767.998
2752.81
2713.908
2699.992
2684.641

Table C.2 (cont’d)
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

8
10
9
7
10
9
7
8
8
9
10
8
9
8
10
9
10
9
10

157

6
5
5
9
6
6
10
7
8
7
7
9
8
10
8
9
9
10
10

2635.93
2608.741
2583.552
2554.182
2532.217
2499.712
2480.305
2306.582
2204.442
2132.18
2078.223
2035.634
1986.855
1974.096
1951.817
1846.959
1823.693
1729.718
1673.866

REFERENCES

158

REFERENCES

[1]

C. Porth, Pathophysiology : co cepts of altered health states, 7th ed. Lippincott Williams
& Wilkins, 2005.

[2]

M. Jessup and S. Brozena, “Heart Failure,” New England Journal of Medicine, vol. 348,
no. 20, pp. 2007–2018, 2003.

[3]

R. W. Schrier and W. T. Abraham, “Hormones and Hemodynamics in Heart Failure,” New
England Journal of Medicine, vol. 341, no. 8, pp. 577–585, 1999.

[4]

W. D. Clouse, J. W. Hallett Jr, H. V. Schaff, M. M. Gayari, D. M. Ilstrup, and L. J. Melton
III, “Improved prognosis of thoracic aortic aneurysms,” JAMA: The Journal of the
American Medical Association, vol. 280, no. 22, pp. 1926–1929, 1998.

[5]

H. W. Kniemeyer, T. Kessler, P. U. Reber, H. B. Ris, H. Hakki, and M. K. Widmer,
“Treatment of Ruptured Abdominal Aortic Aneurysm, a Permanent Challenge or a Waste
of Resources? Prediction of Outcome Using a Multi-organ-dysfunction Score,” European
Journal of Vascular and Endovascular Surgery, vol. 19, no. 2, pp. 190–196, Feb. 2000.

[6]

T. Mussivand, “Mechanical Circulatory Devices for the Treatment of Heart Failure,”
Journal of Cardiac Surgery, vol. 14, no. 3, pp. 218–228, 1999.

[7]

F. A. Lederle, S. E. Wilson, G. R. Johnson, D. B. Reinke, F. N. Littooy, C. W. Acher, D. J.
Ballard, L. M. Messina, I. L. Gordon, E. P. Chute, and others, “Immediate repair compared
with surveillance of small abdominal aortic aneurysms,” New England Journal of Medicine,
vol. 346, no. 19, pp. 1437–1444, 2002.

[8]

T. F. US Preventive Services, “Screening for Abdominal Aortic Aneurysm:
Recommendation Statement,” Annals of Internal Medicine, vol. 142, no. 3, pp. 198–202,
Feb. 2005.

[9]

J. L. Cronenwett, “Variables that affect the expansion rate and rupture of abdominal aortic
aneurysms,” Annals of the New York Academy of Sciences, vol. 800, no. 1, pp. 56–67, 1996.

[10] K. Halliday and A. Al-Kutoubi, “Draped Aorta: CT Sign of Contained Leak of Aortic
Aneurysms,” Radiology, vol. 199, no. 1, pp. 41–43, 1996.
[11] C. Dubost, M. Allary, and N. Oeconomos, “Resection of an aneurysm of the abdominal
aorta: reestablishment of the continuity by a preserved human arterial graft, with result after
five months,” Archives of Surgery, vol. 64, no. 3, p. 405, 1952.

159

[12] J. C. Parodi, J. C. Palmaz, and H. D. Barone, “Transfemoral intraluminal graft implantation
for abdominal aortic aneurysms,” Annals of vascular surgery, vol. 5, no. 6, pp. 491–499,
1991.
[13] G. Martufi, E. S. Di Martino, C. H. Amon, S. C. Muluk, and E. A. Finol, “ThreeDimensional Geometrical Characterization of Abdominal Aortic Aneurysms: Image-Based
Wall Thickness Distribution,” J. Biomech. Eng., vol. 131, no. 6, p. 061015, 2009.
[14] R. Limet, N. Sakalihassan, and A. Albert, “Determination of the expansion rate and
incidence of rupture of abdominal aortic aneurysms,” J Vasc Surg, vol. 14, no. 4, pp. 540–
548, Oct. 1991.
[15] K. Myers, T. Devine, C. Barras, and G. Self, “Endoluminal Versus Open repair for
abdominal aortic aneurysms,” in 2nd Virtual Congress of Cardiology, 2001.
[16] D. Drury, J. A. Michaels, L. Jones, and L. Ayiku, “Systematic review of recent evidence
for the safety and efficacy of elective endovascular repair in the management of infrarenal
abdominal aortic aneurysm,” British Journal of Surgery, vol. 92, no. 8, pp. 937–946, Aug.
2005.
[17] “Long-Term Outcomes of Immediate Repair Compared with Surveillance of Small
Abdominal Aortic Aneurysms,” New England Journal of Medicine, vol. 346, no. 19, pp.
1445–1452, 2002.
[18] J. T. Powell, “Final 12-year follow-up of Surgeryversus Surveillance in the UK Small
Aneurysm Trial,” British Journal of Surgery, vol. 94, no. 6, pp. 702–708, Jun. 2007.
[19] C. D. Karkos, U. Mukhopadhyay, I. Papakostas, J. Ghosh, G. J. L. Thomson, and R.
Hughes, “Abdominal Aortic Aneurysm: the Role of Clinical Examination and
Opportunistic Detection,” European Journal of Vascular and Endovascular Surgery, vol.
19, no. 3, pp. 299 – 303, 2000.
[20] D. H. J. Wang, M. S. Makaroun, M. W. Webster, and D. A. Vorp, “Effect of intraluminal
thrombus on wall stress in patient-specific models of abdominal aortic aneurysm,” Journal
of Vascular Surgery, vol. 36, no. 3, pp. 598–604, Sep. 2002.
[21] U. S. A. T. Participants, “Long-term outcomes of immediate repair compared with
surveillance of small abdominal aortic aneurysms,” N. Engl. J. Med, vol. 346, pp. 1445–
1452, 2002.
[22] R. C. Darling, C. R. Messina, D. C. Brewster, and L. W. Ottinger, “Autopsy study of
unoperated abdominal aortic aneurysms. The case for early resection.,” Circulation, vol. 56,
no. 3 Suppl, 1977.
[23] G. Giannoglou, G. Giannakoulas, J. Soulis, Y. Chatzizisis, T. Perdikides, N. Melas, G.
Parcharidis, and G. Louridas, “Predicting the Risk of Rupture of Abdominal Aortic
160

Aneurysms by Utilizing Various Geometrical Parameters: Revisiting the Diameter
Criterion,” Angiology, vol. 57, no. 4, pp. 487–494, Aug. 2006.
[24] M. L. Raghavan and D. A. Vorp, “Toward a biomechanical tool to evaluate rupture
potential of abdominal aortic aneurysm: identification of a finite strain constitutive model
and evaluation of its applicability,” Journal of biomechanics, vol. 33, no. 4, pp. 475–482,
2000.
[25] B. J. Doyle, A. Callanan, P. E. Burke, P. A. Grace, M. T. Walsh, D. A. Vorp, and T. M.
McGloughlin, “Vessel asymmetry as an additional diagnostic tool in the assessment of
abdominal aortic aneurysms,” Journal of Vascular Surgery, vol. 49, no. 2, pp. 443–454,
2009.
[26] M. L. Raghavan, D. A. Vorp, M. P. Federle, M. S. Makaroun, and M. W. Webster, “Wall
stress distribution on three-dimensionally reconstructed models of human abdominal aortic
aneurysm,” Journal of Vascular Surgery, vol. 31, no. 4, pp. 760–769, 2000.
[27] E. Georgakarakos, C. V. Ioannou, Y. Kamarianakis, Y. Papaharilaou, T. Kostas, E.
Manousaki, and A. N. Katsamouris, “The role of geometric parameters in the prediction of
abdominal aortic aneurysm wall stress,” European Journal of Vascular and Endovascular
Surgery, vol. 39, no. 1, pp. 42–48, 2010.
[28] L. Speelman, A. Bohra, E. M. H. Bosboom, G. W. H. Schurink, F. N. van de Vosse, M. S.
Makaroun, and D. A. Vorp, “Effects of wall calcifications in patient-specific wall stress
analyses of abdominal aortic aneurysms,” Journal of Biomechanical Engineering, vol. 129,
pp. 105–109, 2007.
[29] J. Shum, G. Martufi, E. Martino, C. B. Washington, J. Grisafi, S. C. Muluk, and E. A.
Finol, “Quantitative Assessment of Abdominal Aortic Aneurysm Geometry,” Ann Biomed
Eng, vol. 39, no. 1, pp. 277–286, Oct. 2010.
[30] M. Xenos, S. Rambhia, Y. Alemu, S. Einav, J. J. Ricotta, N. Labropoulos, A. Tassiopoulos,
and D. Bluestein, “Mimics Based Image Reconstruction Augments Diagnosis and
Management of Vascular Pathologies: A Study of Ruptured Abdominal Aortic Aneurisms.”
[31] D. Bluestein, K. Dumont, M. De Beule, J. Ricotta, P. Impellizzeri, B. Verhegghe, and P.
Verdonck, “Intraluminal thrombus and risk of rupture in patient specific abdominal aortic
aneurysm - FSI modelling,” Comp. Methods in Biomechanics & Biomedical Eng., vol. 12,
no. 1, pp. 73–81, Feb. 2009.
[32] M. F. Fillinger, M. L. Raghavan, S. P. Marra, J. L. Cronenwett, and F. E. Kennedy, “In
vivo analysis of mechanical wall stress and abdominal aortic aneurysm rupture risk,” J Vasc
Surg, vol. 36, no. 3, pp. 589–597, 2002.

161

[33] B. Wolters, M. C. M. Rutten, G. W. H. Schurink, U. Kose, J. De Hart, and F. N. Van De
Vosse, “A patient-specific computational model of fluid–structure interaction in abdominal
aortic aneurysms,” Medical engineering & physics, vol. 27, no. 10, pp. 871–883, 2005.
[34] J. H. Leung, A. R. Wright, N. Cheshire, J. Crane, S. A. Thom, A. D. Hughes, and Y. Xu,
“Fluid structure interaction of patient specific abdominal aortic aneurysms: a comparison
with solid stress models,” BioMedical Engineering OnLine, vol. 5, no. 1, p. 33, 2006.
[35] Y. Papaharilaou, J. A. Ekaterinaris, E. Manousaki, and A. N. Katsamouris, “A decoupled
fluid structure approach for estimating wall stress in abdominal aortic aneurysms,” Journal
of biomechanics, vol. 40, no. 2, pp. 367–377, 2007.
[36] A. Borghi, N. B. Wood, R. H. Mohiaddin, and X. Y. Xu, “3D geometric reconstruction of
thoracic aortic aneurysms,” Biomedical engineering online, vol. 5, no. 1, p. 59, 2006.
[37] D. A. Steinman, J. S. Milner, C. J. Norley, S. P. Lownie, and D. W. Holdsworth, “Imagebased computational simulation of flow dynamics in a giant intracranial aneurysm,”
American journal of neuroradiology, vol. 24, no. 4, pp. 559–566, 2003.
[38] J. P. V. Geest, D. H. Wang, S. R. Wisniewski, M. S. Makaroun, and D. A. Vorp, “Towards
a noninvasive method for determination of patient-specific wall strength distribution in
abdominal aortic aneurysms,” Annals of Biomedical Engineering, vol. 34, pp. 1098–1106,
2006.
[39] D. A. Vorp and J. P. Geest, “Biomechanical determinants of abdominal aortic aneurysm
rupture,” Arterioscler. Thromb. Vasc. Biol., vol. 25, pp. 1558–1566, 2005.
[40] B. J. Doyle, T. J. Corbett, A. Callanan, M. T. Walsh, D. A. Vorp, and T. M. McGloughlin,
“An experimental and numerical comparison of the rupture locations of an abdominal aortic
aneurysm,” Journal Information, vol. 16, no. 3, 2009.
[41] J. Swedenborg and P. Eriksson, “The Intraluminal Thrombus as a Source of Proteolytic
Activity,” Annals of the New York Academy of Sciences, vol. 1085, no. 1, pp. 133–138,
Nov. 2006.
[42] Y. G. Wolf, W. S. Thomas, F. J. Brennan, W. G. Goff, M. J. Sise, and E. F. Bernstein,
“Computed tomography scanning findings associated with rapid expansion of abdominal
aortic aneurysms,” J Vasc Surg, vol. 20, no. 4, pp. 529–538, Oct. 1994.
[43] J. W. Hinnen, O. H. . Koning, M. J. . Visser, and H. J. Van Bockel, “Effect of intraluminal
thrombus on pressure transmission in the abdominal aortic aneurysm,” Journal of vascular
surgery, vol. 42, no. 6, pp. 1176–1182, 2005.
[44] D. A. Vorp, M. L. Raghavan, and M. W. Webster, “Mechanical wall stress in abdominal
aortic aneurysm: influence of diameter and asymmetry,” Journal of Vascular Surgery, vol.
27, no. 4, pp. 632–639, 1998.
162

[45] P. N. Watton, N. A. Hill, and M. Heil, “A mathematical model for the growth of the
abdominal aortic aneurysm,” Biomechanics and Modeling in Mechanobiology, vol. 3, no. 2,
pp. 98–113, Sep. 2004.
[46] C. Kleinstreuer and Z. Li, “Analysis and computer program for rupture-risk prediction of
abdominal aortic aneurysms,” Biomed Eng Online, vol. 5, p. 19, 2006.
[47] M. S. Sacks, D. A. Vorp, M. L. Raghavan, M. P. Federle, and M. W. Webster, “In vivo
three-dimensional surface geometry of abdominal aortic aneurysms,” Annals of biomedical
engineering, vol. 27, no. 4, pp. 469–479, 1999.
[48] J. Shum, A. Xu, I. Chatnuntawech, and E. A. Finol, “A Framework for the Automatic
Generation of Surface Topologies for Abdominal Aortic Aneurysm Models,” Annals of
Biomedical Engineering, vol. 39, no. 1, pp. 249–259, Sep. 2010.
[49] G. Choi, C. P. Cheng, N. M. Wilson, and C. A. Taylor, “Methods for Quantifying ThreeDimensional Deformation of Arteries due to Pulsatile and Nonpulsatile Forces:
Implications for the Design of Stents and Stent Grafts,” Annals of Biomedical Engineering,
vol. 37, no. 1, pp. 14–33, Nov. 2008.
[50] C. B. Washington, J. Shum, S. C. Muluk, and E. A. Finol, “The Association of Wall
Mechanics and Morphology: A Case Study of Abdominal Aortic Aneurysm Growth,”
Journal of biomechanical engineering, vol. 133, p. 104501, 2011.
[51] A. Long, L. Rouet, J. S. Lindholt, and E. Allaire, “Measuring the Maximum Diameter of
Native Abdominal Aortic Aneurysms: Review and Critical Analysis,” European Journal of
Vascular and Endovascular Surgery, Feb. 2012.
[52] B. S. Ferket, N. Grootenboer, E. B. Colkesen, J. J. Visser, M. R. H. M. van Sambeek, S.
Spronk, E. W. Steyerberg, and M. G. M. Hunink, “Systematic review of guidelines on
abdominal aortic aneurysm screening,” J Vasc Surg, Feb. 2011.
[53] J. T. Powell, M. J. Sweeting, L. C. Brown, S. M. Gotensparre, F. G. Fowkes, and S. G.
Thompson, “Systematic review and meta-analysis of growth rates of small abdominal aortic
aneurysms,” British Journal of Surgery, vol. 98, no. 5, pp. 609–618, May 2011.
[54] R. Acharya, R. Wasserman, J. Stevens, and C. Hinojosa, “Biomedical imaging modalities:
a tutorial,” Computerized Medical Imaging and Graphics, vol. 19, no. 1, pp. 3–25, 1995.
[55] Olaf T. von Ramm and Stephen W. Smith, “Three-dimensional Imaging System,” U.S.
Patent 469443415-Sep-1987.
[56] J. A. Ten Bosch, E. V. Rouwet, C. T. . Peters, L. Jansen, H. J. . Verhagen, M. H. Prins, and
J. A. . Teijink, “Contrast-enhanced ultrasound versus computed tomographic angiography

163

for surveillance of endovascular abdominal aortic aneurysm repair,” Journal of Vascular
and Interventional Radiology, vol. 21, no. 5, pp. 638–643, 2010.
[57] N. Grøndal, M. B. Bramsen, M. D. Thomsen, C. B. Rasmussen, and J. S. Lindholt, “The
Cardiac Cycle is a Major Contributor to Variability in Size Measurements of Abdominal
Aortic Aneurysms by Ultrasound,” European Journal of Vascular and Endovascular
Surgery, 2011.
[58] A. Abbas, A. Smith, M. Cecelja, and M. Waltham, “Assessment of the accuracy of
AortaScan for detection of abdominal aortic aneurysm (AAA),” European Journal of
Vascular and Endovascular Surgery, 2011.
[59] Mimics Medical Imaging Software. Leuven, Belgium: Materialise, Inc.
[60] S. Zeinali-Davarani, L. G. Raguin, D. A. Vorp, and S. Baek, “Identification of in vivo
material and geometric parameters of a human aorta: toward patient-specific modeling of
abdominal aortic aneurysm,” Biomechanics and Modeling in Mechanobiology, vol. 10, no.
5, pp. 689–699, Nov. 2010.
[61] MATLAB. MathWorks, 2010.
[62] S.-W. Lee, L. Antiga, J. D. Spence, and D. A. Steinman, “Geometry of the Carotid
Bifurcation Predicts Its Exposure to Disturbed Flow,” Stroke, vol. 39, no. 8, pp. 2341–2347,
Jun. 2008.
[63] J. B. Thomas, “Variation in the Carotid Bifurcation Geometry of Young Versus Older
Adults: Implications for Geometric Risk of Atherosclerosis,” Stroke, vol. 36, no. 11, pp.
2450–2456, Oct. 2005.
[64] D. A. Steinman, J. B. Thomas, H. M. Ladak, J. S. Milner, B. K. Rutt, and J. D. Spence,
“Reconstruction of carotid bifurcation hemodynamics and wall thickness using
computational fluid dynamics and MRI,” Magnetic Resonance in Medicine, vol. 47, no. 1,
pp. 149–159, Jan. 2002.

164